Count only reads that "aligned concordantly exactly 1 time."

        :param str aligned_bam: File for which to count mapped reads.
199241
def count_concordant(self, aligned_bam):
        
        cmd = self.tools.samtools + " view " + aligned_bam + " | "
        cmd += "grep 'YT:Z:CP'" + " | uniq -u | wc -l | sed -E 's/^[[:space:]]+//'"
        
        return subprocess.check_output(cmd, shell=True)
863096
def _reads_in_peaks(bam_file, peaks_file, sample):
    """"""
    if not peaks_file:
        return {}
    rip = number_of_mapped_reads(sample, bam_file, bed_file = peaks_file)
    return {"metrics": {"RiP": rip}}
440153
def parse_hisat2_logs(self, f):
        

        # Regexes
        regexes = {
            'unpaired_total': r"Total(?: unpaired)? reads: (\d+)",
            'unpaired_aligned_none': r"Aligned 0 times?: (\d+) \([\d\.]+%\)",
            'unpaired_aligned_one': r"Aligned 1 time: (\d+) \([\d\.]+%\)",
            'unpaired_aligned_multi': r"Aligned >1 times: (\d+) \([\d\.]+%\)",
    		'paired_total': r"Total pairs: (\d+)",
    		'paired_aligned_none': r"Aligned concordantly or discordantly 0 time: (\d+) \([\d\.]+%\)",
    		'paired_aligned_one': r"Aligned concordantly 1 time: (\d+) \([\d\.]+%\)",
    		'paired_aligned_multi': r"Aligned concordantly >1 times: (\d+) \([\d\.]+%\)",
    		'paired_aligned_discord_one': r"Aligned discordantly 1 time: (\d+) \([\d\.]+%\)",
        }

        # Go through log file line by line
        s_name = f['s_name']
        parsed_data = {}

        for l in f['f']:

            # Attempt in vain to find original hisat2 command, logged by another program
            hscmd = re.search(r"hisat2 .+ -[1U] ([^\s,]+)", l)
            if hscmd:
                s_name = self.clean_s_name(hscmd.group(1), f['root'])
                log.debug("Found a HISAT2 command, updating sample name to '{}'".format(s_name))

            # Run through all regexes
            for k, r in regexes.items():
                match = re.search(r, l)
                if match:
                    parsed_data[k] = int(match.group(1))

            # Overall alignment rate
            overall = re.search(r"Overall alignment rate: ([\d\.]+)%", l)
            if overall:
                parsed_data['overall_alignment_rate'] = float(overall.group(1))

                # Save parsed data
                if s_name in self.hisat2_data:
                    log.debug("Duplicate sample name found! Overwriting: {}".format(s_name))
                self.add_data_source(f, s_name)
                self.hisat2_data[s_name] = parsed_data

                # Reset in case we find more in this log file
                s_name = f['s_name']
                parsed_data = {}
1191583
def detect_alignment_annotations(queryalignment, tags=False):
    
    annotations = set()
    for k, v in BARCODEINFO.items():
        if tags:
            if queryalignment.has_tag(v.bamtag):
                annotations.add(k)
        else:
            if v.readprefix in queryalignment.qname:
                annotations.add(k)
    return annotations
187976
def process_alignment(data, alt_input=None):
    
    data = cwlutils.normalize_missing(utils.to_single_data(data))
    data = cwlutils.unpack_tarballs(data, data)
    fastq1, fastq2 = dd.get_input_sequence_files(data)
    if alt_input:
        fastq1, fastq2 = alt_input
    config = data["config"]
    aligner = config["algorithm"].get("aligner", None)
    if fastq1 and objectstore.file_exists_or_remote(fastq1) and aligner:
        logger.info("Aligning lane %s with %s aligner" % (data["rgnames"]["lane"], aligner))
        data = align_to_sort_bam(fastq1, fastq2, aligner, data)
        if dd.get_umi_consensus(data):
            data["umi_bam"] = dd.get_work_bam(data)
            if fastq2:
                f1, f2, avg_cov = postalign.umi_consensus(data)
                data["config"]["algorithm"]["rawumi_avg_cov"] = avg_cov
                del data["config"]["algorithm"]["umi_type"]
                data["config"]["algorithm"]["mark_duplicates"] = False
                data = align_to_sort_bam(f1, f2, aligner, data)
            else:
                raise ValueError("Single fastq input for UMI processing; fgbio needs paired reads: %s" %
                                 dd.get_sample_name(data))
        data = _add_supplemental_bams(data)
    elif fastq1 and objectstore.file_exists_or_remote(fastq1) and fastq1.endswith(".bam"):
        sort_method = config["algorithm"].get("bam_sort")
        bamclean = config["algorithm"].get("bam_clean")
        if bamclean is True or bamclean == "picard":
            if sort_method and sort_method != "coordinate":
                raise ValueError("Cannot specify `bam_clean: picard` with `bam_sort` other than coordinate: %s"
                                 % sort_method)
            out_bam = cleanbam.picard_prep(fastq1, data["rgnames"], dd.get_ref_file(data), data["dirs"],
                                           data)
        elif bamclean == "fixrg":
            out_bam = cleanbam.fixrg(fastq1, data["rgnames"], dd.get_ref_file(data), data["dirs"], data)
        elif bamclean == "remove_extracontigs":
            out_bam = cleanbam.remove_extracontigs(fastq1, data)
        elif sort_method:
            runner = broad.runner_from_path("picard", config)
            out_file = os.path.join(data["dirs"]["work"], "{}-sort.bam".format(
                os.path.splitext(os.path.basename(fastq1))[0]))
            if not utils.file_exists(out_file):
                work_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "bamclean",
                                                           dd.get_sample_name(data)))
                out_file = os.path.join(work_dir, "{}-sort.bam".format(dd.get_sample_name(data)))
            out_bam = runner.run_fn("picard_sort", fastq1, sort_method, out_file)
        else:
            out_bam = _link_bam_file(fastq1, os.path.join(dd.get_work_dir(data), "prealign",
                                                          dd.get_sample_name(data)), data)
        bam.index(out_bam, data["config"])
        bam.check_header(out_bam, data["rgnames"], dd.get_ref_file(data), data["config"])
        dedup_bam = postalign.dedup_bam(out_bam, data)
        bam.index(dedup_bam, data["config"])
        data["work_bam"] = dedup_bam
    elif fastq1 and objectstore.file_exists_or_remote(fastq1) and fastq1.endswith(".cram"):
        data["work_bam"] = fastq1
    elif fastq1 is None and not dd.get_aligner(data):
        data["config"]["algorithm"]["variantcaller"] = False
        data["work_bam"] = None
    elif not fastq1:
        raise ValueError("No 'files' specified for input sample: %s" % dd.get_sample_name(data))
    elif "kraken" in config["algorithm"]:  # kraken doesn's need bam
        pass
    else:
        raise ValueError("Could not process input file from sample configuration. \n" +
                         fastq1 +
                         "\nIs the path to the file correct or is empty?\n" +
                         "If it is a fastq file (not pre-aligned BAM or CRAM), "
                         "is an aligner specified in the input configuration?")
    if data.get("work_bam"):
        # Add stable 'align_bam' target to use for retrieving raw alignment
        data["align_bam"] = data["work_bam"]
        data = _add_hla_files(data)
    return [[data]]
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)
637636
def _check_reads_hit(self, alignment_io, min_aligned_fraction):
        ''''''
        to_return = []
        alignment_length = None
        for s in SeqIO.parse(alignment_io, "fasta"):
            if not alignment_length:
                alignment_length = len(s.seq)
                min_length = int(min_aligned_fraction * alignment_length)
                logging.debug("Determined min number of aligned bases to be %s" % min_length)
            elif len(s.seq) != alignment_length:
                raise Exception("Alignment file appears to not be of uniform length")

            num_unaligned = s.seq.count('-')
            num_aligned = alignment_length-num_unaligned
            logging.debug("Sequence %s has %d aligned positions" % (s.name, alignment_length-num_unaligned))
            if num_aligned <= min_length:
                to_return.append(s.name)
        return to_return
892143
def _extract_split_and_discordants(in_bam, work_dir, data):
    
    sr_file = os.path.join(work_dir, "%s-sr.bam" % os.path.splitext(os.path.basename(in_bam))[0])
    disc_file = os.path.join(work_dir, "%s-disc.bam" % os.path.splitext(os.path.basename(in_bam))[0])
    if not utils.file_exists(sr_file) or not utils.file_exists(disc_file):
        with file_transaction(data, sr_file) as tx_sr_file:
            with file_transaction(data, disc_file) as tx_disc_file:
                cores = dd.get_num_cores(data)
                ref_file = dd.get_ref_file(data)
                cmd = ("extract-sv-reads -e --threads {cores} -T {ref_file} "
                       "-i {in_bam} -s {tx_sr_file} -d {tx_disc_file}")
                do.run(cmd.format(**locals()), "extract split and discordant reads", data)
    for fname in [sr_file, disc_file]:
        bam.index(fname, data["config"])
    return sr_file, disc_file
1337233
def _get_fastq_in_region(region, align_bam, out_base):
    
    out_file = "{0}.fastq".format(out_base)
    if not file_exists(out_file):
        with pysam.Samfile(align_bam, "rb") as in_pysam:
            with file_transaction(out_file) as tx_out_file:
                with open(tx_out_file, "w") as out_handle:
                    contig, start, end = region
                    for read in in_pysam.fetch(contig, int(start), int(end)):
                        seq = Seq.Seq(read.seq)
                        qual = list(read.qual)
                        if read.is_reverse:
                            seq = seq.reverse_complement()
                            qual.reverse()
                        out_handle.write("@{name}\n{seq}\n+\n{qual}\n".format(
                                name=read.qname, seq=str(seq), qual="".join(qual)))
    return out_file
522195
def _counts_at_position(positions, orig_reader, cmp_reader):
    
    pos_counts = collections.defaultdict(lambda:
                 collections.defaultdict(lambda:
                 collections.defaultdict(int)))
    for orig_parts in orig_reader:
        cmp_parts = next(cmp_reader)
        for pos in positions:
            try:
                pos_counts[pos][int(orig_parts[pos+1])][int(cmp_parts[pos+1])] += 1
            except IndexError:
                pass
    for pos, count_dict in pos_counts.iteritems():
        for orig_val, cmp_dict in count_dict.iteritems():
            for cmp_val, count in cmp_dict.iteritems():
                yield pos+1, orig_val, cmp_val, count
429781
def parse_align_file(file_in):
    
    loc_id = 1
    bedfile_clusters = ""
    bamfile = pybedtools.BedTool(file_in)
    bed = pybedtools.BedTool.bam_to_bed(bamfile)
    for c, start, end, name, q, strand in bed:
        loc_id += 1
        bedfile_clusters += "%s\t%s\t%s\t%s\t%s\t%s\n" % \
                            (c, start, end, name, loc_id, strand)
    return bedfile_clusters
636978
def _cluster(bam_file, ma_file, out_dir, reference, annotation_file=None):
    
    seqcluster = op.join(get_bcbio_bin(), "seqcluster")
    # cl = ["cluster", "-o", out_dir, "-m", ma_file, "-a", bam_file, "-r", reference]
    if annotation_file:
        annotation_file = "-g " + annotation_file
    else:
        annotation_file = ""

    if not file_exists(op.join(out_dir, "counts.tsv")):
        cmd = ("{seqcluster} cluster -o {out_dir} -m {ma_file} -a {bam_file} -r {reference} {annotation_file}")
        do.run(cmd.format(**locals()), "Running seqcluster.")
    counts = op.join(out_dir, "counts.tsv")
    stats = op.join(out_dir, "read_stats.tsv")
    json = op.join(out_dir, "seqcluster.json")
    return {'out_dir': out_dir, 'count_file': counts, 'stat_file': stats, 'json': json}
1149178
def build_bam_tags():
    ''''''
    #pylint: disable=unused-argument
    def _combine_filters(fam, paired_align, align):
        filters = [x.filter_value for x in [fam, align] if x and x.filter_value]
        if filters:
            return ";".join(filters).replace('; ', ';')
        return None

    boolean_tag_value = {True:1}
    tags = [
        BamTag("X0", "Z",
               ("filter (why the alignment was excluded)"),
               _combine_filters),
        BamTag("X1", "Z",
               ("leftmost~rightmost matched pair positions"),
               lambda fam, pair, align: pair.positions('{left}~{right}')),
        BamTag("X2", "Z",
               ("L~R CIGARs"),
               lambda fam, pair, align: pair.cigars('{left}~{right}')),
        BamTag("X3", "i",
               "unique identifier for this alignment family",
               lambda fam, pair, align: fam.umi_sequence),
        BamTag("X4", "Z",
               ("L~R UMT barcodes for this alignment family; because "
                "of fuzzy matching the family UMT may be distinct "
                "from the UMT of the original alignment"),
               lambda fam, pair, align: fam.umt('{left}~{right}')),
        BamTag("X5", "i",
               "family size (number of align pairs in this family)",
               lambda fam, pair, align: fam.included_pair_count),
        BamTag("X6", "i",
               ("presence of this tag signals that this alignment "
                "would be the template for the consensus alignment"),
               lambda fam, pair, align: boolean_tag_value.get(fam.is_consensus_template(align), None))]
    return tags
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)
290402
def run_align(*data):
    
    work_dir = dd.get_work_dir(data[0][0])
    out_dir = op.join(work_dir, "seqcluster", "prepare")
    seq_out = op.join(out_dir, "seqs.fastq")
    bam_dir = op.join(work_dir, "align")
    new_bam_file = op.join(bam_dir, "seqs.bam")
    tools = dd.get_expression_caller(data[0][0])
    if not file_exists(new_bam_file):
        sample = process_alignment(data[0][0], [seq_out, None])
        bam_file = dd.get_work_bam(sample[0][0])
        shutil.move(bam_file, new_bam_file)
        shutil.move(bam_file + ".bai", new_bam_file + ".bai")
        shutil.rmtree(op.join(bam_dir, sample[0][0]["rgnames"]['sample']))
    for sample in data:
        # sample[0]["align_bam"] = sample[0]["clean_fastq"]
        sample[0]["cluster_bam"] = new_bam_file

    if "mirdeep2" in tools:
        novel_db = mirdeep.run(data)
    return data
1184277
def get_seqstarts(bamfile, N):
    
    import pysam
    bamfile = pysam.AlignmentFile(bamfile, "rb")
    seqsize = {}
    for kv in bamfile.header["SQ"]:
        if kv["LN"] < 10 * N:
            continue
        seqsize[kv["SN"]] = kv["LN"] / N + 1

    allseqs = natsorted(seqsize.keys())
    allseqsizes = np.array([seqsize[x] for x in allseqs])
    seqstarts = np.cumsum(allseqsizes)
    seqstarts = np.roll(seqstarts, 1)
    total_bins = seqstarts[0]
    seqstarts[0] = 0
    seqstarts = dict(zip(allseqs, seqstarts))

    return seqstarts, seqsize, total_bins
623274
def clean_bam_file(bam_in, mask=None):
    
    seq_obj = defaultdict(int)
    if mask:
        mask_file = op.splitext(bam_in)[0] + "_mask.bam"
        if not file_exists(mask_file):
            pybedtools.BedTool(bam_file).intersect(b=mask, v=True).saveas(mask_file)
        bam_in = mask_file
    out_file = op.splitext(bam_in)[0] + "_rmlw.bam"
    # bam.index(bam_in, {'algorithm':{}})
    run("samtools index %s" % bam_in)
    if not file_exists(bam_in + ".bai"):
        raise IOError("Failed to created bam index of %s. Try to do it manually" % bam_in)
    bam_handle = pysam.AlignmentFile(bam_in, "rb")
    with pysam.AlignmentFile(out_file, "wb", template=bam_handle) as out_handle:
        for read in bam_handle.fetch():
            seq_name = int(read.query_name.replace('seq_', ''))
            match_size = [nts for oper, nts in read.cigartuples if oper == 0]
            subs_size = [nts for oper, nts in read.cigartuples if oper == 4]
            if match_size[0] < 17:
                continue
            if subs_size:
                if subs_size[0] > 3:
                    continue
            try:
                nh = read.get_tag('NH')
            except KeyError:
                nh = 1
            seq_obj[seq_name] = sequence(seq_name)
            seq_obj[seq_name].align = nh
            out_handle.write(read)
    return out_file, seq_obj