Prepare a sample to be run, potentially converting from BAM to
    FASTQ and/or downsampling the number of reads for a test run
761331
def prepare_sample(data):
    
    data = utils.to_single_data(data)
    logger.debug("Preparing %s" % data["rgnames"]["sample"])
    data["files"] = get_fastq_files(data)
    # get_fastq_files swaps over quality scores to standard, unless trimming
    if not(dd.get_trim_reads(data)):
        data = dd.set_quality_format(data, "standard")
    return [[data]]
1240585
def prep_bam_inputs(out_dir, sample, call_file, bam_file):
    
    base = utils.splitext_plus(os.path.basename(bam_file))[0]
    with open(call_file) as in_handle:
        for cur_hla in (x.strip() for x in in_handle):
            out_file = os.path.join(utils.safe_makedir(os.path.join(out_dir, base)),
                                    "%s.type.%s.filtered.bam" % (base, cur_hla))
            if not os.path.exists(out_file):
                cmd = ["samtools", "view", "-b","-o", out_file, bam_file, cur_hla]
                subprocess.check_call(cmd)
1062248
def bam_to_fastq(bam_file, is_paired):
    
    out_files, out_handles = _get_fastq_handles(bam_file,
            is_paired)
    if len(out_handles) > 0:
        in_bam = pysam.Samfile(bam_file, mode='rb')
        for read in in_bam:
            num = 1 if (not read.is_paired or read.is_read1) else 2
            # reverse the sequence and quality if mapped to opposite strand
            if read.is_reverse:
                seq = str(Seq.reverse_complement(Seq.Seq(read.seq)))
                qual = "".join(reversed(read.qual))
            else:
                seq = read.seq
                qual = read.qual
            out_handles[num].write("@%s\n%s\n+\n%s\n" % (read.qname,
                seq, qual))
    [h.close() for h in out_handles.values()]
    return out_files
451861
def _run_germline(bam_file, data, ref_file, region, out_file):
    
    work_dir = utils.safe_makedir("%s-work" % utils.splitext_plus(out_file)[0])
    region_bed = strelka2.get_region_bed(region, [data], out_file, want_gzip=False)
    example_dir = _make_examples(bam_file, data, ref_file, region_bed, out_file, work_dir)
    if _has_candidate_variants(example_dir):
        tfrecord_file = _call_variants(example_dir, region_bed, data, out_file)
        return _postprocess_variants(tfrecord_file, data, ref_file, out_file)
    else:
        return vcfutils.write_empty_vcf(out_file, data["config"], [dd.get_sample_name(data)])
1141350
def main():
    '''
    
    '''
    args = get_args()
    try:
        utils.make_output_dir(args.outdir)
        utils.init_logs(args)
        args.format = nanoplotter.check_valid_format(args.format)
        settings = vars(args)
        settings["path"] = path.join(args.outdir, args.prefix)
        sources = {
            "fastq": args.fastq,
            "bam": args.bam,
            "cram": args.cram,
            "fastq_rich": args.fastq_rich,
            "fastq_minimal": args.fastq_minimal,
            "summary": args.summary,
            "fasta": args.fasta,
            "ubam": args.ubam,
        }

        if args.pickle:
            datadf = pickle.load(open(args.pickle, 'rb'))
        else:
            datadf = get_input(
                source=[n for n, s in sources.items() if s][0],
                files=[f for f in sources.values() if f][0],
                threads=args.threads,
                readtype=args.readtype,
                combine="simple",
                barcoded=args.barcoded)
        if args.store:
            pickle.dump(
                obj=datadf,
                file=open(settings["path"] + "NanoPlot-data.pickle", 'wb'))
        if args.raw:
            datadf.to_csv(settings["path"] + "NanoPlot-data.tsv.gz",
                          sep="\t",
                          index=False,
                          compression="gzip")

        settings["statsfile"] = [make_stats(datadf, settings, suffix="")]
        datadf, settings = filter_and_transform_data(datadf, settings)
        if settings["filtered"]:  # Bool set when filter was applied in filter_and_transform_data()
            settings["statsfile"].append(make_stats(datadf, settings, suffix="_post_filtering"))

        if args.barcoded:
            barcodes = list(datadf["barcode"].unique())
            plots = []
            for barc in barcodes:
                logging.info("Processing {}".format(barc))
                settings["path"] = path.join(args.outdir, args.prefix + barc + "_")
                dfbarc = datadf[datadf["barcode"] == barc]
                if len(dfbarc) > 5:
                    settings["title"] = barc
                    plots.extend(
                        make_plots(dfbarc, settings)
                    )
                else:
                    sys.stderr.write("Found barcode {} less than 5x, ignoring...\n".format(barc))
                    logging.info("Found barcode {} less than 5 times, ignoring".format(barc))
        else:
            plots = make_plots(datadf, settings)
        make_report(plots, settings)
        logging.info("Finished!")
    except Exception as e:
        logging.error(e, exc_info=True)
        print("\n\n\nIf you read this then NanoPlot has crashed :-(")
        print("Please try updating NanoPlot and see if that helps...\n")
        print("If not, please report this issue at https://github.com/wdecoster/NanoPlot/issues")
        print("If you could include the log file that would be really helpful.")
        print("Thanks!\n\n\n")
        raise
1286730
def downsample(f1, f2, data, N, quick=False):
    
    if quick:
        rand_records = range(N)
    else:
        records = sum(1 for _ in open(f1)) / 4
        N = records if N > records else N
        rand_records = sorted(random.sample(xrange(records), N))

    fh1 = open_possible_gzip(f1)
    fh2 = open_possible_gzip(f2) if f2 else None
    outf1 = os.path.splitext(f1)[0] + ".subset" + os.path.splitext(f1)[1]
    outf2 = os.path.splitext(f2)[0] + ".subset" + os.path.splitext(f2)[1] if f2 else None

    if utils.file_exists(outf1):
        if not outf2:
            return outf1, outf2
        elif utils.file_exists(outf2):
            return outf1, outf2

    out_files = (outf1, outf2) if outf2 else (outf1)

    with file_transaction(out_files) as tx_out_files:
        if isinstance(tx_out_files, six.string_types):
            tx_out_f1 = tx_out_files
        else:
            tx_out_f1, tx_out_f2 = tx_out_files
        sub1 = open_possible_gzip(tx_out_f1, "w")
        sub2 = open_possible_gzip(tx_out_f2, "w") if outf2 else None
        rec_no = - 1
        for rr in rand_records:
            while rec_no < rr:
                rec_no += 1
                for i in range(4): fh1.readline()
                if fh2:
                    for i in range(4): fh2.readline()
            for i in range(4):
                sub1.write(fh1.readline())
                if sub2:
                    sub2.write(fh2.readline())
            rec_no += 1
        fh1.close()
        sub1.close()
        if f2:
            fh2.close()
            sub2.close()

    return outf1, outf2
634899
def _merge_fastqc(samples):
    
    fastqc_list = collections.defaultdict(list)
    seen = set()
    for data in samples:
        name = dd.get_sample_name(data)
        if name in seen:
            continue
        seen.add(name)
        fns = glob.glob(os.path.join(dd.get_work_dir(data), "qc", dd.get_sample_name(data), "fastqc") + "/*")
        for fn in fns:
            if fn.endswith("tsv"):
                metric = os.path.basename(fn)
                fastqc_list[metric].append([name, fn])

    for metric in fastqc_list:
        dt_by_sample = []
        for fn in fastqc_list[metric]:
            dt = pd.read_csv(fn[1], sep="\t")
            dt['sample'] = fn[0]
            dt_by_sample.append(dt)
        dt = utils.rbind(dt_by_sample)
        dt.to_csv(metric, sep="\t", index=False, mode ='w')
    return samples
474788
def _new_sample(self, sink):
        
        if self.running:
            # New data is available from the pipeline! Dump it into our
            # queue (or possibly block if we're full).
            buf = sink.emit('pull-sample').get_buffer()

            # We can't use Gst.Buffer.extract() to read the data as it crashes
            # when called through PyGObject. We also can't use
            # Gst.Buffer.extract_dup() because we have no way in Python to free
            # the memory that it returns. Instead we get access to the actual
            # data via Gst.Memory.map().
            mem = buf.get_all_memory()
            success, info = mem.map(Gst.MapFlags.READ)
            if success:
                data = info.data
                mem.unmap(info)
                self.queue.put(data)
            else:
                raise GStreamerError("Unable to map buffer memory while reading the file.")
        return Gst.FlowReturn.OK
1357394
def _create_sample_list(in_bams, vcf_file):
    
    out_file = "%s-sample_list.txt" % os.path.splitext(vcf_file)[0]
    with open(out_file, "w") as out_handle:
        for in_bam in in_bams:
            with pysam.Samfile(in_bam, "rb") as work_bam:
                for rg in work_bam.header.get("RG", []):
                    out_handle.write("%s\n" % rg["SM"])
    return out_file
1050825
def SampleStart(self):
    """"""
    self._start_cpu_time = time.clock()
    self.start_sample_time = time.time()
    self.total_cpu_time = 0
1312137
def compress(samples, run_parallel):
    
    to_cram = []
    finished = []
    for data in [x[0] for x in samples]:
        if "cram" in dd.get_archive(data) or "cram-lossless" in dd.get_archive(data):
            to_cram.append([data])
        else:
            finished.append([data])
    crammed = run_parallel("archive_to_cram", to_cram)
    return finished + crammed
1211310
def _create_sample_file(data, out_dir):
    """"""
    sample_file = os.path.join(out_dir, "sample_file.txt")
    with open(sample_file, 'w') as outh:
        for sample in data:
            outh.write(sample[0]["clean_fastq"] + "\n")
    return sample_file
1009590
def _sample_template(sample, out_dir):
    """"""
    bam_fn = dd.get_work_bam(sample)
    genome = dd.get_genome_build(sample)
    if genome in supported:
        peaks = sample.get("peaks_files", []).get("main")
        if peaks:
            r_code = ("library(ChIPQC);\n"
                      "sample = ChIPQCsample(\"{bam_fn}\","
                      "\"{peaks}\", "
                      "annotation = \"{genome}\","
                      ");\n"
                      "ChIPQCreport(sample);\n")
            r_code_fn = os.path.join(out_dir, "chipqc.r")
            with open(r_code_fn, 'w') as inh:
                inh.write(r_code.format(**locals()))
            return r_code_fn
801867
def __prepare(self, filestem, b, examples=None, pos=None, neg=None):
        
        if examples:
            examplesFile = open('%s/%s.pl' % (self.tmpdir, filestem), 'w')
            examplesFile.write(examples)
            examplesFile.close()
        elif pos and neg:
            posFile = open('%s/%s.f' % (self.tmpdir, filestem), 'w')
            negFile = open('%s/%s.n' % (self.tmpdir, filestem), 'w')
            posFile.write(pos)
            negFile.write(neg)
            posFile.close()
            negFile.close()
        else:
            raise Exception('You need to provide either a single file of classified examples or \
                two files, positive and negative examples.')
        bFile = open('%s/%s.b' % (self.tmpdir, filestem), 'w')
        # Write settings.
        for setting, val in self.settings.items():
            bFile.write(':- set(%s,%s).\n' % (setting, val))
        bFile.write(b)
        bFile.close()
147115
function() {
      let sample = clock(start)

      if (test.after) test.after()
      bench.after.forEach(call)

      if (warmups == -1) {
        warmups = Math.max(minWarmups, sample > 100 ? 1 : sample > 10 ? 5 : 50)
      }
      if (warmups > 0) {
        warmups--
      } else {
        samples[n++] = sample
        onSample(samples[n - 1], cycle)

        if (minTime <= clock(t) && minSamples <= n) {
          return cycle.done() // all done!
        }
      }

      // wait then repeat
      wait(delay).then(next)
    }
92032
def _make_examples(bam_file, data, ref_file, region_bed, out_file, work_dir):
    
    log_dir = utils.safe_makedir(os.path.join(work_dir, "log"))
    example_dir = utils.safe_makedir(os.path.join(work_dir, "examples"))
    if len(glob.glob(os.path.join(example_dir, "%s.tfrecord*.gz" % dd.get_sample_name(data)))) == 0:
        with tx_tmpdir(data) as tx_example_dir:
            cmd = ["dv_make_examples.py", "--cores", dd.get_num_cores(data), "--ref", ref_file,
                   "--reads", bam_file, "--regions", region_bed, "--logdir", log_dir,
                   "--examples", tx_example_dir, "--sample", dd.get_sample_name(data)]
            do.run(cmd, "DeepVariant make_examples %s" % dd.get_sample_name(data))
            for fname in glob.glob(os.path.join(tx_example_dir, "%s.tfrecord*.gz" % dd.get_sample_name(data))):
                utils.copy_plus(fname, os.path.join(example_dir, os.path.basename(fname)))
    return example_dir
812614
def run_mosdepth(data, target_name, bed_file, per_base=False, quantize=None, thresholds=None):
    
    MosdepthCov = collections.namedtuple("MosdepthCov", ("dist", "per_base", "regions", "quantize", "thresholds"))
    bam_file = dd.get_align_bam(data) or dd.get_work_bam(data)
    work_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "coverage", dd.get_sample_name(data)))
    prefix = os.path.join(work_dir, "%s-%s" % (dd.get_sample_name(data), target_name))
    old_dist_file = "%s.mosdepth.dist.txt" % (prefix)
    out = MosdepthCov((old_dist_file if utils.file_uptodate(old_dist_file, bam_file) else
                       "%s.mosdepth.%s.dist.txt" % (prefix, "region" if bed_file else "global")),
                      ("%s.per-base.bed.gz" % prefix) if per_base else None,
                      ("%s.regions.bed.gz" % prefix) if bed_file else None,
                      ("%s.quantized.bed.gz" % prefix) if quantize else None,
                      ("%s.thresholds.bed.gz" % prefix) if thresholds else None)
    if not utils.file_uptodate(out.dist, bam_file):
        with file_transaction(data, out.dist) as tx_out_file:
            tx_prefix = os.path.join(os.path.dirname(tx_out_file), os.path.basename(prefix))
            num_cores = dd.get_cores(data)
            bed_arg = ("--by %s" % bed_file) if bed_file else ""
            perbase_arg = "" if per_base else "--no-per-base"
            mapq_arg = "-Q 1" if (per_base or quantize) else ""
            if quantize:
                quant_arg = "--quantize %s" % quantize[0]
                quant_export = " && ".join(["export MOSDEPTH_Q%s=%s" % (i, x) for (i, x) in enumerate(quantize[1])])
                quant_export += " && "
            else:
                quant_arg, quant_export = "", ""

            thresholds_cmdl = ("-T " + ",".join([str(t) for t in thresholds])) if out.thresholds else ""
            cmd = ("{quant_export}mosdepth -t {num_cores} -F 1804 {mapq_arg} {perbase_arg} {bed_arg} {quant_arg} "
                   "{tx_prefix} {bam_file} {thresholds_cmdl}")
            message = "Calculating coverage: %s %s" % (dd.get_sample_name(data), target_name)
            do.run(cmd.format(**locals()), message.format(**locals()))
            if out.per_base:
                shutil.move(os.path.join(os.path.dirname(tx_out_file), os.path.basename(out.per_base)), out.per_base)
            if out.regions:
                shutil.move(os.path.join(os.path.dirname(tx_out_file), os.path.basename(out.regions)), out.regions)
            if out.quantize:
                shutil.move(os.path.join(os.path.dirname(tx_out_file), os.path.basename(out.quantize)), out.quantize)
            if out.thresholds:
                shutil.move(os.path.join(os.path.dirname(tx_out_file), os.path.basename(out.thresholds)), out.thresholds)
    return out