Convert CRAM file into BAM.
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
1109902
def ref_file_from_bam(bam_file, data):
    
    new_ref = os.path.join(utils.safe_makedir(os.path.join(dd.get_work_dir(data), "inputs", "ref")),
                           "%s-subset.fa" % dd.get_genome_build(data))
    if not utils.file_exists(new_ref):
        with file_transaction(data, new_ref) as tx_out_file:
            contig_file = "%s-contigs.txt" % utils.splitext_plus(new_ref)[0]
            with open(contig_file, "w") as out_handle:
                for contig in [x.contig for x in idxstats(bam_file, data) if x.contig != "*"]:
                    out_handle.write("%s\n" % contig)
            cmd = "seqtk subseq -l 100 %s %s > %s" % (dd.get_ref_file(data), contig_file, tx_out_file)
            do.run(cmd, "Subset %s to BAM file contigs" % dd.get_genome_build(data))
    ref.fasta_idx(new_ref, data["config"])
    runner = broad.runner_from_path("picard", data["config"])
    runner.run_fn("picard_index_ref", new_ref)
    return {"base": new_ref}
277520
function parseBamHeader(r) {
        if (!r) {
            return callback(null, "Couldn't access BAM");
        }

        var unc = unbgzf(r, r.byteLength);
        var uncba = new Uint8Array(unc);

        var magic = readInt(uncba, 0);
        if (magic != BAM_MAGIC) {
            return callback(null, "Not a BAM file, magic=0x" + magic.toString(16));
        }
        var headLen = readInt(uncba, 4);
        var header = '';
        for (var i = 0; i < headLen; ++i) {
            header += String.fromCharCode(uncba[i + 8]);
        }

        var nRef = readInt(uncba, headLen + 8);
        var p = headLen + 12;

        bam.chrToIndex = {};
        bam.indexToChr = [];
        for (var i = 0; i < nRef; ++i) {
            var lName = readInt(uncba, p);
            var name = '';
            for (var j = 0; j < lName-1; ++j) {
                name += String.fromCharCode(uncba[p + 4 + j]);
            }
            var lRef = readInt(uncba, p + lName + 4);
            bam.chrToIndex[name] = i;
            if (name.indexOf('chr') == 0) {
                bam.chrToIndex[name.substring(3)] = i;
            } else {
                bam.chrToIndex['chr' + name] = i;
            }
            bam.indexToChr.push(name);

            p = p + 8 + lName;
        }

        if (bam.indices) {
            return callback(bam);
        }
    }
29086
def get_chrs_from_bam(self, file_name):
        
        x = subprocess.check_output(self.tools.samtools + " view -H " + file_name + " | grep '^@SQ' | cut -f2| sed s'/SN://'", shell=True)
        # Chromosomes will be separated by newlines; split into list to return
        return x.split()
105474
def validate_bam(self, input_bam):
        
        cmd = self.tools.java + " -Xmx" + self.pm.javamem
        cmd += " -jar " + self.tools.picard + " ValidateSamFile"
        cmd += " INPUT=" + input_bam
        return cmd
79060
public static void checkBamOrCramFile(InputStream is, String bamFileName, boolean checkSort) throws IOException {
        SamReaderFactory srf = SamReaderFactory.make();
        srf.validationStringency(ValidationStringency.LENIENT);

        SamReader reader = srf.open(SamInputResource.of(is));
        SAMFileHeader fileHeader = reader.getFileHeader();
        SAMFileHeader.SortOrder sortOrder = fileHeader.getSortOrder();
        reader.close();

        if (reader.type().equals(SamReader.Type.SAM_TYPE)) {
            throw new IOException("Expected binary SAM file. File " + bamFileName + " is not binary.");
        }

        if (checkSort) {
            switch (sortOrder) {
                case coordinate:
                    break;
                case queryname:
                case unsorted:
                default:
                    throw new IOException("Expected sorted file. File '" + bamFileName + "' is not sorted by coordinates("
                            + sortOrder.name() + ")");
            }
        }
    }
1092272
def parsebam(self):
        
        # Threading is actually the worst - need multiprocessing to make this work at all
        logging.info('Parsing BAM files')
        # The sample objects are too big to get pickled. To hack our way around this, try to dump the sample object to
        # json, and have the processing function turn the object into a dictionary.
        json_files = list()
        with tempfile.TemporaryDirectory() as tmpdir:
            best_assemblies = list()
            sample_names = list()
            for sample in self.runmetadata:
                json_name = os.path.join(tmpdir, '{sn}.json'.format(sn=sample.name))
                best_assemblies.append(sample.general.bestassemblyfile)
                sample_names.append(sample.name)
                with open(json_name, 'w') as f:
                    json.dump(sample[self.analysistype].dump(), f, sort_keys=True, indent=4)
                json_files.append(json_name)
            p = multiprocessing.Pool(processes=self.cpus)
            analysis_type_list = [self.analysistype] * len(self.runmetadata)
            iupac_list = [self.iupac] * len(self.runmetadata)
            cutoff_list = [self.cutoff] * len(self.runmetadata)
            depth_list = [self.averagedepth] * len(self.runmetadata)
            allow_soft_clip_list = [self.allow_soft_clips] * len(self.runmetadata)
            sample_results = p.starmap(Sippr.parse_one_sample,
                                       zip(json_files, sample_names, best_assemblies, analysis_type_list,
                                           iupac_list, cutoff_list, depth_list, allow_soft_clip_list))
            p.close()
            p.join()
        # Since we had to json-ize the sample objects, we now need to update the metadata for everything.
        for sample in self.runmetadata:
                sample[self.analysistype].faidict = dict()
                sample[self.analysistype].results = dict()
                sample[self.analysistype].avgdepth = dict()
                sample[self.analysistype].resultssnp = dict()
                sample[self.analysistype].snplocations = dict()
                sample[self.analysistype].resultsgap = dict()
                sample[self.analysistype].gaplocations = dict()
                sample[self.analysistype].sequences = dict()
                sample[self.analysistype].maxcoverage = dict()
                sample[self.analysistype].mincoverage = dict()
                sample[self.analysistype].standarddev = dict()
                # Figure out which of the sample results to use.
                for sample_result in sample_results:
                    if sample_result['name'] == sample.name:
                        sample[self.analysistype].faidict = sample_result['faidict']
                        sample[self.analysistype].results = sample_result['results']
                        sample[self.analysistype].avgdepth = sample_result['avgdepth']
                        sample[self.analysistype].resultssnp = sample_result['resultssnp']
                        sample[self.analysistype].snplocations = sample_result['snplocations']
                        sample[self.analysistype].resultsgap = sample_result['resultsgap']
                        sample[self.analysistype].gaplocations = sample_result['gaplocations']
                        sample[self.analysistype].sequences = sample_result['sequences']
                        sample[self.analysistype].maxcoverage = sample_result['maxcoverage']
                        sample[self.analysistype].mincoverage = sample_result['mincoverage']
                        sample[self.analysistype].standarddev = sample_result['standarddev']
        logging.info('Done parsing BAM files')
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
648352
def rm(self, filename):
        """"""
        command = """
            try:
                import os
            except ImportError:
                import uos as os
            os.remove('{0}')
        """.format(
            filename
        )
        self._pyboard.enter_raw_repl()
        try:
            out = self._pyboard.exec_(textwrap.dedent(command))
        except PyboardError as ex:
            message = ex.args[2].decode("utf-8")
            # Check if this is an OSError #2, i.e. file/directory doesn't exist
            # and rethrow it as something more descriptive.
            if message.find("OSError: [Errno 2] ENOENT") != -1:
                raise RuntimeError("No such file/directory: {0}".format(filename))
            # Check for OSError #13, the directory isn't empty.
            if message.find("OSError: [Errno 13] EACCES") != -1:
                raise RuntimeError("Directory is not empty: {0}".format(filename))
            else:
                raise ex
        self._pyboard.exit_raw_repl()
653719
def _extract(archive, compression, cmd, format, verbosity, outdir):
    """"""
    targetname = util.get_single_outfile(outdir, archive)
    try:
        with lzma.LZMAFile(archive, **_get_lzma_options(format)) as lzmafile:
            with open(targetname, 'wb') as targetfile:
                data = lzmafile.read(READ_SIZE_BYTES)
                while data:
                    targetfile.write(data)
                    data = lzmafile.read(READ_SIZE_BYTES)
    except Exception as err:
        msg = "error extracting %s to %s: %s" % (archive, targetname, err)
        raise util.PatoolError(msg)
    return None
356819
def to_sdf(files, data):
    
    # BAM
    if len(files) == 1 and files[0].endswith(".bam"):
        qual = []
        format = ["-f", "sam-pe" if bam.is_paired(files[0]) else "sam-se"]
        inputs = [files[0]]
    # fastq
    else:
        qual = ["-q", "illumina" if dd.get_quality_format(data).lower() == "illumina" else "sanger"]
        format = ["-f", "fastq"]
        if len(files) == 2:
            inputs = ["-l", files[0], "-r", files[1]]
        else:
            assert len(files) == 1
            inputs = [files[0]]
    work_dir = utils.safe_makedir(os.path.join(data["dirs"]["work"], "align_prep"))
    out_file = os.path.join(work_dir,
                            "%s.sdf" % utils.splitext_plus(os.path.basename(os.path.commonprefix(files)))[0])
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            cmd = _rtg_cmd(["rtg", "format", "-o", tx_out_file] + format + qual + inputs)
            do.run(cmd, "Format inputs to indexed SDF")
    return out_file
558370
def md5sum(filename):
    """"""
    hash_md5 = hashlib.md5()
    with open(filename, 'rb') as handle:
        for chunk in iter(lambda: handle.read(4096), b''):
            hash_md5.update(chunk)
    return hash_md5.hexdigest()
173480
def convert(self, value, param, ctx):  # pylint: disable=inconsistent-return-statements
        """"""
        matches = MEMORY_RE.match(value.lower())
        if matches is None:
            self.fail('%s is not a valid value for memory amount' % value, param, ctx)
        amount_str, unit = matches.groups()
        amount = int(amount_str)
        if unit in [None, 'm', 'mb']:
            # Assume the user intends gigabytes if they specify a number < 1024
            if amount < 1024:
                return amount * 1024
            else:
                if amount % 1024 != 0:
                    self.fail('%s is not an integer that is divisable by 1024' % value, param, ctx)
                return amount
        elif unit in ['g', 'gb']:
            return amount * 1024
802724
def convert_bam_to_fastq(in_file, work_dir, data, dirs, config):
    
    return alignprep.prep_fastq_inputs([in_file], data)
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
976711
def open_pysam_file(fname, ftype):
    """"""
    try:
        if ftype == 'bam':
            fpysam = pysam.AlignmentFile(fname, 'rb')
        elif ftype == 'fasta':
            fpysam = pysam.FastaFile(fname)
        yield fpysam
    except:
        raise
    else:
        fpysam.close()
1020339
def bed(args):
    
    p = OptionParser(bed.__doc__)
    opts, args = p.parse_args(args)

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

    bedfile = args[0]
    bamfiles = args[1:]
    for bamfile in bamfiles:
        cmd = "bamToBed -i {0}".format(bamfile)
        sh(cmd, outfile=bedfile, append=True)