1212
1313from Bio import SeqIO
1414
15- from . import samtools , picard # was: from viral_ngs import tools
1615from . import samtools
1716from . import picard
1817from . import file as util_file , misc as util_misc # was: from viral_ngs import util
@@ -54,7 +53,7 @@ def align_bam(self, inBam, refDb, outBam, options=None,
5453 threads = util_misc .sanitize_thread_count (threads )
5554
5655 # fetch list of RGs
57- rgs = list (samtools .getReadGroups (inBam ).keys ())
56+ rgs = list (samtools_tool .getReadGroups (inBam ).keys ())
5857
5958 if len (rgs ) == 0 :
6059 # Can't do this
@@ -79,16 +78,16 @@ def align_bam(self, inBam, refDb, outBam, options=None,
7978 threads = threads ,
8079 should_index = False # Don't index intermediate BAMs that will be merged
8180 )
82- if not samtools .isEmpty (tmp_bam ):
81+ if not samtools_tool .isEmpty (tmp_bam ):
8382 align_bams .append (tmp_bam )
8483 else :
8584 log .warning ("No alignment output for RG %s in file %s against %s" , rg , inBam , refDb )
8685
8786 if len (align_bams )== 0 :
8887 log .warning ("All read groups in file %s appear to be empty." , inBam )
8988 with util_file .tempfname ('.empty.sam' ) as empty_sam :
90- samtools .dumpHeader (inBam , empty_sam )
91- samtools .sort (empty_sam , outBam )
89+ samtools_tool .dumpHeader (inBam , empty_sam )
90+ samtools_tool .sort (empty_sam , outBam )
9291 else :
9392 # Merge BAMs, sort, and index
9493 picardOptions = ['SORT_ORDER=coordinate' , 'USE_THREADING=true' , 'CREATE_INDEX=true' ]
@@ -117,7 +116,7 @@ def align_one_rg(self, inBam, refDb, outBam, rgid=None, preset=None, options=Non
117116 samtools_tool = samtools .SamtoolsTool ()
118117
119118 # Require exactly one RG
120- rgs = samtools .getReadGroups (inBam )
119+ rgs = samtools_tool .getReadGroups (inBam )
121120 if len (rgs ) == 0 :
122121 raise InvalidBamHeaderError ("{} lacks read groups" .format (inBam ))
123122 elif len (rgs ) == 1 :
@@ -133,13 +132,13 @@ def align_one_rg(self, inBam, refDb, outBam, rgid=None, preset=None, options=Non
133132 removeInput = False
134133 if len (rgs ) == 1 :
135134 one_rg_inBam = inBam
136- samtools . SamtoolsTool () .dumpHeader (one_rg_inBam , headerFile )
135+ samtools_tool .dumpHeader (one_rg_inBam , headerFile )
137136 else :
138137 # strip inBam to one read group
139138 with util_file .tempfname ('.onebam.bam' ) as tmp_bam :
140- samtools .view (['-1' , '-r' , rgid ], inBam , tmp_bam )
139+ samtools_tool .view (['-1' , '-r' , rgid ], inBam , tmp_bam )
141140 # special exit if this file is empty
142- if samtools .isEmpty (tmp_bam ):
141+ if samtools_tool .isEmpty (tmp_bam ):
143142 log .warning ("No reads present for RG %s in file: %s" , rgid , inBam )
144143 shutil .copyfile (tmp_bam , outBam )
145144 return
@@ -148,13 +147,13 @@ def align_one_rg(self, inBam, refDb, outBam, rgid=None, preset=None, options=Non
148147 removeInput = True
149148
150149 with open (headerFile , 'wt' ) as outf :
151- for row in samtools .getHeader (inBam ):
150+ for row in samtools_tool .getHeader (inBam ):
152151 if len (row ) > 0 and row [0 ] == '@RG' :
153152 if rgid != list (x [3 :] for x in row if x .startswith ('ID:' ))[0 ]:
154153 # skip all read groups that are not rgid
155154 continue
156155 outf .write ('\t ' .join (row ) + '\n ' )
157- samtools .reheader (tmp_bam , headerFile , one_rg_inBam )
156+ samtools_tool .reheader (tmp_bam , headerFile , one_rg_inBam )
158157
159158 # get the read group line to give to mm2
160159 readgroup_line = ""
@@ -187,7 +186,7 @@ def align_one_rg(self, inBam, refDb, outBam, rgid=None, preset=None, options=Non
187186 options .extend (('-x' , preset ))
188187
189188 # perform actual alignment
190- if samtools .isEmpty (one_rg_inBam ):
189+ if samtools_tool .isEmpty (one_rg_inBam ):
191190 log .warning ("Input file %s appears to lack reads for RG '%s'" , inBam , rgid )
192191 # minimap doesn't like empty inputs, so copy empty bam through
193192 # samtools.sort(one_rg_inBam, outBam)
@@ -213,20 +212,21 @@ def align_cmd(self, inReads, refDb, outAlign, options=None, threads=None, should
213212 samtools_tool = samtools .SamtoolsTool ()
214213
215214 with util_file .tempfname ('.aligned.sam' ) as aln_sam :
216- fastq_pipe = samtools .bam2fq_pipe (inReads )
215+ fastq_pipe = samtools_tool .bam2fq_pipe (inReads )
217216 options .extend (('-a' , refDb , '-' , '-o' , aln_sam ))
218217 self .execute (options , stdin = fastq_pipe .stdout )
219218 if fastq_pipe .wait ():
220219 raise subprocess .CalledProcessError (fastq_pipe .returncode , "samtools.bam2fq_pipe() for {}" .format (inReads ))
221- samtools .sort (aln_sam , outAlign , threads = threads )
220+ samtools_tool .sort (aln_sam , outAlign , threads = threads )
222221
223222 # cannot index sam files; only do so if a bam/cram is desired
224223 if should_index and (outAlign .endswith (".bam" ) or outAlign .endswith (".cram" )):
225- samtools .index (outAlign , threads = threads )
224+ samtools_tool .index (outAlign , threads = threads )
226225
227226 def scaffold (self , contigs_fasta , ref_fasta , outAlign , divergence = 20 , options = None , threads = None ):
228227 options = [] if not options else options
229228
229+ samtools_tool = samtools .SamtoolsTool ()
230230 threads = util_misc .sanitize_thread_count (threads )
231231 if '-t' not in options :
232232 options .extend (('-t' , str (threads )))
@@ -243,11 +243,11 @@ def scaffold(self, contigs_fasta, ref_fasta, outAlign, divergence=20, options=No
243243 with util_file .tempfname ('.aligned.sam' ) as aln_sam :
244244 options .extend (('-a' , ref_fasta , contigs_fasta , '-o' , aln_sam ))
245245 self .execute (options )
246- samtools .sort (aln_sam , outAlign , threads = threads )
246+ samtools_tool .sort (aln_sam , outAlign , threads = threads )
247247
248248 # cannot index sam files; only do so if a bam/cram is desired
249249 if (outAlign .endswith (".bam" ) or outAlign .endswith (".cram" )):
250- samtools .index (outAlign )
250+ samtools_tool .index (outAlign )
251251
252252 def idxstats (self , inReads , refDb , outIdxstats , outReadlist = None , threads = None ):
253253 """
@@ -302,7 +302,7 @@ def idxstats(self, inReads, refDb, outIdxstats, outReadlist=None, threads=None):
302302
303303 try :
304304 # Start samtools bam2fq pipe for input (use 4 threads for BAM decompression)
305- fastq_pipe = samtools .bam2fq_pipe (inReads , threads = 4 )
305+ fastq_pipe = samtools_tool .bam2fq_pipe (inReads , threads = 4 )
306306
307307 # Start minimap2 process with PAF output to stdout
308308 tool_cmd = [self .install_and_get_path ()] + options
0 commit comments