#!/usr/bin/env python3 # -*- coding: UTF-8 -*- import sys import pysam # 提取fr+rr>=2的点 if len(sys.argv) != 3: print(" ".join(['usage:python', sys.argv[0], 'outputDir', 'tumor'])) sys.exit() # output_dir = sys.argv[1] # tumor = sys.argv[2] # infile = "".join([output_dir, '/alignment/', tumor, '.rmdup.bam']) # outfile = "".join([output_dir, '/alignment/', tumor, '.2r.rmdup.bam']) samfile = pysam.AlignmentFile(sys.argv[1], "rb") bamfile = pysam.AlignmentFile(sys.argv[2], "wb", template=samfile) for s in samfile: fr = 0 rr = 0 if s.has_tag('FR'): fr = s.get_tag('FR') if s.has_tag('RR'): rr = s.get_tag('RR') if fr + rr >= 2: bamfile.write(s) # os.system('samtools index ' + outfile)