%prog mini bamfile region

    Extract mini-bam for a single region.
175768
def mini(args):
    
    p = OptionParser(mini.__doc__)
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    bamfile, region = args
    get_minibam(bamfile, region)
1300715
def _bgzip_from_bam(bam_file, dirs, data, is_retry=False, output_infix=''):
    
    # tools
    config = data["config"]
    bamtofastq = config_utils.get_program("bamtofastq", config)
    resources = config_utils.get_resources("bamtofastq", config)
    cores = config["algorithm"].get("num_cores", 1)
    max_mem = config_utils.convert_to_bytes(resources.get("memory", "1G")) * cores
    bgzip = tools.get_bgzip_cmd(config, is_retry)
    # files
    work_dir = utils.safe_makedir(os.path.join(dirs["work"], "align_prep"))
    out_file_1 = os.path.join(work_dir, "%s%s-1.fq.gz" % (os.path.splitext(os.path.basename(bam_file))[0], output_infix))
    out_file_2 = out_file_1.replace("-1.fq.gz", "-2.fq.gz")
    needs_retry = False
    if is_retry or not utils.file_exists(out_file_1):
        if not bam.is_paired(bam_file):
            out_file_2 = None
        with file_transaction(config, out_file_1) as tx_out_file:
            for f in [tx_out_file, out_file_1, out_file_2]:
                if f and os.path.exists(f):
                    os.remove(f)
            fq1_bgzip_cmd = "%s -c /dev/stdin > %s" % (bgzip, tx_out_file)
            prep_cmd = _seqtk_fastq_prep_cl(data, read_num=0)
            if prep_cmd:
                fq1_bgzip_cmd = prep_cmd + " | " + fq1_bgzip_cmd
            sortprefix = "%s-sort" % os.path.splitext(tx_out_file)[0]
            if bam.is_paired(bam_file):
                prep_cmd = _seqtk_fastq_prep_cl(data, read_num=1)
                fq2_bgzip_cmd = "%s -c /dev/stdin > %s" % (bgzip, out_file_2)
                if prep_cmd:
                    fq2_bgzip_cmd = prep_cmd + " | " + fq2_bgzip_cmd
                out_str = ("F=>({fq1_bgzip_cmd}) F2=>({fq2_bgzip_cmd}) S=/dev/null O=/dev/null "
                           "O2=/dev/null collate=1 colsbs={max_mem}")
            else:
                out_str = "S=>({fq1_bgzip_cmd})"
            bam_file = objectstore.cl_input(bam_file)
            extra_opts = " ".join([str(x) for x in resources.get("options", [])])
            cmd = "{bamtofastq} filename={bam_file} T={sortprefix} {extra_opts} " + out_str
            try:
                do.run(cmd.format(**locals()), "BAM to bgzipped fastq",
                       checks=[do.file_reasonable_size(tx_out_file, bam_file)],
                       log_error=False)
            except subprocess.CalledProcessError as msg:
                if not is_retry and "deflate failed" in str(msg):
                    logger.info("bamtofastq deflate IO failure preparing %s. Retrying with single core."
                                % (bam_file))
                    needs_retry = True
                else:
                    logger.exception()
                    raise
    if needs_retry:
        return _bgzip_from_bam(bam_file, dirs, data, is_retry=True)
    else:
        return [x for x in [out_file_1, out_file_2] if x is not None and utils.file_exists(x)]
1358520
def run(align_bams, items, ref_file, assoc_files, region, out_file):
    
    call_file = "%s-raw.vcf.gz" % utils.splitext_plus(out_file)[0]
    strelka_work_dir = "%s-work" % utils.splitext_plus(out_file)[0]
    paired = vcfutils.get_paired_bams(align_bams, items)
    if paired:
        assert paired.normal_bam, "Strelka2 requires a normal sample"
        call_file = _run_somatic(paired, ref_file, assoc_files, region, call_file, strelka_work_dir)
    else:
        call_file = _run_germline(align_bams, items, ref_file,
                                  assoc_files, region, call_file, strelka_work_dir)
    return _af_annotate_and_filter(paired, items, call_file, out_file)
401443
def _piped_bamprep_region(data, region, out_file, tmp_dir):
    
    if _need_prep(data):
        prep_params = _get_prep_params(data)
        _piped_bamprep_region_gatk(data, region, prep_params, out_file, tmp_dir)
    else:
        raise ValueError("No realignment specified")
289331
def _check_sample(in_bam, rgnames):
    
    with pysam.Samfile(in_bam, "rb") as bamfile:
        rg = bamfile.header.get("RG", [{}])
    msgs = []
    warnings = []
    if len(rg) > 1:
        warnings.append("Multiple read groups found in input BAM. Expect single RG per BAM.")
    if len(rg) == 0:
        msgs.append("No read groups found in input BAM. Expect single RG per BAM.")
    if len(rg) > 0 and any(x.get("SM") != rgnames["sample"] for x in rg):
        msgs.append("Read group sample name (SM) does not match configuration `description`: %s vs %s"
                    % (rg[0].get("SM"), rgnames["sample"]))
    if len(msgs) > 0:
        raise ValueError("Problems with pre-aligned input BAM file: %s\n" % (in_bam)
                         + "\n".join(msgs) +
                         "\nSetting `bam_clean: fixrg`\n"
                         "in the configuration can often fix this issue.")
    if warnings:
        print("*** Potential problems in input BAM compared to reference:\n%s\n" %
              "\n".join(warnings))
299119
def detect_bam_features(bamfile, n_entries=1000):
    '''  '''

    inbam = pysam.Samfile(bamfile)
    inbam = inbam.fetch(until_eof=True)

    tags = ["NH", "X0", "XT"]
    available_tags = {x: 1 for x in tags}

    for n, read in enumerate(inbam):
        if n > n_entries:
            break

        if read.is_unmapped:
            continue

        else:
            for tag in tags:
                if not read.has_tag(tag):
                    available_tags[tag] = 0

    return available_tags
578539
def run_tnscope(align_bams, items, ref_file, assoc_files,
                     region=None, out_file=None):
    
    if out_file is None:
        out_file = "%s-variants.vcf.gz" % utils.splitext_plus(align_bams[0])[0]
    if not utils.file_exists(out_file):
        variant_regions = bedutils.population_variant_regions(items, merged=True)
        interval = _get_interval(variant_regions, region, out_file, items)
        with file_transaction(items[0], out_file) as tx_out_file:
            paired = vcfutils.get_paired_bams(align_bams, items)
            assert paired and paired.normal_bam, "Require normal BAM for Sentieon TNscope"
            dbsnp = "--dbsnp %s" % (assoc_files.get("dbsnp")) if "dbsnp" in assoc_files else ""
            license = license_export(items[0])
            cores = dd.get_num_cores(items[0])
            cmd = ("{license}sentieon driver -t {cores} -r {ref_file} "
                   "-i {paired.tumor_bam} -i {paired.normal_bam} {interval} "
                   "--algo TNscope "
                   "--tumor_sample {paired.tumor_name} --normal_sample {paired.normal_name} "
                   "{dbsnp} {tx_out_file}")
            do.run(cmd.format(**locals()), "Sentieon TNscope")
    return out_file
673922
def _get_nblock_regions(in_file, min_n_size, ref_regions):
    
    out_lines = []
    called_contigs = set([])
    with utils.open_gzipsafe(in_file) as in_handle:
        for line in in_handle:
            contig, start, end, ctype = line.rstrip().split()
            called_contigs.add(contig)
            if (ctype in ["REF_N", "NO_COVERAGE", "EXCESSIVE_COVERAGE", "LOW_COVERAGE"] and
                  int(end) - int(start) > min_n_size):
                out_lines.append("%s\t%s\t%s\n" % (contig, start, end))
    for refr in ref_regions:
        if refr.chrom not in called_contigs:
            out_lines.append("%s\t%s\t%s\n" % (refr.chrom, 0, refr.stop))
    return pybedtools.BedTool("\n".join(out_lines), from_string=True)
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)
290681
def is_empty(bam_file):
    
    bam_file = objectstore.cl_input(bam_file)
    cmd = ("set -o pipefail; "
           "samtools view {bam_file} | head -1 | wc -l")
    p = subprocess.Popen(cmd.format(**locals()), shell=True,
                         executable=do.find_bash(),
                         stdout=subprocess.PIPE, stderr=subprocess.PIPE,
                         preexec_fn=lambda: signal.signal(signal.SIGPIPE, signal.SIG_DFL))
    stdout, stderr = p.communicate()
    stdout = stdout.decode()
    stderr = stderr.decode()
    if ((p.returncode == 0 or p.returncode == 141) and
         (stderr == "" or (stderr.startswith("gof3r") and stderr.endswith("broken pipe")))):
        return int(stdout) == 0
    else:
        raise ValueError("Failed to check empty status of BAM file: %s" % str(stderr))
651403
def _prepare_smoove_bams(full_bams, sr_bams, disc_bams, items, tx_work_dir):
    
    input_dir = utils.safe_makedir(tx_work_dir)
    out = []
    for full, sr, disc, data in zip(full_bams, sr_bams, disc_bams, items):
        if sr and disc:
            new_full = os.path.join(input_dir, "%s.bam" % dd.get_sample_name(data))
            new_sr = os.path.join(input_dir, "%s.split.bam" % dd.get_sample_name(data))
            new_disc = os.path.join(input_dir, "%s.disc.bam" % dd.get_sample_name(data))
            utils.symlink_plus(full, new_full)
            utils.symlink_plus(sr, new_sr)
            utils.symlink_plus(disc, new_disc)
            out.append(new_full)
        else:
            out.append(full)
    return out
539486
def extract_ape (archive, compression, cmd, verbosity, interactive, outdir):
    """"""
    outfile = util.get_single_outfile(outdir, archive, extension=".wav")
    return [cmd, archive, outfile, '-d']
68580
def _split_by_ready_regions(ext, file_key, dir_ext_fn):
    
    def _sort_by_size(region_w_bams):
        region, _ = region_w_bams
        _, start, end = region
        return end - start
    def _assign_bams_to_regions(data):
        """Ensure BAMs aligned with input regions, either global or individual.
        """
        for i, region in enumerate(data["region"]):
            work_bams = []
            for xs in data["region_bams"]:
                if len(xs) == 1:
                    work_bams.append(xs[0])
                else:
                    work_bams.append(xs[i])
            for work_bam in work_bams:
                assert os.path.exists(work_bam), work_bam
            yield region, work_bams
    def _do_work(data):
        if "region" in data:
            name = data["group"][0] if "group" in data else data["description"]
            out_dir = os.path.join(data["dirs"]["work"], dir_ext_fn(data))
            out_file = os.path.join(out_dir, "%s%s" % (name, ext))
            assert isinstance(data["region"], (list, tuple))
            out_parts = []
            for r, work_bams in sorted(_assign_bams_to_regions(data), key=_sort_by_size, reverse=True):
                out_region_dir = os.path.join(out_dir, r[0])
                out_region_file = os.path.join(out_region_dir,
                                               "%s-%s%s" % (name, pregion.to_safestr(r), ext))
                out_parts.append((r, work_bams, out_region_file))
            return out_file, out_parts
        else:
            return None, []
    return _do_work
1145637
def minimac(args):
    
    p = OptionParser(minimac.__doc__)
    p.set_home("shapeit")
    p.set_home("minimac")
    p.set_outfile()
    p.set_chr()
    p.set_ref()
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    txtfile, = args
    ref = opts.ref
    mm = MakeManager()
    pf = txtfile.split(".")[0]
    allrawvcf = []
    alloutvcf = []
    chrs = opts.chr.split(",")
    for x in chrs:
        px = CM[x]
        chrvcf = pf + ".{0}.vcf".format(px)
        if txtfile.endswith(".vcf"):
            cmd = "vcftools --vcf {0} --chr {1}".format(txtfile, x)
            cmd += " --out {0}.{1} --recode".format(pf, px)
            cmd += " && mv {0}.{1}.recode.vcf {2}".format(pf, px, chrvcf)
        else:  # 23andme
            cmd = "python -m jcvi.formats.vcf from23andme {0} {1}".format(txtfile, x)
            cmd += " --ref {0}".format(ref)
        mm.add(txtfile, chrvcf, cmd)

        chrvcf_hg38 = pf + ".{0}.23andme.hg38.vcf".format(px)
        minimac_liftover(mm, chrvcf, chrvcf_hg38, opts)
        allrawvcf.append(chrvcf_hg38)

        minimacvcf = "{0}.{1}.minimac.dose.vcf".format(pf, px)
        if x == "X":
            minimac_X(mm, x, chrvcf, opts)
        elif x in ["Y", "MT"]:
            cmd = "python -m jcvi.variation.impute passthrough"
            cmd += " {0} {1}".format(chrvcf, minimacvcf)
            mm.add(chrvcf, minimacvcf, cmd)
        else:
            minimac_autosome(mm, x, chrvcf, opts)

        # keep the best line for multi-allelic markers
        uniqvcf= "{0}.{1}.minimac.uniq.vcf".format(pf, px)
        cmd = "python -m jcvi.formats.vcf uniq {0} > {1}".\
                    format(minimacvcf, uniqvcf)
        mm.add(minimacvcf, uniqvcf, cmd)

        minimacvcf_hg38 = "{0}.{1}.minimac.hg38.vcf".format(pf, px)
        minimac_liftover(mm, uniqvcf, minimacvcf_hg38, opts)
        alloutvcf.append(minimacvcf_hg38)

    if len(allrawvcf) > 1:
        rawhg38vcfgz = pf + ".all.23andme.hg38.vcf.gz"
        cmd = "vcf-concat {0} | bgzip > {1}".format(" ".join(allrawvcf), rawhg38vcfgz)
        mm.add(allrawvcf, rawhg38vcfgz, cmd)

    if len(alloutvcf) > 1:
        outhg38vcfgz = pf + ".all.minimac.hg38.vcf.gz"
        cmd = "vcf-concat {0} | bgzip > {1}".format(" ".join(alloutvcf), outhg38vcfgz)
        mm.add(alloutvcf, outhg38vcfgz, cmd)

    mm.write()
156792
def sample_name(in_bam):
    
    with pysam.AlignmentFile(in_bam, "rb", check_sq=False) as in_pysam:
        try:
            if "RG" in in_pysam.header:
                return in_pysam.header["RG"][0]["SM"]
        except ValueError:
            return None
533524
def isolate_region(sequences, start, end, gap_char='-'):
    
    # Check arguments
    if end <= start:
        raise ValueError("start of slice must precede end ({0} !> {1})".format(
            end, start))

    for sequence in sequences:
        seq = sequence.seq
        start_gap = gap_char * start
        end_gap = gap_char * (len(seq) - end)
        seq = Seq(start_gap + str(seq[start:end]) + end_gap,
                alphabet=seq.alphabet)
        sequence.seq = seq
        yield sequence
1315553
def run(bam_file, data, out_dir):
    
    out = dict()

    out_dir = utils.safe_makedir(out_dir)
    if dd.get_coverage(data) and dd.get_coverage(data) not in ["None"]:
        merged_bed_file = bedutils.clean_file(dd.get_coverage_merged(data), data, prefix="cov-", simple=True)
        target_name = "coverage"
    elif dd.get_coverage_interval(data) != "genome":
        merged_bed_file = dd.get_variant_regions_merged(data) or dd.get_sample_callable(data)
        target_name = "variant_regions"
    else:
        merged_bed_file = None
        target_name = "genome"

    avg_depth = cov.get_average_coverage(target_name, merged_bed_file, data)
    if target_name == "coverage":
        out_files = cov.coverage_region_detailed_stats(target_name, merged_bed_file, data, out_dir)
    else:
        out_files = []

    out['Avg_coverage'] = avg_depth

    samtools_stats_dir = os.path.join(out_dir, os.path.pardir, 'samtools')
    from bcbio.qc import samtools
    samtools_stats = samtools.run(bam_file, data, samtools_stats_dir)["metrics"]

    out["Total_reads"] = total_reads = int(samtools_stats["Total_reads"])
    out["Mapped_reads"] = mapped = int(samtools_stats["Mapped_reads"])
    out["Mapped_paired_reads"] = int(samtools_stats["Mapped_paired_reads"])
    out['Duplicates'] = dups = int(samtools_stats["Duplicates"])

    if total_reads:
        out["Mapped_reads_pct"] = 100.0 * mapped / total_reads
    if mapped:
        out['Duplicates_pct'] = 100.0 * dups / mapped

    if dd.get_coverage_interval(data) == "genome":
        mapped_unique = mapped - dups
    else:
        mapped_unique = readstats.number_of_mapped_reads(data, bam_file, keep_dups=False)
    out['Mapped_unique_reads'] = mapped_unique

    if merged_bed_file:
        ontarget = readstats.number_of_mapped_reads(
            data, bam_file, keep_dups=False, bed_file=merged_bed_file, target_name=target_name)
        out["Ontarget_unique_reads"] = ontarget
        if mapped_unique:
            out["Ontarget_pct"] = 100.0 * ontarget / mapped_unique
            out['Offtarget_pct'] = 100.0 * (mapped_unique - ontarget) / mapped_unique
            if dd.get_coverage_interval(data) != "genome":
                # Skip padded calculation for WGS even if the "coverage" file is specified
                # the padded statistic makes only sense for exomes and panels
                padded_bed_file = bedutils.get_padded_bed_file(out_dir, merged_bed_file, 200, data)
                ontarget_padded = readstats.number_of_mapped_reads(
                    data, bam_file, keep_dups=False, bed_file=padded_bed_file, target_name=target_name + "_padded")
                out["Ontarget_padded_pct"] = 100.0 * ontarget_padded / mapped_unique
        if total_reads:
            out['Usable_pct'] = 100.0 * ontarget / total_reads

    indexcov_files = _goleft_indexcov(bam_file, data, out_dir)
    out_files += [x for x in indexcov_files if x and utils.file_exists(x)]
    out = {"metrics": out}
    if len(out_files) > 0:
        out["base"] = out_files[0]
        out["secondary"] = out_files[1:]
    return out