Perform full BAM preparation using pipes to avoid intermediate disk IO.

    Handles realignment of original BAMs.
351511
def piped_bamprep(data, region=None, out_file=None):
    
    data["region"] = region
    if not _need_prep(data):
        return [data]
    else:
        utils.safe_makedir(os.path.dirname(out_file))
        if region[0] == "nochrom":
            prep_bam = shared.write_nochr_reads(data["work_bam"], out_file, data["config"])
        elif region[0] == "noanalysis":
            prep_bam = shared.write_noanalysis_reads(data["work_bam"], region[1], out_file,
                                                     data["config"])
        else:
            if not utils.file_exists(out_file):
                with tx_tmpdir(data) as tmp_dir:
                    _piped_bamprep_region(data, region, out_file, tmp_dir)
            prep_bam = out_file
        bam.index(prep_bam, data["config"])
        data["work_bam"] = prep_bam
        return [data]
1206064
def _finalize_merge(out_file, bam_files, config):
    
    # Ensure timestamps are up to date on output file and index
    # Works around issues on systems with inconsistent times
    for ext in ["", ".bai"]:
        if os.path.exists(out_file + ext):
            subprocess.check_call(["touch", out_file + ext])
    for b in bam_files:
        utils.save_diskspace(b, "BAM merged to %s" % out_file, config)
338452
def _run_wham(inputs, background_bams):
    
    out_file = os.path.join(_sv_workdir(inputs[0]), "%s-wham.vcf.gz" % dd.get_sample_name(inputs[0]))
    if not utils.file_exists(out_file):
        with file_transaction(inputs[0], out_file) as tx_out_file:
            cores = dd.get_cores(inputs[0])
            ref_file = dd.get_ref_file(inputs[0])
            include_chroms = ",".join([c.name for c in ref.file_contigs(ref_file)
                                       if chromhacks.is_autosomal_or_x(c.name)])
            all_bams = ",".join([x["align_bam"] for x in inputs] + background_bams)
            cmd = ("whamg -x {cores} -a {ref_file} -f {all_bams} -c {include_chroms} "
                   "| bgzip -c > {tx_out_file}")
            do.run(cmd.format(**locals()), "WHAM SV caller: %s" % ", ".join(dd.get_sample_name(d) for d in inputs))
    return vcfutils.bgzip_and_index(out_file, inputs[0]["config"])
173819
def prepend_to_file(path, data, bufsize=1<<15):
    
    # Backup the file #
    backupname = path + os.extsep + 'bak'
    # Remove previous backup if it exists #
    try: os.unlink(backupname)
    except OSError: pass
    os.rename(path, backupname)
    # Open input/output files,  note: outputfile's permissions lost #
    with open(backupname) as inputfile:
        with open(path, 'w') as outputfile:
            outputfile.write(data)
            buf = inputfile.read(bufsize)
            while buf:
                outputfile.write(buf)
                buf = inputfile.read(bufsize)
    # Remove backup on success #
    os.remove(backupname)
1063851
def fai_from_bam(ref_file, bam_file, out_file, data):
    
    contigs = set([x.contig for x in idxstats(bam_file, data)])
    if not utils.file_uptodate(out_file, bam_file):
        with open(ref.fasta_idx(ref_file, data["config"])) as in_handle:
            with file_transaction(data, out_file) as tx_out_file:
                with open(tx_out_file, "w") as out_handle:
                    for line in (l for l in in_handle if l.strip()):
                        if line.split()[0] in contigs:
                            out_handle.write(line)
    return out_file
1115619
def _mutect_call_prep(align_bams, items, ref_file, assoc_files,
                       region=None, out_file=None):
    
    base_config = items[0]["config"]
    broad_runner = broad.runner_from_path("picard", base_config)
    broad_runner.run_fn("picard_index_ref", ref_file)

    broad_runner = broad.runner_from_config(base_config, "mutect")
    _check_mutect_version(broad_runner)
    for x in align_bams:
        bam.index(x, base_config)

    paired = vcfutils.get_paired_bams(align_bams, items)
    if not paired:
        raise ValueError("Specified MuTect calling but 'tumor' phenotype not present in batch\n"
                         "https://bcbio-nextgen.readthedocs.org/en/latest/contents/"
                         "pipelines.html#cancer-variant-calling\n"
                         "for samples: %s" % ", " .join([dd.get_sample_name(x) for x in items]))
    params = ["-R", ref_file, "-T", "MuTect", "-U", "ALLOW_N_CIGAR_READS"]
    params += ["--read_filter", "NotPrimaryAlignment"]
    params += ["-I:tumor", paired.tumor_bam]
    params += ["--tumor_sample_name", paired.tumor_name]
    if paired.normal_bam is not None:
        params += ["-I:normal", paired.normal_bam]
        params += ["--normal_sample_name", paired.normal_name]
    if paired.normal_panel is not None:
        params += ["--normal_panel", paired.normal_panel]
    params += _config_params(base_config, assoc_files, region, out_file, items)
    return broad_runner, params
202958
def align_to_sort_bam(fastq1, fastq2, aligner, data):
    
    names = data["rgnames"]
    align_dir_parts = [data["dirs"]["work"], "align", names["sample"]]
    if data.get("disambiguate"):
        align_dir_parts.append(data["disambiguate"]["genome_build"])
    aligner_index = _get_aligner_index(aligner, data)
    align_dir = utils.safe_makedir(os.path.join(*align_dir_parts))
    ref_file = tz.get_in(("reference", "fasta", "base"), data)
    if fastq1.endswith(".bam"):
        data = _align_from_bam(fastq1, aligner, aligner_index, ref_file,
                               names, align_dir, data)
    else:
        data = _align_from_fastq(fastq1, fastq2, aligner, aligner_index, ref_file,
                                 names, align_dir, data)
    if data["work_bam"] and utils.file_exists(data["work_bam"]):
        if data.get("align_split") and dd.get_mark_duplicates(data):
            # If merging later with with bamsormadup need query sorted inputs
            # but CWL requires a bai file. Create a fake one to make it happy.
            bam.fake_index(data["work_bam"], data)
        else:
            bam.index(data["work_bam"], data["config"])
        for extra in ["-sr", "-disc"]:
            extra_bam = utils.append_stem(data['work_bam'], extra)
            if utils.file_exists(extra_bam):
                bam.index(extra_bam, data["config"])
    return data
1369847
def prepare(self):
        
        for s in self.get_all_sources():
            if s.missing():
                msg = "Source `%s' not found, needed by target `%s'."
                raise SCons.Errors.StopError(msg % (s, self.batches[0].targets[0]))
801867
def __prepare(self, filestem, b, examples=None, pos=None, neg=None):
        
        if examples:
            examplesFile = open('%s/%s.pl' % (self.tmpdir, filestem), 'w')
            examplesFile.write(examples)
            examplesFile.close()
        elif pos and neg:
            posFile = open('%s/%s.f' % (self.tmpdir, filestem), 'w')
            negFile = open('%s/%s.n' % (self.tmpdir, filestem), 'w')
            posFile.write(pos)
            negFile.write(neg)
            posFile.close()
            negFile.close()
        else:
            raise Exception('You need to provide either a single file of classified examples or \
                two files, positive and negative examples.')
        bFile = open('%s/%s.b' % (self.tmpdir, filestem), 'w')
        # Write settings.
        for setting, val in self.settings.items():
            bFile.write(':- set(%s,%s).\n' % (setting, val))
        bFile.write(b)
        bFile.close()
1240585
def prep_bam_inputs(out_dir, sample, call_file, bam_file):
    
    base = utils.splitext_plus(os.path.basename(bam_file))[0]
    with open(call_file) as in_handle:
        for cur_hla in (x.strip() for x in in_handle):
            out_file = os.path.join(utils.safe_makedir(os.path.join(out_dir, base)),
                                    "%s.type.%s.filtered.bam" % (base, cur_hla))
            if not os.path.exists(out_file):
                cmd = ["samtools", "view", "-b","-o", out_file, bam_file, cur_hla]
                subprocess.check_call(cmd)
499756
def picard_bam_to_fastq(picard, in_bam, fastq_one, fastq_two=None):
    
    if not file_exists(fastq_one):
        with tx_tmpdir(picard._config) as tmp_dir:
            with file_transaction(picard._config, fastq_one) as tx_out1:
                opts = [("INPUT", in_bam),
                        ("FASTQ", tx_out1),
                        ("TMP_DIR", tmp_dir)]
                if fastq_two is not None:
                    opts += [("SECOND_END_FASTQ", fastq_two)]
                picard.run("SamToFastq", opts)
    return (fastq_one, fastq_two)
1136328
def picard_reorder(picard, in_bam, ref_file, out_file):
    
    if not file_exists(out_file):
        with tx_tmpdir(picard._config) as tmp_dir:
            with file_transaction(picard._config, out_file) as tx_out_file:
                opts = [("INPUT", in_bam),
                        ("OUTPUT", tx_out_file),
                        ("REFERENCE", ref_file),
                        ("ALLOW_INCOMPLETE_DICT_CONCORDANCE", "true"),
                        ("TMP_DIR", tmp_dir)]
                picard.run("ReorderSam", opts)
    return out_file
326810
def process_ubam(bam, **kwargs):
    
    logging.info("Nanoget: Starting to collect statistics from ubam file {}.".format(bam))
    samfile = pysam.AlignmentFile(bam, "rb", check_sq=False)
    if not samfile.has_index():
        pysam.index(bam)
        # Need to reload the samfile after creating index
        samfile = pysam.AlignmentFile(bam, "rb")
        logging.info("Nanoget: No index for bam file could be found, created index.")
    datadf = pd.DataFrame(
        data=[(read.query_name, nanomath.ave_qual(read.query_qualities), read.query_length)
              for read in samfile.fetch(until_eof=True)],
        columns=["readIDs", "quals", "lengths"]) \
        .dropna(axis='columns', how='all') \
        .dropna(axis='index', how='any')
    logging.info("Nanoget: ubam {} contains {} reads.".format(
        bam, datadf["lengths"].size))
    return ut.reduce_memory_usage(datadf)
76734
def _setup_output(self, path, force):
        ''''''
        # Allow for special case of output being a pipe
        if os.path.isdir(path) or os.path.isfile(path):
            if force:
                logging.warn("Deleting previous file/directory '%s'" % path)
                if os.path.isfile(path):
                    os.remove(path)
                else:
                    shutil.rmtree(path)
            else:
                raise Exception("Cowardly refusing to overwrite already existing path at %s" % path)
405846
def _read_bam(bam_fn, precursors):
    
    mode = "r" if bam_fn.endswith("sam") else "rb"
    handle = pysam.Samfile(bam_fn, mode)
    reads = defaultdict(realign)
    for line in handle:
        chrom = handle.getrname(line.reference_id)
        # print("%s %s %s %s" % (line.query_name, line.reference_start, line.query_sequence, chrom))
        query_name = line.query_name
        if query_name not in reads:
            reads[query_name].sequence = line.query_sequence
        iso = isomir()
        iso.align = line
        iso.start = line.reference_start
        iso.subs, iso.add = _realign(reads[query_name].sequence, precursors[chrom], line.reference_start)
        reads[query_name].set_precursor(chrom, iso)

    reads = _clean_hits(reads)
    return reads
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
568743
def aligner_from_header(in_bam):
    
    from bcbio.pipeline.alignment import TOOLS
    with pysam.Samfile(in_bam, "rb") as bamfile:
        for pg in bamfile.header.get("PG", []):
            for ka in TOOLS.keys():
                if pg.get("PN", "").lower().find(ka) >= 0:
                    return ka