Convert SAM to BAM, including unmapped reads from fastq BAM file.
3524
def picard_sam_to_bam(picard, align_sam, fastq_bam, ref_file,
                      is_paired=False):
    
    to_retain = ["XS", "XG", "XM", "XN", "XO", "YT"]
    if align_sam.endswith(".sam"):
        out_bam = "%s.bam" % os.path.splitext(align_sam)[0]
    elif align_sam.endswith("-align.bam"):
        out_bam = "%s.bam" % align_sam.replace("-align.bam", "")
    else:
        raise NotImplementedError("Input format not recognized")
    if not file_exists(out_bam):
        with tx_tmpdir(picard._config) as tmp_dir:
            with file_transaction(picard._config, out_bam) as tx_out_bam:
                opts = [("UNMAPPED", fastq_bam),
                        ("ALIGNED", align_sam),
                        ("OUTPUT", tx_out_bam),
                        ("REFERENCE_SEQUENCE", ref_file),
                        ("TMP_DIR", tmp_dir),
                        ("PAIRED_RUN", ("true" if is_paired else "false")),
                        ]
                opts += [("ATTRIBUTES_TO_RETAIN", x) for x in to_retain]
                picard.run("MergeBamAlignment", opts)
    return out_bam
544667
def sampe(args, opts):
    
    dbfile, read1file, read2file = args
    dbfile = check_index(dbfile)
    sai1file = check_aln(dbfile, read1file, cpus=opts.cpus)
    sai2file = check_aln(dbfile, read2file, cpus=opts.cpus)

    samfile, _, unmapped = get_samfile(read1file, dbfile,
                                       bam=opts.bam, unmapped=opts.unmapped)
    if not need_update((dbfile, sai1file, sai2file), samfile):
        logging.error("`{0}` exists. `bwa samse` already run.".format(samfile))
        return "", samfile

    cmd = "bwa sampe " + " ".join((dbfile, sai1file, sai2file,
                                   read1file, read2file))
    cmd += " " + opts.extra
    if opts.cutoff:
        cmd += " -a {0}".format(opts.cutoff)
    if opts.uniq:
        cmd += " -n 1"

    return cmd, samfile
1277275
def _slice_bam(in_bam, region, tmp_dir, config):
    """"""
    name_file = os.path.splitext(os.path.basename(in_bam))[0]
    out_file = os.path.join(tmp_dir, os.path.join(tmp_dir, name_file + _to_str(region) + ".bam"))
    sambamba = config_utils.get_program("sambamba", config)
    region = _to_sambamba(region)
    with file_transaction(out_file) as tx_out_file:
        cmd = ("{sambamba} slice {in_bam} {region} -o {tx_out_file}")
        do.run(cmd.format(**locals()), "Slice region", {})
    return out_file
1185087
def get_reads(sam, \
        contigs = False, mismatches = False, mm_option = False, \
        sort_sam = True, req_map = False, region = False, sbuffer = False):
    
    tempdir = '%s/' % (os.path.abspath(sam).rsplit('/', 1)[0])
    if sort_sam is True:
        mapping = '%s.sorted.sam' % (sam.rsplit('.', 1)[0])
        if sam != '-':
            if os.path.exists(mapping) is False:
                os.system("\
                    sort -k1 --buffer-size=%sG -T %s -o %s %s\
                    " % (sbuffer, tempdir, mapping, sam))
        else:
            mapping = 'stdin-sam.sorted.sam'
            p = Popen("sort -k1 --buffer-size=%sG -T %s -o %s" \
                    % (sbuffer, tempdir, mapping), stdin = sys.stdin, shell = True)
            p.communicate()
        mapping = open(mapping)
    else:
        if sam == '-':
            mapping = sys.stdin
        else:
            mapping = open(sam)
    for read in reads_from_mapping(mapping, contigs, mismatches, mm_option, req_map, region):
        yield read
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
553544
def bam_to_fastq_bedtools(self, bam_file, out_fastq_pre, paired_end):
        
        self.make_sure_path_exists(os.path.dirname(out_fastq_pre))
        fq1 = out_fastq_pre + "_R1.fastq"
        fq2 = None
        cmd = self.tools.bedtools + " bamtofastq -i " + bam_file + " -fq " + fq1 + ".fastq"
        if paired_end:
            fq2 = out_fastq_pre + "_R2.fastq"
            cmd += " -fq2 " + fq2

        return cmd, fq1, fq2
1071523
def _get_all_unmapped_reads(self, fout):
        ''''''
        sam_reader = pysam.Samfile(self.bam, "rb")
        for read in sam_reader.fetch(until_eof=True):
            if read.is_unmapped:
                print(mapping.aligned_read_to_read(read, ignore_quality=not self.fastq_out), file=fout)
238083
def _exclude_region(self, contig, start, end, fout):
        ''''''
        sam_reader = pysam.Samfile(self.bam, "rb")
        exclude_interval = pyfastaq.intervals.Interval(start, end - 1)
        for read in sam_reader.fetch(contig):
            read_interval = pyfastaq.intervals.Interval(read.pos, read.reference_end - 1)
            if not read_interval.intersects(exclude_interval):
                print(mapping.aligned_read_to_read(read, ignore_quality=not self.fastq_out), file=fout)
917668
def _update_pysam_files(cf, destdir):
    ''''''
    basename = os.path.basename(destdir)
    for filename in cf:
        if not filename:
            continue
        dest = filename + ".pysam.c"
        with open(filename, encoding="utf-8") as infile:
            lines = "".join(infile.readlines())

            with open(dest, "w", encoding="utf-8") as outfile:
                outfile.write('#include "{}.pysam.h"\n\n'.format(basename))
                subname, _ = os.path.splitext(os.path.basename(filename))
                if subname in MAIN.get(basename, []):
                    lines = re.sub(r"int main\(", "int {}_main(".format(
                        basename), lines)
                else:
                    lines = re.sub(r"int main\(", "int {}_{}_main(".format(
                        basename, subname), lines)
                lines = re.sub("stderr", "{}_stderr".format(basename), lines)
                lines = re.sub("stdout", "{}_stdout".format(basename), lines)
                lines = re.sub(r" printf\(", " fprintf({}_stdout, ".format(basename), lines)
                lines = re.sub(r"([^kf])puts\(", r"\1{}_puts(".format(basename), lines)
                lines = re.sub(r"putchar\(([^)]+)\)",
                               r"fputc(\1, {}_stdout)".format(basename), lines)

                fn = os.path.basename(filename)
                # some specific fixes:
                SPECIFIC_SUBSTITUTIONS = {
                    "bam_md.c": (
                        'sam_open_format("-", mode_w',
                        'sam_open_format({}_stdout_fn, mode_w'.format(basename)),
                    "phase.c": (
                        'putc("ACGT"[f->seq[j] == 1? (c&3, {}_stdout) : (c>>16&3)]);'.format(basename),
                        'putc("ACGT"[f->seq[j] == 1? (c&3) : (c>>16&3)], {}_stdout);'.format(basename)),
                    "cut_target.c": (
                        'putc(33 + (cns[j]>>8>>2, {}_stdout));'.format(basename),
                        'putc(33 + (cns[j]>>8>>2), {}_stdout);'.format(basename))
                    }
                if fn in SPECIFIC_SUBSTITUTIONS:
                    lines = lines.replace(
                        SPECIFIC_SUBSTITUTIONS[fn][0],
                        SPECIFIC_SUBSTITUTIONS[fn][1])
                outfile.write(lines)

    with open(os.path.join("import", "pysam.h")) as inf, \
         open(os.path.join(destdir, "{}.pysam.h".format(basename)), "w") as outf:
        outf.write(re.sub("@pysam@", basename, inf.read()))

    with open(os.path.join("import", "pysam.c")) as inf, \
         open(os.path.join(destdir, "{}.pysam.c".format(basename)), "w") as outf:
        outf.write(re.sub("@pysam@", basename, inf.read()))
1340134
def _keep_assembled_chrom(bam_file, genome, config):
    """"""
    fai = "%s.fai" % genome
    chrom = []
    with open(fai) as inh:
        for line in inh:
            c = line.split("\t")[0]
            if c.find("_") < 0:
                chrom.append(c)
    chroms = " ".join(chrom)
    out_file = utils.append_stem(bam_file, '_chrom')
    samtools = config_utils.get_program("samtools", config)
    if not utils.file_exists(out_file):
        with file_transaction(out_file) as tx_out:
            cmd = "{samtools} view -b {bam_file} {chroms} > {tx_out}"
            do.run(cmd.format(**locals()), "Remove contigs from %s" % bam_file)
        bam.index(out_file, config)
    return out_file
254349
def sam_conversions(self, sam_file, depth=True):
        
        cmd = self.tools.samtools + " view -bS " + sam_file + " > " + sam_file.replace(".sam", ".bam") + "\n"
        cmd += self.tools.samtools + " sort " + sam_file.replace(".sam", ".bam") + " -o " + sam_file.replace(".sam", "_sorted.bam") + "\n"
        cmd += self.tools.samtools + " index " + sam_file.replace(".sam", "_sorted.bam") + "\n"
        if depth:
            cmd += self.tools.samtools + " depth " + sam_file.replace(".sam", "_sorted.bam") + " > " + sam_file.replace(".sam", "_sorted.depth") + "\n"
        return cmd
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
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
991145
def bamtag(sam):
    ''' 
    '''
    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
    if is_python3():
        queryalignment = next(track)
    else:
        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)

    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:
            aln.tags += [('XC', match.group('CB'))]
        if "molecular" in annotations:
            aln.tags += [('RX', match.group('MB'))]
        if "sample" in annotations:
            aln.tags += [('XS', match.group('SB'))]

        out_file.write(aln)

    total_time = time.time() - start_time
    logger.info('BAM tag conversion done - {:.3}s, {:,} alns/min'.format(total_time, int(60. * count / total_time)))
    logger.info("Processed %d alignments." % count)
251426
def _cram_to_fastq_region(cram_file, work_dir, base_name, region, data):
    
    ref_file = tz.get_in(["reference", "fasta", "base"], data)
    resources = config_utils.get_resources("bamtofastq", data["config"])
    cores = tz.get_in(["config", "algorithm", "num_cores"], data, 1)
    max_mem = config_utils.convert_to_bytes(resources.get("memory", "1G")) * cores
    rext = "-%s" % region.replace(":", "_").replace("-", "_") if region else "full"
    out_s, out_p1, out_p2, out_o1, out_o2 = [os.path.join(work_dir, "%s%s-%s.fq.gz" %
                                                          (base_name, rext, fext))
                                             for fext in ["s1", "p1", "p2", "o1", "o2"]]
    if not utils.file_exists(out_p1):
        with file_transaction(data, out_s, out_p1, out_p2, out_o1, out_o2) as \
             (tx_out_s, tx_out_p1, tx_out_p2, tx_out_o1, tx_out_o2):
            cram_file = objectstore.cl_input(cram_file)
            sortprefix = "%s-sort" % utils.splitext_plus(tx_out_s)[0]
            cmd = ("bamtofastq filename={cram_file} inputformat=cram T={sortprefix} "
                   "gz=1 collate=1 colsbs={max_mem} exclude=SECONDARY,SUPPLEMENTARY "
                   "F={tx_out_p1} F2={tx_out_p2} S={tx_out_s} O={tx_out_o1} O2={tx_out_o2} "
                   "reference={ref_file}")
            if region:
                cmd += " ranges='{region}'"
            do.run(cmd.format(**locals()), "CRAM to fastq %s" % region if region else "")
    return [[out_p1, out_p2, out_s]]
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
600522
def bam(args):
    
    from jcvi.formats.sizes import Sizes
    from jcvi.formats.sam import index

    p = OptionParser(bam.__doc__)
    p.set_home("eddyyeh")
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    gsnapfile, fastafile = args
    EYHOME = opts.eddyyeh_home
    pf = gsnapfile.rsplit(".", 1)[0]
    uniqsam = pf + ".unique.sam"
    samstats = uniqsam + ".stats"
    sizesfile = Sizes(fastafile).filename
    if need_update((gsnapfile, sizesfile), samstats):
        cmd = op.join(EYHOME, "gsnap2gff3.pl")
        cmd += " --format sam -i {0} -o {1}".format(gsnapfile, uniqsam)
        cmd += " -u -l {0} -p {1}".format(sizesfile, opts.cpus)
        sh(cmd)

    index([uniqsam])

    return uniqsam