Use GATK to extract reads from full BAM file.
7422
def _gatk_extract_reads_cl(data, region, prep_params, tmp_dir):
    
    args = ["PrintReads",
            "-L", region_to_gatk(region),
            "-R", dd.get_ref_file(data),
            "-I", data["work_bam"]]
    # GATK3 back compatibility, need to specify analysis type
    if "gatk4" in dd.get_tools_off(data):
        args = ["--analysis_type"] + args
    runner = broad.runner_from_config(data["config"])
    return runner.cl_gatk(args, tmp_dir)
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
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))
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
641892
def _get_region(self, contig, start, end, fout, min_length=250):
        ''''''
        sam_reader = pysam.Samfile(self.bam, "rb")
        trimming_end = (start == 0)
        for read in sam_reader.fetch(contig, start, end):
            read_interval = pyfastaq.intervals.Interval(read.pos, read.reference_end - 1)
            seq = mapping.aligned_read_to_read(read, ignore_quality=not self.fastq_out, revcomp=False)

            if trimming_end:
                bases_off_start = 0
                bases_off_end = max(0, read.reference_end - 1 - end)
                #seq.seq = seq.seq[:read.query_alignment_end - bases_off_end]
                seq = seq.subseq(0, read.query_alignment_end - bases_off_end)
            else:
                bases_off_start = max(0, start - read.pos + 1)
                #seq.seq = seq.seq[bases_off_start  + read.query_alignment_start:]
                seq = seq.subseq(bases_off_start  + read.query_alignment_start, len(seq))

            if read.is_reverse:
                seq.revcomp()

            if len(seq) >= min_length:
                print(seq, file=fout)
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
518751
def filter_reads(self, input_bam, output_bam, metrics_file, paired=False, cpus=16, Q=30):
        
        nodups = re.sub("\.bam$", "", output_bam) + ".nodups.nofilter.bam"
        cmd1 = self.tools.sambamba + " markdup -t {0} -r --compression-level=0 {1} {2} 2> {3}".format(cpus, input_bam, nodups, metrics_file)
        cmd2 = self.tools.sambamba + ' view -t {0} -f bam --valid'.format(cpus)
        if paired:
            cmd2 += ' -F "not (unmapped or mate_is_unmapped) and proper_pair'
        else:
            cmd2 += ' -F "not unmapped'
        cmd2 += ' and not (secondary_alignment or supplementary) and mapping_quality >= {0}"'.format(Q)
        cmd2 += ' {0} |'.format(nodups)
        cmd2 += self.tools.sambamba + " sort -t {0} /dev/stdin -o {1}".format(cpus, output_bam)
        cmd3 = "if [[ -s {0} ]]; then rm {0}; fi".format(nodups)
        cmd4 = "if [[ -s {0} ]]; then rm {0}; fi".format(nodups + ".bai")
        return [cmd1, cmd2, cmd3, cmd4]
848059
def main(self):
        
        # Find the target files
        self.targets()
        kmer = 15 if self.analysistype == 'GDCS' else 17
        # Use bbduk to bait the FASTQ reads matching the target sequences
        self.bait(maskmiddle='t', k=kmer)
        # If desired, use bbduk to bait the target sequences with the previously baited FASTQ files
        if self.revbait:
            self.reversebait(maskmiddle='t', k=kmer)
        # Run the bowtie2 read mapping module
        self.mapping()
        # Use samtools to index the sorted bam file
        self.indexing()
        # Parse the results
        # self.parsing()
        self.parsebam()
673621
def _goleft_indexcov(bam_file, data, out_dir):
    
    if not dd.get_coverage_interval(data) == "genome":
        return []
    out_dir = utils.safe_makedir(os.path.join(out_dir, "indexcov"))
    out_files = [os.path.join(out_dir, "%s-indexcov.%s" % (dd.get_sample_name(data), ext))
                 for ext in ["roc", "ped", "bed.gz"]]
    if not utils.file_uptodate(out_files[-1], bam_file):
        with transaction.tx_tmpdir(data) as tmp_dir:
            tmp_dir = utils.safe_makedir(os.path.join(tmp_dir, dd.get_sample_name(data)))
            gender_chroms = [x.name for x in ref.file_contigs(dd.get_ref_file(data)) if chromhacks.is_sex(x.name)]
            gender_args = "--sex %s" % (",".join(gender_chroms)) if gender_chroms else ""
            cmd = "goleft indexcov --directory {tmp_dir} {gender_args} -- {bam_file}"
            try:
                do.run(cmd.format(**locals()), "QC: goleft indexcov")
            except subprocess.CalledProcessError as msg:
                if not ("indexcov: no usable" in str(msg) or
                        ("indexcov: expected" in str(msg) and "sex chromosomes, found:" in str(msg))):
                    raise
            for out_file in out_files:
                orig_file = os.path.join(tmp_dir, os.path.basename(out_file))
                if utils.file_exists(orig_file):
                    utils.copy_plus(orig_file, out_file)
    # MultiQC needs non-gzipped/BED inputs so unpack the file
    out_bed = out_files[-1].replace(".bed.gz", ".tsv")
    if utils.file_exists(out_files[-1]) and not utils.file_exists(out_bed):
        with transaction.file_transaction(data, out_bed) as tx_out_bed:
            cmd = "gunzip -c %s > %s" % (out_files[-1], tx_out_bed)
            do.run(cmd, "Unpack indexcov BED file")
    out_files[-1] = out_bed
    return [x for x in out_files if utils.file_exists(x)]
580998
def bam_to_fastq_pair(in_file, target_region, pair):
    
    space, start, end = target_region
    bam_file = pysam.Samfile(in_file, "rb")
    for read in bam_file:
        if (not read.is_unmapped and not read.mate_is_unmapped
                and bam_file.getrname(read.tid) == space
                and bam_file.getrname(read.mrnm) == space
                and read.pos >= start and read.pos <= end
                and read.mpos >= start and read.mpos <= end
                and not read.is_secondary
                and read.is_paired and getattr(read, "is_read%s" % pair)):
            seq = Seq.Seq(read.seq)
            qual = list(read.qual)
            if read.is_reverse:
                seq = seq.reverse_complement()
                qual.reverse()
            yield read.qname, str(seq), "".join(qual)
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)
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)
1018157
def mutect2_caller(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):
        paired = vcfutils.get_paired_bams(align_bams, items)
        broad_runner = broad.runner_from_config(items[0]["config"])
        gatk_type = broad_runner.gatk_type()
        _prep_inputs(align_bams, ref_file, items)
        with file_transaction(items[0], out_file) as tx_out_file:
            params = ["-T", "Mutect2" if gatk_type == "gatk4" else "MuTect2",
                      "--annotation", "ClippingRankSumTest",
                      "--annotation", "DepthPerSampleHC"]
            if gatk_type == "gatk4":
                params += ["--reference", ref_file]
            else:
                params += ["-R", ref_file]
            for a in annotation.get_gatk_annotations(items[0]["config"], include_baseqranksum=False):
                params += ["--annotation", a]
            # Avoid issues with BAM CIGAR reads that GATK doesn't like
            if gatk_type == "gatk4":
                params += ["--read-validation-stringency", "LENIENT"]
            params += _add_tumor_params(paired, items, gatk_type)
            params += _add_region_params(region, out_file, items, gatk_type)
            # Avoid adding dbSNP/Cosmic so they do not get fed to variant filtering algorithm
            # Not yet clear how this helps or hurts in a general case.
            #params += _add_assoc_params(assoc_files)
            resources = config_utils.get_resources("mutect2", items[0]["config"])
            if "options" in resources:
                params += [str(x) for x in resources.get("options", [])]
            assert LooseVersion(broad_runner.gatk_major_version()) >= LooseVersion("3.5"), \
                "Require full version of GATK 3.5+ for mutect2 calling"
            broad_runner.new_resources("mutect2")
            gatk_cmd = broad_runner.cl_gatk(params, os.path.dirname(tx_out_file))
            if gatk_type == "gatk4":
                tx_raw_prefilt_file = "%s-raw%s" % utils.splitext_plus(tx_out_file)
                tx_raw_file = "%s-raw-filt%s" % utils.splitext_plus(tx_out_file)
                filter_cmd = _mutect2_filter(broad_runner, tx_raw_prefilt_file, tx_raw_file, ref_file)
                cmd = "{gatk_cmd} -O {tx_raw_prefilt_file} && {filter_cmd}"
            else:
                tx_raw_file = "%s-raw%s" % utils.splitext_plus(tx_out_file)
                cmd = "{gatk_cmd} > {tx_raw_file}"
            do.run(cmd.format(**locals()), "MuTect2")
            out_file = _af_filter(paired.tumor_data, tx_raw_file, out_file)
    return vcfutils.bgzip_and_index(out_file, items[0]["config"])
160448
def to_bam(in_file, out_file, data):
    
    if not utils.file_uptodate(out_file, in_file):
        with file_transaction(data, out_file) as tx_out_file:
            cmd = ["samtools", "view", "-O", "BAM", "-o", tx_out_file, in_file]
            do.run(cmd, "Convert CRAM to BAM")
    bam.index(out_file, data["config"])
    return out_file
947891
def read_alignment(out_sam, loci, seqs, out_file):
    """"""
    hits = defaultdict(list)
    with open(out_file, "w") as out_handle:
        samfile = pysam.Samfile(out_sam, "r")
        for a in samfile.fetch():
            if not a.is_unmapped:
                nm = int([t[1] for t in a.tags if t[0] == "NM"][0])
                a = makeBED(a)
                if not a:
                    continue
                ref, locus = get_loci(samfile.getrname(int(a.chr)), loci)
                hits[a.name].append((nm, "%s %s %s %s %s %s" % (a.name, a.name.split("-")[0], locus, ref, a.start, a.end)))
        for hit in hits.values():
            nm = hit[0][0]
            for l in hit:
                if nm == l[0]:
                    print(l[1], file=out_handle)
    return out_file
179517
def _clean_alignment(args):
    
    logger.info("Clean bam file with highly repetitive reads with low counts. sum(counts)/n_hits > 1%")
    bam_file, seq_obj = clean_bam_file(args.afile, args.mask)
    logger.info("Using %s file" % bam_file)
    detect_complexity(bam_file, args.ref, args.out)
    return bam_file, seq_obj
1054963
def extract_dms (archive, compression, cmd, verbosity, interactive, outdir):
    """"""
    check_archive_ext(archive)
    cmdlist = [cmd, '-d', outdir]
    if verbosity > 1:
        cmdlist.append('-v')
    cmdlist.extend(['u', archive])
    return cmdlist