@@ -129,23 +129,16 @@ def get_bwt_sequences(self):
129129
130130 os .system ("samtools index {input_bam}" .format (
131131 input_bam = self .input_bam_file ))
132- aligners = []
133- bam_file = pysam .AlignmentFile (self .input_bam_file , "rb" )
134- header = bam_file .text .split ("\n " )
135-
136- for h in header :
137- if "@PG" in h :
138- fields = h .strip ().split ("\t " )
139- for f in fields :
140- if f .startswith ("ID:" ):
141- aligners .append (f [3 :])
142- # print(">>>>", aligners)
143- if "KMA" in aligners :
144- os .system ("""samtools view -F 4 -F 2048 {bam} | while read line; do awk '{cmd}'; done > {out}"""
145- .format (bam = self .input_bam_file , cmd = """{print ">"$1"__"$4"__"$3"__"$5"\\ n"$11}""" , out = self .fasta_file ))
146- else :
147- os .system ("""samtools view -F 4 -F 2048 {bam} | while read line; do awk '{cmd}'; done > {out}"""
148- .format (bam = self .input_bam_file , cmd = """{print ">"$1"__"$3"__"$2"__"$5"\\ n"$10}""" , out = self .fasta_file ))
132+
133+ exclude_flag = 4 | 2048
134+
135+ with open (os .path .join (self .fasta_file ), 'w' ) as fout :
136+ with pysam .AlignmentFile (self .input_bam_file , "rb" ) as bam :
137+ for read in bam .fetch (until_eof = True ):
138+ if read .flag & exclude_flag :
139+ continue
140+ fout .write ('>%s__%s__%s__%s\n %s\n ' % (
141+ read .query_name .replace (" " , "_" ), read .reference_name , read .flag , read .mapping_quality , read .seq ))
149142
150143 def get_bwt_alignment_data (self , header ):
151144 """
0 commit comments