How to use the pomoxis.util.intervaltrees_from_bed function in pomoxis

To help you get started, we’ve selected a few pomoxis examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github nanoporetech / pomoxis / pomoxis / stats_from_bam.py View on Github external
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
github nanoporetech / pomoxis / pomoxis / catalogue_errors.py View on Github external
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