return only one BAM file with only chromosome 21
973400
def _slice_bam_chr21(in_bam, data):
    
    sambamba = config_utils.get_program("sambamba", data["config"])
    out_file = "%s-chr%s" % os.path.splitext(in_bam)
    if not utils.file_exists(out_file):
        bam.index(in_bam, data['config'])
        with pysam.Samfile(in_bam, "rb") as bamfile:
            bam_contigs = [c["SN"] for c in bamfile.header["SQ"]]
        chromosome = "21"
        if "chr21" in bam_contigs:
            chromosome = "chr21"
        with file_transaction(data, out_file) as tx_out_file:
            cmd = ("{sambamba} slice -o {tx_out_file} {in_bam} {chromosome}").format(**locals())
            out = subprocess.check_output(cmd, shell=True)
    return out_file
686150
def check_bam(bam, o):
    
    try:
        p = sp.Popen(['samtools', 'view', bam], stdout=sp.PIPE)
        # Count paired alignments
        paired = 0
        read_lengths = defaultdict(int)
        while o > 0:  # Count down number of lines
            line = p.stdout.readline().decode().split("\t")
            flag = int(line[1])
            read_lengths[len(line[9])] += 1
            if 1 & flag:  # check decimal flag contains 1 (paired)
                paired += 1
            o -= 1
        p.kill()
    except OSError:
        reason = "Note (samtools not in path): For NGS inputs, " \
                 "pep needs samtools to auto-populate " \
                 "'read_length' and 'read_type' attributes; " \
                 "these attributes were not populated."
        raise OSError(reason)

    _LOGGER.debug("Read lengths: {}".format(read_lengths))
    _LOGGER.debug("paired: {}".format(paired))
    return read_lengths, paired
523320
def tobam_cl(data, out_file, is_paired=False):
    
    do_dedup = _check_dedup(data)
    umi_consensus = dd.get_umi_consensus(data)
    with file_transaction(data, out_file) as tx_out_file:
        if not do_dedup:
            yield (sam_to_sortbam_cl(data, tx_out_file), tx_out_file)
        elif umi_consensus:
            yield (_sam_to_grouped_umi_cl(data, umi_consensus, tx_out_file), tx_out_file)
        elif is_paired and _need_sr_disc_reads(data) and not _too_many_contigs(dd.get_ref_file(data)):
            sr_file = "%s-sr.bam" % os.path.splitext(out_file)[0]
            disc_file = "%s-disc.bam" % os.path.splitext(out_file)[0]
            with file_transaction(data, sr_file) as tx_sr_file:
                with file_transaction(data, disc_file) as tx_disc_file:
                    yield (samblaster_dedup_sort(data, tx_out_file, tx_sr_file, tx_disc_file),
                           tx_out_file)
        else:
            yield (_biobambam_dedup_sort(data, tx_out_file), tx_out_file)
954107
def _check_bam_contigs(in_bam, ref_file, config):
    
    # GATK allows chromosome M to be in multiple locations, skip checking it
    allowed_outoforder = ["chrM", "MT"]
    ref_contigs = [c.name for c in ref.file_contigs(ref_file, config)]
    with pysam.Samfile(in_bam, "rb") as bamfile:
        bam_contigs = [c["SN"] for c in bamfile.header["SQ"]]
    extra_bcs = [x for x in bam_contigs if x not in ref_contigs]
    extra_rcs = [x for x in ref_contigs if x not in bam_contigs]
    problems = []
    warnings = []
    for bc, rc in zip_longest([x for x in bam_contigs if (x not in extra_bcs and
                                                                     x not in allowed_outoforder)],
                                         [x for x in ref_contigs if (x not in extra_rcs and
                                                                     x not in allowed_outoforder)]):
        if bc != rc:
            if bc and rc:
                problems.append("Reference mismatch. BAM: %s Reference: %s" % (bc, rc))
            elif bc:
                warnings.append("Extra BAM chromosomes: %s" % bc)
            elif rc:
                warnings.append("Extra reference chromosomes: %s" % rc)
    for bc in extra_bcs:
        warnings.append("Extra BAM chromosomes: %s" % bc)
    for rc in extra_rcs:
        warnings.append("Extra reference chromosomes: %s" % rc)
    if problems:
        raise ValueError("Unexpected order, name or contig mismatches between input BAM and reference file:\n%s\n"
                         "Setting `bam_clean: remove_extracontigs` in the configuration can often fix this issue."
                         % "\n".join(problems))
    if warnings:
        print("*** Potential problems in input BAM compared to reference:\n%s\n" %
              "\n".join(warnings))
922210
def _check_umi_type(bam_file):
    
    with pysam.Samfile(bam_file, "rb") as in_bam:
        for read in in_bam:
            cur_umi = None
            for tag in ["RX", "XC"]:
                try:
                    cur_umi = read.get_tag(tag)
                    break
                except KeyError:
                    pass
            if cur_umi:
                if "-" in cur_umi and len(cur_umi.split("-")) == 2:
                    return "paired", tag
                else:
                    return "adjacency", tag
585055
def run(name, chip_bam, input_bam, genome_build, out_dir, method, resources, data):
    
    # output file name need to have the caller name
    config = dd.get_config(data)
    out_file = os.path.join(out_dir, name + "_peaks_macs2.xls")
    macs2_file = os.path.join(out_dir, name + "_peaks.xls")
    if utils.file_exists(out_file):
        _compres_bdg_files(out_dir)
        return _get_output_files(out_dir)
    macs2 = config_utils.get_program("macs2", config)
    options = " ".join(resources.get("macs2", {}).get("options", ""))
    genome_size = bam.fasta.total_sequence_length(dd.get_ref_file(data))
    genome_size = "" if options.find("-g") > -1 else "-g %s" % genome_size
    paired = "-f BAMPE" if bam.is_paired(chip_bam) else ""
    with utils.chdir(out_dir):
        cmd = _macs2_cmd(method)
        try:
            do.run(cmd.format(**locals()), "macs2 for %s" % name)
            utils.move_safe(macs2_file, out_file)
        except subprocess.CalledProcessError:
            raise RuntimeWarning("macs2 terminated with an error.\n"
                                 "Please, check the message and report "
                                 "error if it is related to bcbio.\n"
                                 "You can add specific options for the sample "
                                 "setting resources as explained in docs: "
                                 "https://bcbio-nextgen.readthedocs.org/en/latest/contents/configuration.html#sample-specific-resources")
    _compres_bdg_files(out_dir)
    return _get_output_files(out_dir)
334006
def get_bams():
    
    for size, md5, url in bams:
        bam = os.path.join(
            args.data_dir,
            os.path.basename(url).replace('.bam', '_%s.bam' % CHROM))
        if not _up_to_date(md5, bam):
            logger.info(
                'Downloading reads on chromosome %s from %s to %s'
                % (CHROM, url, bam))
            cmds = ['samtools', 'view', '-b', url, COORD, '>', bam]
            logged_command(cmds)
        bai = bam + '.bai'
        if not os.path.exists(bai):
            logger.info('indexing %s' % bam)
            logger.info(' '.join(cmds))
            cmds = [
                'samtools',
                'index',
                bam]
            logged_command(cmds)
        if os.path.exists(os.path.basename(url) + '.bai'):
            os.unlink(os.path.basename(url) + '.bai')

    for size, md5, fn in bais:
        if not _up_to_date(md5, fn):
            cmds = [
                'samtools', 'index', bai.replace('.bai', '')]
            logged_command(cmds)
3608
def _target_chroms_and_header(bam_file, data):
    
    special_remaps = {"chrM": "MT", "MT": "chrM"}
    target_chroms = dict([(x.name, i) for i, x in enumerate(ref.file_contigs(dd.get_ref_file(data)))
                          if chromhacks.is_autosomal_or_sex(x.name)])
    out_chroms = []
    with pysam.Samfile(bam_file, "rb") as bamfile:
        for bami, bam_contig in enumerate([c["SN"] for c in bamfile.header["SQ"]]):
            if bam_contig in target_chroms:
                target_chrom = bam_contig
            elif bam_contig in special_remaps and special_remaps[bam_contig] in target_chroms:
                target_chrom = special_remaps[bam_contig]
            elif bam_contig.startswith("chr") and bam_contig.replace("chr", "") in target_chroms:
                target_chrom = bam_contig.replace("chr", "")
            elif "chr%s" % bam_contig in target_chroms:
                target_chrom = "chr%s" % bam_contig
            else:
                target_chrom = None
            # target_chrom == bam_contig ensures we don't try chr1 -> 1 style remapping
            if target_chrom and target_chrom == bam_contig:
                # Order not required if dealing with SAM file header fixing
                #assert bami == target_chroms[target_chrom], \
                #    ("remove_extracontigs: Non-matching order of standard contig: %s %s (%s vs %s)" %
                #     (bam_file, target_chrom, bami, target_chroms[target_chrom]))
                out_chroms.append(target_chrom)
    assert out_chroms, ("remove_extracontigs: Did not find any chromosomes in reference file: %s %s" %
                        (bam_file, target_chroms))
    return out_chroms
305931
def batchlobstr(args):
    
    p = OptionParser(batchlobstr.__doc__)
    p.add_option("--haploid", default="chrY,chrM",
                 help="Use haploid model for these chromosomes")
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    bamlist, = args
    cmd = "python -m jcvi.variation.str lobstr TREDs"
    cmd += " --input_bam_path {}"
    cmd += " --haploid {}".format(opts.haploid)
    cmd += " --simulation"
    cmds = [cmd.format(x.strip()) for x in open(bamlist).readlines()]
    p = Parallel(cmds, cpus=opts.cpus)
    p.run()
991154
def bam2es(
        bam_fn,
        es_fo,
        allowed_delta,
    ):
        

        es_fo.write("# RN:   read name" + os.linesep)
        es_fo.write("# Q:    is mapped with quality" + os.linesep)
        es_fo.write("# Chr:  chr id" + os.linesep)
        es_fo.write("# D:    direction" + os.linesep)
        es_fo.write("# L:    leftmost nucleotide" + os.linesep)
        es_fo.write("# R:    rightmost nucleotide" + os.linesep)
        es_fo.write("# Cat:  category of alignment assigned by LAVEnder" + os.linesep)
        es_fo.write("#         M_i    i-th segment is correctly mapped" + os.linesep)
        es_fo.write("#         m      segment should be unmapped but it is mapped" + os.linesep)
        es_fo.write("#         w      segment is mapped to a wrong location" + os.linesep)
        es_fo.write("#         U      segment is unmapped and should be unmapped" + os.linesep)
        es_fo.write("#         u      segment is unmapped and should be mapped" + os.linesep)
        es_fo.write("# Segs: number of segments" + os.linesep)
        es_fo.write("# " + os.linesep)
        es_fo.write("# RN\tQ\tChr\tD\tL\tR\tCat\tSegs" + os.linesep)

        with pysam.AlignmentFile(bam_fn, "rb") as sam:
            references_dict = {}

            for i in range(len(sam.references)):
                references_dict[sam.references[i]] = i + 1

            for read in sam:
                rnf_read_tuple = rnftools.rnfformat.ReadTuple()
                rnf_read_tuple.destringize(read.query_name)

                left = read.reference_start + 1
                right = read.reference_end
                chrom_id = references_dict[sam.references[read.reference_id]]

                nb_of_segments = len(rnf_read_tuple.segments)

                if rnf_read_tuple.segments[0].genome_id == 1:
                    should_be_mapped = True
                else:
                    should_be_mapped = False

                # read is unmapped
                if read.is_unmapped:
                    # read should be mapped
                    if should_be_mapped:
                        category = "u"
                    # read should be unmapped
                    else:
                        category = "U"
                # read is mapped
                else:
                    # read should be mapped
                    if should_be_mapped:
                        exists_corresponding_segment = False

                        for j in range(len(rnf_read_tuple.segments)):
                            segment = rnf_read_tuple.segments[j]
                            if (
                                (segment.left == 0 or abs(segment.left - left) <= allowed_delta)
                                and (segment.right == 0 or abs(segment.right - right) <= allowed_delta)
                                and (segment.left != 0 or segment.right == 0)
                                and (chrom_id == 0 or chrom_id == segment.chr_id)
                            ):
                                exists_corresponding_segment = True
                                segment = str(j + 1)
                                break

                        # read was mapped to correct location
                        if exists_corresponding_segment:  # exists ok location?
                            category = "M_" + segment
                        # read was mapped to incorrect location
                        else:
                            category = "w"
                    # read should be unmapped
                    else:
                        category = "m"

                es_fo.write(
                    "\t".join(
                        map(
                            str,
                            [
                                # read name
                                read.query_name,
                                # aligned?
                                "unmapped" if read.is_unmapped else "mapped_" + str(read.mapping_quality),
                                # reference id
                                chrom_id,
                                # direction
                                "R" if read.is_reverse else "F",
                                # left
                                left,
                                # right
                                right,
                                # assigned category
                                category,
                                # count of segments
                                nb_of_segments
                            ]
                        )
                    ) + os.linesep
                )
1025750
def prepare_exclude_file(items, base_file, chrom=None):
    
    items = shared.add_highdepth_genome_exclusion(items)
    out_file = "%s-exclude%s.bed" % (utils.splitext_plus(base_file)[0], "-%s" % chrom if chrom else "")
    if not utils.file_exists(out_file) and not utils.file_exists(out_file + ".gz"):
        with shared.bedtools_tmpdir(items[0]):
            with file_transaction(items[0], out_file) as tx_out_file:
                # Get a bedtool for the full region if no variant regions
                want_bedtool = callable.get_ref_bedtool(tz.get_in(["reference", "fasta", "base"], items[0]),
                                                        items[0]["config"], chrom)
                want_bedtool = pybedtools.BedTool(shared.subset_variant_regions(want_bedtool.saveas().fn,
                                                                                chrom, tx_out_file, items))
                sv_exclude_bed = _get_sv_exclude_file(items)
                if sv_exclude_bed and len(want_bedtool) > 0:
                    want_bedtool = want_bedtool.subtract(sv_exclude_bed, nonamecheck=True).saveas()
                full_bedtool = callable.get_ref_bedtool(tz.get_in(["reference", "fasta", "base"], items[0]),
                                                        items[0]["config"])
                if len(want_bedtool) > 0:
                    full_bedtool.subtract(want_bedtool, nonamecheck=True).saveas(tx_out_file)
                else:
                    full_bedtool.saveas(tx_out_file)
    return out_file
511859
def run(bam_file, data, out_dir):
    """"""
    m = {"base": None, "secondary": []}
    m.update(_mirbase_stats(data, out_dir))
    m["secondary"].append(_seqcluster_stats(data, out_dir))
1192527
def get_hla_choice(h, hlas, normal_bam, tumor_bam):
    
    def get_counts(bam_file):
        counts = {}
        for line in subprocess.check_output(["samtools", "idxstats", bam_file]).split("\n"):
            if line.startswith(h):
                name, _, count, _ = line.split()
                counts[name] = int(count)
        return counts
    tcounts = get_counts(tumor_bam)
    ncounts = get_counts(normal_bam)
    check_hlas = [x for x in hlas if x.startswith(h) and tcounts.get(x, 0) > 0 and ncounts.get(x, 0) > 0]
    cur_hlas = sorted(check_hlas, key=lambda x: ncounts[x], reverse=True)
    #print(cur_hlas[0], tcounts.get(cur_hlas[0]), ncounts.get(cur_hlas[0]))
    return cur_hlas[0]
1326332
def subset_bam_by_region(in_file, region, config, out_file_base=None):
    
    if out_file_base is not None:
        base, ext = os.path.splitext(out_file_base)
    else:
        base, ext = os.path.splitext(in_file)
    out_file = "%s-subset%s%s" % (base, region, ext)
    if not file_exists(out_file):
        with pysam.Samfile(in_file, "rb") as in_bam:
            target_tid = in_bam.gettid(region)
            assert region is not None, \
                   "Did not find reference region %s in %s" % \
                   (region, in_file)
            with file_transaction(config, out_file) as tx_out_file:
                with pysam.Samfile(tx_out_file, "wb", template=in_bam) as out_bam:
                    for read in in_bam:
                        if read.tid == target_tid:
                            out_bam.write(read)
    return out_file
727207
def get_minibam_bed(bamfile, bedfile, minibam=None):
    
    pf = op.basename(bedfile).split(".")[0]
    minibamfile = minibam or op.basename(bamfile).replace(".bam", ".{}.bam".format(pf))
    minisamfile = minibam.replace(".bam", ".sam")
    baifile = minibamfile + ".bai"
    if op.exists(baifile):
        sh("rm {}".format(baifile))

    cmd = "samtools view -H {} > {}".format(bamfile, minisamfile)
    sh(cmd)

    cmd = "cat {}".format(bedfile)
    cmd += " | perl -lane 'print \"$F[0]:$F[1]-$F[2]\"'"
    cmd += " | xargs -n1 -t -I \{\}"
    cmd += " samtools view {}".format(bamfile)
    cmd += " \{\} >> " + minisamfile
    sh(cmd)

    cmd = "samtools view {} -b".format(minisamfile)
    cmd += " | samtools sort -"
    cmd += " -o {0}".format(minibamfile)

    sh(cmd)
    sh("samtools index {0}".format(minibamfile))
    return minibamfile
1132967
def filter_bam(job, job_vars):
    
    input_args, ids = job_vars
    work_dir = job.fileStore.getLocalTempDir()
    cores = input_args['cpu_count']
    sudo = input_args['sudo']
    # I/O
    transcriptome_bam = return_input_paths(job, work_dir, ids, 'transcriptome.bam')
    output = os.path.join(work_dir, 'filtered.bam')
    # Command
    parameters = ['sam-filter',
                  '--strip-indels',
                  '--max-insert', '1000',
                  '--mapq', '1',
                  '--in', docker_path(transcriptome_bam),
                  '--out', docker_path(output)]
    docker_call(tool='quay.io/ucsc_cgl/ubu:1.2--02806964cdf74bf5c39411b236b4c4e36d026843',
                tool_parameters=parameters, work_dir=os.path.dirname(output), java_opts='-Xmx30g', sudo=sudo)
    # Write to FileStore
    ids['filtered.bam'] = job.fileStore.writeGlobalFile(output)
    # Run child job
    return job.addChildJobFn(rsem, job_vars, cores=cores, disk='30 G').rv()
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