Traverse a BAM file to find all the reads overlapping a specified locus.

    Parameters are the same as those for read_locus_generator.
335389
def locus_reads_dataframe(*args, **kwargs):
    
    df_builder = DataFrameBuilder(
        LocusRead,
        variant_columns=False,
        converters={
            "reference_positions": list_to_string,
            "quality_scores": list_to_string,
        })
    for locus_read in locus_read_generator(*args, **kwargs):
        df_builder.add(variant=None, element=locus_read)
    return df_builder.to_dataframe()
151076
def alts(args):
    
    import pysam
    from jcvi.utils.iter import pairwise
    from jcvi.utils.grouper import Grouper

    p = OptionParser(alts.__doc__)
    p.set_outfile(outfile="TREDs.alts.csv")
    opts, args = p.parse_args(args)

    if len(args) < 1:
        sys.exit(not p.print_help())

    treds = args
    repo = TREDsRepo()
    if "all" in treds:
        treds = repo.names

    pad_left, pad_right = 1000, 10000
    READLEN = 150
    fw = must_open(opts.outfile, "w")
    print("TRED,alts,alts.hg19", file=fw)  # Header
    for tred in treds:
        ref_regions = []

        # Simulate a depth 1000 BAM with 300 repeats
        for ref in ("hg38", "hg19"):

            # This is the region that involves the TRED locus
            repo = TREDsRepo(ref=ref)
            t = repo[tred]
            chr, start, end = t.chr, t.repeat_start, t.repeat_end
            start -= pad_left
            end += pad_right

            tred_ref = "{}_{}".format(tred, ref)
            if not op.isdir(tred_ref):
                simulate([tred_ref, "300", "300", "--depth=1000",
                          "--ref={}".format(ref), "--tred={}".format(tred)])
            bamfile = op.join(tred_ref, "300.bam")

            # Parse the BAM file, retrieve all regions
            bamfile = pysam.AlignmentFile(bamfile, "rb")
            nreads = altreads = 0
            alt_points = set()
            for read in bamfile.fetch():
                fname, fstart = bamfile.getrname(read.reference_id), read.reference_start
                rname, rstart = bamfile.getrname(read.next_reference_id), read.next_reference_start
                f_in_region = in_region(fname, fstart, chr, start, end)
                r_in_region = in_region(rname, rstart, chr, start, end)
                if (not f_in_region) and r_in_region:
                    alt_points.add((fname, fstart))
                    altreads += 1
                if (not r_in_region) and f_in_region:
                    alt_points.add((rname, rstart))
                    altreads += 1
                nreads += 1

            logging.debug("A total of {} reads ({} alts) processed".\
                        format(nreads, altreads))
            alt_points = natsorted(alt_points)

            # Chain these points together into regions
            g = Grouper()
            for a in alt_points:
                g.join(a)
            for a, b in pairwise(alt_points):
                achr, apos = a
                bchr, bpos = b
                if achr != bchr:
                    continue
                if (bpos - apos) > READLEN:
                    continue
                g.join(a, b)

            # All regions that contain ALT
            alt_sum = 0
            regions = []
            for c in g:
                chr_min, pos_min = min(c)
                chr_max, pos_max = max(c)
                assert chr_min, chr_max
                pos_min -= READLEN
                pos_max += READLEN
                regions.append((chr_min, pos_min, pos_max))
                alt_sum += pos_max - pos_min

            regions = "|".join(["{}:{}-{}".format(c, start, end) \
                            for c, start, end in natsorted(regions)])
            ref_regions.append(regions)

        line = ",".join([tred] + ref_regions)
        print(line, file=sys.stderr)
        print(line, file=fw)
        logging.debug("Alternative region sum: {} bp".format(alt_sum))

    fw.close()
876452
def find_features(seqs, locus_tag="all", utr_len=200):
    """"""
    found_features = []

    for seq_i in seqs:
        for feature in seq_i.features:
            if feature.type == "CDS" and (locus_tag == "all" or \
                    ('locus_tag' in feature.qualifiers and \
                    feature.qualifiers['locus_tag'][0] == locus_tag)):
                start = max(0, feature.location.nofuzzy_start - utr_len)
                stop  = max(0, feature.location.nofuzzy_end + utr_len)
                feature_seq = seq_i.seq[start:stop]
                f_match = FeatureMatch(feature, feature_seq, feature.strand,
                                       utr_len)
                found_features.append(f_match)

    return found_features
1213820
def greedy_merge_helper(
        variant_sequences,
        min_overlap_size=MIN_VARIANT_SEQUENCE_ASSEMBLY_OVERLAP_SIZE):
    
    merged_variant_sequences = {}
    merged_any = False

    # here we'll keep track of sequences that haven't been merged yet, and add them in at the end
    unmerged_variant_sequences = set(variant_sequences)
    for i in range(len(variant_sequences)):
        sequence1 = variant_sequences[i]
        # it works to loop over the triangle (i+1 onwards) because combine() tries flipping the
        # arguments if sequence1 is on the right of sequence2
        for j in range(i + 1, len(variant_sequences)):
            sequence2 = variant_sequences[j]
            combined = sequence1.combine(sequence2)
            if combined is None:
                continue
            if combined.sequence in merged_variant_sequences:
                existing = merged_variant_sequences[combined.sequence]
                # the existing VariantSequence and the newly merged
                # VariantSequence should differ only in which reads support them
                combined = combined.add_reads(existing.reads)
            merged_variant_sequences[combined.sequence] = combined
            unmerged_variant_sequences.discard(sequence1)
            unmerged_variant_sequences.discard(sequence2)
            merged_any = True
    result = list(merged_variant_sequences.values()) + list(unmerged_variant_sequences)
    return result, merged_any
831755
def overlaps(self,in_genomic_range,padding=0):
    
    if padding > 0:
      in_genomic_range = GenomicRange(in_genomic_range.chr,max([1,in_genomic_range.start-padding]),in_genomic_range.end+padding)
    if self.chr != in_genomic_range.chr:
      return False
    if self.end < in_genomic_range.start:
      return False
    if in_genomic_range.end < self.start:
      return False
    if self.start > in_genomic_range.end:
      return False
    if in_genomic_range.start > self.end:
      return False
    if self.start <= in_genomic_range.start and self.end >= in_genomic_range.start:
      return True
    if self.start <= in_genomic_range.end and self.end >= in_genomic_range.end:
      return True
    if self.start >= in_genomic_range.start and self.end <= in_genomic_range.end:
      return True
    if self.start <= in_genomic_range.start and self.end >= in_genomic_range.end:
      return True
    if in_genomic_range.start <= self.start and in_genomic_range.end >= self.start:
      return True
    if in_genomic_range.start <= self.end and in_genomic_range.end >= self.end:
      return True
    sys.stderr.write("overlaps: unprogrammed error\n")
    return False
975020
def write_map(self, path):
        '''

        '''
        starts = [index[0] for index in self.overlap_indices]
        features = []
        for i, start in enumerate(starts):
            stop = start + len(self.overlaps[i])
            name = 'overlap {}'.format(i + 1)
            feature_type = 'misc'
            strand = 0
            features.append(coral.Feature(name, start, stop, feature_type,
                                          strand=strand))
        seq_map = coral.DNA(self.template, features=features)
        coral.seqio.write_dna(seq_map, path)
724786
def paired_reader_from_bamfile(args,
                               log,
                               usage_logger,
                               annotated_writer):
    ''''''
    total_aligns = pysamwrapper.total_align_count(args.input_bam)
    bamfile_generator  = _bamfile_generator(args.input_bam)
    return _paired_reader(args.umt_length,
                          bamfile_generator,
                          total_aligns,
                          log,
                          usage_logger,
                          annotated_writer)
906323
def trim_by_coverage(self, min_reads):
        
        read_count_array = self.coverage()
        logger.info("Coverage: %s (len=%d)" % (
            read_count_array, len(read_count_array)))
        sufficient_coverage_mask = read_count_array >= min_reads
        sufficient_coverage_indices = np.argwhere(sufficient_coverage_mask)
        if len(sufficient_coverage_indices) == 0:
            logger.debug("No bases in %s have coverage >= %d" % (self, min_reads))
            return VariantSequence(prefix="", alt="", suffix="", reads=self.reads)
        variant_start_index, variant_end_index = self.variant_indices()
        # assuming that coverage drops off monotonically away from
        # variant nucleotides
        first_covered_index = sufficient_coverage_indices.min()
        last_covered_index = sufficient_coverage_indices.max()
        # adding 1 to last_covered_index since it's an inclusive index
        # whereas variant_end_index is the end of a half-open interval
        if (first_covered_index > variant_start_index or
                last_covered_index + 1 < variant_end_index):
            # Example:
            #   Nucleotide sequence:
            #       ACCCTTTT|AA|GGCGCGCC
            #   Coverage:
            #       12222333|44|33333211
            # Then the mask for bases covered >= 4x would be:
            #       ________|**|________
            # with indices:
            #       first_covered_index = 9
            #       last_covered_index = 10
            #       variant_start_index = 9
            #       variant_end_index = 11
            logger.debug("Some variant bases in %s don't have coverage >= %d" % (
                self, min_reads))
            return VariantSequence(prefix="", alt="", suffix="", reads=self.reads)
        return VariantSequence(
            prefix=self.prefix[first_covered_index:],
            alt=self.alt,
            suffix=self.suffix[:last_covered_index - variant_end_index + 1],
            reads=self.reads)
231986
def epcrparse(self):
        
        from Bio.Blast.Applications import NcbiblastnCommandline
        while True:
            sample, record, line = self.epcrparsequeue.get()
            # Split the data on tabs
            gene, chromosome, strand, start, end, m_match, gaps, act_len_exp_len = line.split('\t')
            # Extract the gene sequence from the contigs
            # The record dictionary has the contig name, and the sequence. Splice out the data using the start and
            # end coordinates specified by ePCR
            genesequence = record[chromosome][int(start) - 1:int(end)]
            # Set up BLASTn using blastn-short, as the probe sequences tend to be very short
            blastn = NcbiblastnCommandline(db=sample[self.analysistype].probes.split('.')[0],
                                           num_threads=12,
                                           task='blastn-short',
                                           num_alignments=1,
                                           outfmt="'6 qseqid sseqid positive mismatch gaps "
                                                  "evalue bitscore slen length'")
            # Run the BLASTn, with the gene sequence as stdin
            out, err = blastn(stdin=genesequence)
            # Split the output string on tabs
            results = out.rstrip().split('\t')
            # Populate the raw blast results
            sample[self.analysistype].rawblastresults[gene] = results
            # Create named variables from the list
            positives = float(results[2])
            mismatches = float(results[3])
            gaps = float(results[4])
            subjectlength = float(results[7])
            # Calculate the percent identity
            percentidentity = float('{:0.2f}'.format((positives - gaps) / subjectlength * 100))
            # Create a dictionary with the desired values to store in the metadata object
            resultdict = {
                'matches': positives,
                'mismatches': mismatches,
                'gaps': gaps,
                'subject_length': subjectlength,
                'percent_identity': percentidentity,
                'match_length': results[8].split('\n')[0]
            }
            # Populate the metadata object with the dictionary
            sample[self.analysistype].blastresults[gene] = resultdict
            self.epcrparsequeue.task_done()
834895
def _overlapping(files):
    
    segments = set()
    for path in files:
        seg = file_segment(path)
        for s in segments:
            if seg.intersects(s):
                return True
        segments.add(seg)
    return False
363088
def _read_geneset(self, study_fn, pop_fn):
        """"""
        pop = set(_.strip() for _ in open(pop_fn) if _.strip())
        study = frozenset(_.strip() for _ in open(study_fn) if _.strip())
        if next(iter(pop)).isdigit():
            pop = set(int(g) for g in pop)
            study = frozenset(int(g) for g in study)
        # some times the pop is a second group to compare, rather than the
        # population in that case, we need to make sure the overlapping terms
        # are removed first
        if self.args.compare:
            common = pop & study
            pop |= study
            pop -= common
            study -= common
            sys.stderr.write("removed %d overlapping items\n" % (len(common)))
            sys.stderr.write("Set 1: {0}, Set 2: {1}\n".format(
                len(study), len(pop)))
        return study, pop
790320
def subset_bamfile(sam, barcodes):
    
    from pysam import AlignmentFile

    start_time = time.time()

    sam_file = open_bamfile(sam)
    out_file = AlignmentFile("-", "wh", template=sam_file)
    track = sam_file.fetch(until_eof=True)

    # peek at first alignment to determine the annotations
    queryalignment = track.next()
    annotations = detect_alignment_annotations(queryalignment)
    track = itertools.chain([queryalignment], track)

    re_string = construct_transformed_regex(annotations)
    parser_re = re.compile(re_string)
    barcodes = set(barcode.strip() for barcode in barcodes)

    for count, aln in enumerate(track, start=1):
        if count and not count % 1000000:
            logger.info("Processed %d alignments." % count)

        match = parser_re.match(aln.qname)
        tags = aln.tags

        if "cellular" in annotations:
            cb = match.group('CB')
            if cb in barcodes:
                out_file.write(aln)
1119661
def _create_clusters(seqL, bam_file, args):
    
    clus_obj = []
    cluster_file = op.join(args.out, "cluster.bed")
    if not os.path.exists(op.join(args.out, 'list_obj.pk')):
        if not file_exists(cluster_file):
            logger.info("Parsing aligned file")
            logger.info("Merging sequences")
            bedtools = os.path.join(os.path.dirname(sys.executable), "bedtools")
            bedtools = bedtools if os.path.exists(bedtools) else "bedtools"
            parse_cmd = "awk '{i=i+1;print $1\"\\t\"$2\"\\t\"$3\"\\t\"$4\"\\t\"i\"\\t\"$6}'"
            cmd = "{bedtools} bamtobed -i {bam_file} | {parse_cmd} | {bedtools} cluster -s -d 20 -i - > {cluster_file}"
            do.run(cmd.format(**locals()))
        c = pybedtools.BedTool(cluster_file)
        logger.info("Creating clusters")
        clus_obj = detect_clusters(c, seqL, args.min_seqs, args.non_un_gl)
        with open(op.join(args.out, 'list_obj.pk'), 'wb') as output:
            pickle.dump(clus_obj, output, pickle.HIGHEST_PROTOCOL)
    else:
        logger.info("Loading previous clusters")
        with open(op.join(args.out, 'list_obj.pk'), 'rb') as input:
            clus_obj = pickle.load(input)
    # bedfile = pybedtools.BedTool(generate_position_bed(clus_obj), from_string=True)
    # seqs_2_loci = bedfile.intersect(pybedtools.BedTool(aligned_bed, from_string=True), wo=True, s=True)
    # seqs_2_position = add_seqs_position_to_loci(seqs_2_loci, seqL)
    logger.info("%s clusters found" % (len(clus_obj.clusid)))
    return clus_obj
1327429
def get_variants_in_region(self, chrom, start, end):
        """"""
        if self.chrom is not None and chrom == self.chrom:
            # We are going to search for 'NA' since the chromosome was set
            chrom = "NA"

        iterator = self._bgen.iter_variants_in_region(
            CHROM_STR_DECODE.get(chrom, chrom), start, end,
        )
        for info, dosage in iterator:
            yield Genotypes(
                Variant(
                    info.name, CHROM_STR_ENCODE.get(info.chrom, info.chrom),
                    info.pos, [info.a1, info.a2],
                ),
                dosage,
                reference=info.a1,
                coded=info.a2,
                multiallelic=True,
            )
1026404
def _find_gene_list_from_bed(bed_file, base_file, data):
    
    # Check for a gene list, we can just return that.
    if is_gene_list(bed_file):
        return bed_file
    out_file = "%s-genes.txt" % utils.splitext_plus(base_file)[0]
    if not os.path.exists(out_file):
        genes = set([])
        import pybedtools
        with utils.open_gzipsafe(bed_file) as in_handle:
            for r in pybedtools.BedTool(in_handle):
                if r.name:
                    if not r.name.startswith("{"):
                        genes.add(r.name)
        with file_transaction(data, out_file) as tx_out_file:
            with open(tx_out_file, "w") as out_handle:
                if len(genes) > 0:
                    out_handle.write("\n".join(sorted(list(genes))) + "\n")
    if utils.file_exists(out_file):
        return out_file
330209
def _parse_gene_coordinates(file):
        
        id_map = {}
        col = ['gene_curie', 'start', 'end', 'strand', 'build']
        if os.path.exists(os.path.join(os.path.dirname(__file__), file)):
            with open(os.path.join(os.path.dirname(__file__), file)) as tsvfile:
                reader = csv.reader(tsvfile, delimiter="\t")
                for row in reader:
                    id_map[row[col.index('gene_curie')]] = {
                        'start': row[col.index('start')],
                        'end': row[col.index('end')],
                        'strand': row[col.index('strand')],
                        'build': row[col.index('build')]
                    }
        return id_map
723333
def synteny_scan(points, xdist, ydist, N):
    
    clusters = Grouper()
    n = len(points)
    points.sort()
    for i in range(n):
        for j in range(i - 1, -1, -1):
            # x-axis distance
            del_x = points[i][0] - points[j][0]
            if del_x > xdist:
                break
            # y-axis distance
            del_y = points[i][1] - points[j][1]
            if abs(del_y) > ydist:
                continue
            # otherwise join
            clusters.join(points[i], points[j])

    # select clusters that are at least >=N
    clusters = [sorted(cluster) for cluster in list(clusters) \
            if _score(cluster) >= N]

    return clusters