Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def _process_reads(bam_fp, start_stop, all_alignments=False, min_length=None, bed_file=None):
start, stop = start_stop
counts = Counter()
results = []
with pysam.AlignmentFile(bam_fp, 'rb') as bam:
if bed_file is not None:
trees = intervaltrees_from_bed(bed_file)
for i, read in enumerate(bam.fetch(until_eof=True)):
if i < start:
continue
elif i == stop:
break
if read.is_unmapped:
counts['unmapped'] += 1
continue
if not all_alignments and (read.is_secondary or read.is_supplementary):
continue
if bed_file is not None:
if not trees[read.reference_name].overlaps(read.reference_start, read.reference_end):
sys.stderr.write('read {} does not overlap with any regions in bedfile\n'.format(read.query_name))
counts['masked'] += 1
continue
def _process_read(bam, read_num, bed_file=None):
"""Load an alignment from bam and return result of `_process_seg`.
:param bam: str, bam file.
:param read_num: int, index of alignment to process.
:param bed_file: path to .bed file of regions to include in analysis.
:returns: result of `_process_seg`.
"""
trees = None
if bed_file is not None:
trees = intervaltrees_from_bed(bed_file)
with pysam.AlignmentFile(bam, 'rb') as bam_obj:
gen = (r for r in bam_obj)
for i in range(read_num + 1):
rec = next(gen)
if rec.is_unmapped or rec.is_supplementary or rec.is_secondary:
return
if bed_file is not None:
tree = trees[rec.reference_name]
if not tree.overlaps(rec.reference_start, rec.reference_end):
#sys.stderr.write('read {} does not overlap with any regions in bedfile\n'.format(rec.query_name))
return
else:
tree = None