Create a final output file with merged CNVkit and THetA copy and population estimates.
183568
def _merge_theta_calls(bounds_file, result_file, cnv_file, data):
    
    out_file = "%s-merged.txt" % (result_file.replace(".BEST.results", ""))
    if not utils.file_uptodate(out_file, result_file):
        with file_transaction(data, out_file) as tx_out_file:
            updater = _update_with_calls(result_file, cnv_file)
            with open(bounds_file) as in_handle:
                with open(tx_out_file, "w") as out_handle:
                    i = 0
                    for line in in_handle:
                        if line.startswith("#"):
                            parts = line.rstrip().split("\t")
                            parts += ["cnv", "pop_cnvs", "pop_pvals"]
                            out_handle.write("\t".join(parts) + "\n")
                        else:
                            out_handle.write(updater(i, line))
                            i += 1
    return out_file
236239
def _finalize_sv(solution_file, data):
    
    out = {"variantcaller": "titancna"}
    with open(solution_file) as in_handle:
        solution = dict(zip(in_handle.readline().strip("\r\n").split("\t"),
                            in_handle.readline().strip("\r\n").split("\t")))
    if solution.get("path"):
        out["purity"] = solution["purity"]
        out["ploidy"] = solution["ploidy"]
        out["cellular_prevalence"] = [x.strip() for x in solution["cellPrev"].split(",")]
        base = os.path.basename(solution["path"])
        out["plot"] = dict([(n, solution["path"] + ext) for (n, ext) in [("rplots", ".Rplots.pdf"),
                                                                         ("cf", "/%s_CF.pdf" % base),
                                                                         ("cna", "/%s_CNA.pdf" % base),
                                                                         ("loh", "/%s_LOH.pdf" % base)]
                            if os.path.exists(solution["path"] + ext)])
        out["subclones"] = "%s.segs.txt" % solution["path"]
        out["hetsummary"] = solution_file
        out["vrn_file"] = to_vcf(out["subclones"], "TitanCNA", _get_header, _seg_to_vcf, data)
        out["lohsummary"] = loh.summary_status(out, data)
    return out
1371849
def _prepare_summary(evolve_file, ssm_file, cnv_file, work_dir, somatic_info):
    
    out_file = os.path.join(work_dir, "%s-phylowgs.txt" % somatic_info.tumor_name)
    if not utils.file_uptodate(out_file, evolve_file):
        with file_transaction(somatic_info.tumor_data, out_file) as tx_out_file:
            with open(tx_out_file, "w") as out_handle:
                ssm_locs = _read_ssm_locs(ssm_file)
                cnv_ssms = _read_cnv_ssms(cnv_file)
                for i, (ids, tree) in enumerate(_evolve_reader(evolve_file)):
                    out_handle.write("* Tree %s\n" % (i + 1))
                    out_handle.write("\n" + "\n".join(tree) + "\n\n")
                    for nid, freq, gids in ids:
                        genes = _gids_to_genes(gids, ssm_locs, cnv_ssms, somatic_info.tumor_data)
                        out_handle.write("%s\t%s\t%s\n" % (nid, freq, ",".join(genes)))
                    out_handle.write("\n")
    return out_file
327823
def _organize_output(output, combine_map, file_key, combine_arg_keys):
    
    out_map = collections.defaultdict(list)
    extra_args = collections.defaultdict(list)
    final_args = collections.OrderedDict()
    extras = []
    for data in output:
        cur_file = data.get(file_key)
        if not cur_file:
            extras.append([data])
        else:
            cur_out = combine_map[cur_file]
            out_map[cur_out].append(cur_file)
            extra_args[cur_out].append([data[x] for x in combine_arg_keys])
            data[file_key] = cur_out
            if cur_out not in final_args:
                final_args[cur_out] = [data]
            else:
                extras.append([data])
    combine_args = [[v, k] + _get_extra_args(extra_args[k], combine_arg_keys)
                    for (k, v) in out_map.items()]
    return combine_args, list(final_args.values()) + extras
469384
def create_downsampled_data(params):
    '''
    
    '''
    maxsamples = 1

    for data_type in ['LFP','CSD']:

        if RANK == 0:
            if not os.path.isdir(os.path.join(params.savefolder, 'populations', 'subsamples')):
                os.mkdir((os.path.join(params.savefolder,'populations','subsamples')))
        COMM.Barrier()

        try:
            assert(ana_params.scaling <= params.recordSingleContribFrac)
        except AssertionError as ae:
            raise ae, 'scaling parameter must be less than simulation recordSingleContribFrac'
        
        samples = int(1. / ana_params.scaling)
        if samples > maxsamples:
            samples = maxsamples
        
        COUNTER = 0
        for j, layer in enumerate(params.y_in_Y): # loop over layers
            for k, pop in enumerate(layer): # loop over populations
                for i, y in enumerate(pop): # loop over cell types
                    if COUNTER % SIZE == RANK:
                        # Load data
                        fname = os.path.join(params.savefolder, 'populations', '%s_%ss.h5' \
                                             % (y, data_type))
                        f = h5py.File(fname)
                        print 'Load %s' % str(f.filename)
                        raw_data = f['data'].value
                        srate = f['srate'].value
                        f.close()

                        ## shuffle data
                        #np.random.shuffle(raw_data)

                        # sample size
                        N = int(params.N_y[np.array(params.y) == y]*ana_params.scaling)
                        try:
                            assert(N <= raw_data.shape[0])
                        except AssertionError as ae:
                            raise ae, 'shape mismatch with sample size'


                        for sample in range(samples): # loop over samples

                            # slice data
                            data = raw_data[sample*N:(sample+1)*N]

                            # create cell resolved file
                            fname = os.path.join(params.savefolder,'populations','subsamples',
                                                 '%s_%ss_%i_%i.h5' \
                                                 % (y, data_type, ana_params.scaling*100, sample))
                            f = h5py.File(fname, 'w')
                            print 'Write %s' % str(f.filename)
                            f['data'] = data
                            f['srate'] = srate
                            f.close()

                            # create cell type resolved file
                            fname = os.path.join(params.savefolder,'populations','subsamples',
                                                 '%s_population_%s_%i_%i.h5' \
                                                 % (y,data_type,ana_params.scaling*100,sample))
                            f = h5py.File(fname, 'w')
                            print 'Write %s' % str(f.filename)
                            f['data'] = data.sum(axis=0)
                            f['srate'] = srate
                            f.close()

                    COUNTER += 1

        COMM.Barrier()
        f = h5py.File(os.path.join(params.savefolder,'populations', '%s_%ss.h5' % (y,data_type)), 'r')
        datashape = f['data'].shape
        f.close()

        COUNTER = 0
        for sample in range(samples): # loop over samples
            if COUNTER % SIZE == RANK:
                # initialize full sum signal
                data_full = np.zeros(datashape[1:])
                for j,layer in enumerate(params.y_in_Y): # loop over layers
                    for k,pop in enumerate(layer): # loop over populations

                        # initialize population resolved sum signal
                        data_Y = np.zeros(datashape[1:])

                        for i,y in enumerate(pop): # loop over cell types

                            # Load data
                            fname = os.path.join(params.savefolder, 'populations',
                                                 'subsamples', '%s_population_%s_%i_%i.h5' \
                                                 % (y, data_type, ana_params.scaling*100,
                                                    sample))
                            f = h5py.File(fname, 'r')    
                            # Update population sum:
                            data_Y += f['data'].value
                            srate = f['srate'].value
                            f.close()

                        # write population sum
                        fname = os.path.join(params.savefolder,'populations','subsamples',
                                             '%s_population_%s_%i_%i.h5' \
                                             % (params.Y[2*j+k], data_type,
                                                ana_params.scaling*100, sample))
                        f = h5py.File(fname,'w')
                        print 'Write %s' % str(f.filename)
                        f['data'] = data_Y
                        f['srate'] = srate
                        f.close() 

                        # update full sum
                        data_full += data_Y

                # write sum
                fname = os.path.join(params.savefolder,'populations','subsamples',
                                     '%ssum_%i_%i.h5' % (data_type,ana_params.scaling*100,sample))
                f = h5py.File(fname,'w')
                print 'Write %s' % str(f.filename)
                f['data'] = data_full
                f['srate'] = srate
                f.close()

            COUNTER += 1

        COMM.Barrier()
640589
def output_to_csv(sources, csv_writeable_file):
    """"""
    file_exists = os.path.isfile(csv_writeable_file)

    with open(csv_writeable_file, 'a') as csvfile:
        csv_dict = OrderedDict()
        csv_dict.update({'Time': time.strftime("%Y-%m-%d_%H:%M:%S")})
        summaries = [val for key, val in sources.items()]
        for summarie in summaries:
            csv_dict.update(summarie.source.get_sensors_summary())

        fieldnames = [key for key, val in csv_dict.items()]
        writer = csv.DictWriter(csvfile, fieldnames=fieldnames)

        if not file_exists:
            writer.writeheader()  # file doesn't exist yet, write a header
        writer.writerow(csv_dict)
75010
def generate_data(self, data_dir, tmp_dir=None, task_id=-1):
    """"""
    if not self._rollouts_by_epoch_and_split[self.current_epoch]:
      # Data not loaded from disk.
      self._split_current_epoch()

    rollouts_by_split = self._rollouts_by_epoch_and_split[self.current_epoch]
    splits_and_paths = self.splits_and_paths(data_dir)

    for (split, paths) in splits_and_paths:
      rollouts = rollouts_by_split[split]
      num_frames = self._calc_num_frames(rollouts)
      shard_size = num_frames // len(paths)

      frame_gen = self._generate_frames(rollouts)
      for (path_index, path) in enumerate(paths):
        limit = shard_size
        # Put the remainder in the last shard to preserve the ordering.
        if path_index == len(paths) - 1:
          limit = None
        generator_utils.generate_files(
            itertools.islice(frame_gen, limit), [path],
            cycle_every_n=float("inf")
        )
783759
def create_report(outdirname, report_filename, **kwargs):
    
    # Checking the required variables
    if "steps" in kwargs:
        assert "descriptions" in kwargs
        assert "long_descriptions" in kwargs
        assert "steps_filename" not in kwargs
    else:
        assert "steps_filename" in kwargs
        assert "descriptions" not in kwargs
        assert "long_descriptions" not in kwargs
    assert "summaries" in kwargs
    assert "background" in kwargs
    assert "project_name" in kwargs
    assert "summary_fn" in kwargs
    assert "report_title" in kwargs
    assert "report_author" in kwargs
    assert "initial_files" in kwargs
    assert "final_nb_markers" in kwargs
    assert "final_nb_samples" in kwargs
    assert "final_files" in kwargs
    assert "plink_version" in kwargs
    assert "graphic_paths_fn" in kwargs

    # Formatting the background section
    background_section = _format_background(kwargs["background"])

    # Writing the method steps to a separate file (for access later)
    steps_filename = None
    if "steps_filename" in kwargs:
        steps_filename = kwargs["steps_filename"]
    else:
        steps_filename = os.path.join(outdirname, "steps_summary.tex")
        with open(steps_filename, "w") as o_file:
            zipped = zip(kwargs["steps"], kwargs["descriptions"],
                         kwargs["long_descriptions"])
            for step, desc, long_desc in zipped:
                if desc.endswith("."):
                    desc = desc[:-1]
                step = step.replace("_", r"\_")
                to_print = latex.item(desc)
                to_print += " [{}].".format(latex.texttt(step))
                if long_desc is not None:
                    to_print += " " + long_desc
                print >>o_file, latex.wrap_lines(to_print) + "\n"

    # Adding the content of the results section
    result_summaries = []
    for name in kwargs["summaries"]:
        full_path = os.path.abspath(name)
        if os.path.isfile(full_path):
            rel_path = os.path.relpath(full_path, outdirname)
            result_summaries.append(re.sub(r"\\", "/", rel_path))

    # Reading the initial_files file
    initial_files = None
    with open(kwargs["initial_files"], "r") as i_file:
        initial_files = i_file.read().splitlines()

    # Reading the final_files file
    final_files = None
    with open(kwargs["final_files"], "r") as i_file:
        final_files = [i.split("\t")[0] for i in i_file.read().splitlines()]

    # Adding the bibliography content
    biblio_entry = latex.bib_entry(
        name="pyGenClean",
        authors="Lemieux Perreault LP, Provost S, Legault MA, Barhdadi A, "
                r"Dub\'e MP",
        title="pyGenClean: efficient tool for genetic data clean up before "
              "association testing",
        journal="Bioinformatics",
        year="2013",
        volume="29",
        number="13",
        pages="1704--1705",
    ) + "\n" * 2 + latex.bib_entry(
        name="plink",
        authors="Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, "
                "Bender D, Maller J, Sklar P, de Bakker PIW, Daly MJ, Sham PC",
        title="PLINK: a tool set for whole-genome association and "
              "population-based linkage analyses",
        journal="American Journal of Human Genetics",
        year="2007",
        volume="81",
        number="3",
        pages="559--575",
    ) + "\n" * 2 + latex.bib_entry(
        name="bafRegress",
        authors=r"Goo J, Matthew F, Kurt NH, Jane MR, Kimberly FD, "
                r"Gon{\c{c}}alo RA, Michael B, Hyun Min K",
        title="Detecting and estimating contamination of human DNA samples in "
              "sequencing and array-based genotype data",
        journal="The American Journal of Human Genetics",
        year="2012",
        volume="91",
        number="5",
        pages="839--848",
    )

    # Getting the template
    main_template = latex.jinja2_env.get_template("main_document.tex")

    # Getting the data
    today = datetime.today()

    # Reading the graphics path
    graphic_paths = []
    if kwargs["graphic_paths_fn"] is not None:
        with open(kwargs["graphic_paths_fn"], "r") as i_file:
            graphic_paths = [
                re.sub(r"\\", "/", path) + ("" if path.endswith("/") else "/")
                for path in i_file.read().splitlines()
            ]

    try:
        with open(report_filename, "w") as i_file:
            # Rendering the template
            print >>i_file, main_template.render(
                project_name=latex.sanitize_tex(kwargs["project_name"]),
                month=today.strftime("%B"),
                day=today.day,
                year=today.year,
                background_content=background_section,
                result_summaries=result_summaries,
                bibliography_content=biblio_entry,
                pygenclean_version=pygenclean_version,
                plink_version=kwargs["plink_version"],
                steps_filename=os.path.basename(steps_filename),
                final_results=_create_summary_table(
                    kwargs["summary_fn"],
                    latex.jinja2_env.get_template("summary_table.tex"),
                    nb_samples=kwargs["final_nb_samples"],
                    nb_markers=kwargs["final_nb_markers"],
                ),
                report_title=latex.sanitize_tex(kwargs["report_title"]),
                report_author=latex.sanitize_tex(kwargs["report_author"]),
                initial_files=initial_files,
                final_files=final_files,
                final_nb_samples=kwargs["final_nb_samples"],
                final_nb_markers=kwargs["final_nb_markers"],
                graphic_paths=graphic_paths,
            )

    except IOError:
        msg = "{}: could not create report".format(report_filename)
        raise ProgramError(msg)
1030031
def merge(args):
    
    from jcvi.formats.gff import merge as gmerge

    p = OptionParser(merge.__doc__)
    p.set_home("maker")
    opts, args = p.parse_args(args)

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

    outdir, outputgff = args
    fsnames, suffix = get_fsnames(outdir)
    nfs = len(fsnames)
    cmd = op.join(opts.maker_home, "bin/gff3_merge")

    outfile = "merge.sh"
    write_file(outfile, mergesh.format(suffix, cmd))

    # Generate per split directory
    # Note that gff3_merge write to /tmp, so I limit processes here to avoid
    # filling up disk space
    sh("parallel -j 8 merge.sh {} ::: " + " ".join(fsnames))

    # One final output
    gffnames = glob("*.all.gff")
    assert len(gffnames) == nfs

    # Again, DO NOT USE gff3_merge to merge with a smallish /tmp/ area
    gfflist = "gfflist"
    fw = open(gfflist, "w")
    print("\n".join(gffnames), file=fw)
    fw.close()

    nlines = sum(1 for x in open(gfflist))
    assert nlines == nfs  # Be extra, extra careful to include all results
    gmerge([gfflist, "-o", outputgff])
    logging.debug("Merged GFF file written to `{0}`".format(outputgff))
730870
def _prepare_output_multi(self, model):
        """"""
        model_name = model.__name__
        current_path = os.path.join(self._output_path, '{model}.{extension}'.format(
            model=model_name,
            extension=self.EXTENSION,
        ))
        self._outfile = codecs.open(current_path, 'w', encoding='utf-8')
        print('Dumping {model} to {file}'.format(model=model_name, file=current_path))
569094
def _add_segmetrics_to_output(out, data):
    
    out_file = "%s-segmetrics.txt" % os.path.splitext(out["cns"])[0]
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            cmd = [os.path.join(os.path.dirname(sys.executable), "cnvkit.py"), "segmetrics",
                   "--median", "--iqr", "--ci", "--pi",
                   "-s", out["cns"], "-o", tx_out_file, out["cnr"]]
            # Use less fine grained bootstrapping intervals for whole genome runs
            if dd.get_coverage_interval(data) == "genome":
                cmd += ["--alpha", "0.1", "--bootstrap", "50"]
            else:
                cmd += ["--alpha", "0.01", "--bootstrap", "500"]
            do.run(cmd, "CNVkit segmetrics")
    out["segmetrics"] = out_file
    return out
1244698
def write_output(self):
        """"""
        for data in self.output_data.values():
            self.create_output(data.get('key'), data.get('value'), data.get('type'))
628274
def batchcn(args):
    
    p = OptionParser(batchcn.__doc__)
    p.add_option("--upload", default="s3://hli-mv-data-science/htang/ccn",
                 help="Upload cn and seg results to s3")
    opts, args = p.parse_args(args)

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

    workdir, samples = args
    upload = opts.upload
    store = upload + "/{}/*.seg".format(workdir)
    computed = [op.basename(x).split(".")[0] for x in glob_s3(store)]
    computed = set(computed)

    # Generate a bunch of cn commands
    fp = open(samples)
    nskipped = ntotal = 0
    cmd = "python -m jcvi.variation.cnv cn --hmm --cleanup {}".format(workdir)
    for row in fp:
        samplekey, path = row.strip().split(",")
        ntotal += 1
        if samplekey in computed:
            nskipped += 1
            continue
        print(" ".join((cmd, samplekey, path)))

    logging.debug("Skipped: {}".format(percentage(nskipped, ntotal)))
897338
def writeMNIST(sc, input_images, input_labels, output, format, num_partitions):
  """"""
  # load MNIST gzip into memory
  with open(input_images, 'rb') as f:
    images = numpy.array(mnist.extract_images(f))

  with open(input_labels, 'rb') as f:
    if format == "csv2":
      labels = numpy.array(mnist.extract_labels(f, one_hot=False))
    else:
      labels = numpy.array(mnist.extract_labels(f, one_hot=True))

  shape = images.shape
  print("images.shape: {0}".format(shape))          # 60000 x 28 x 28
  print("labels.shape: {0}".format(labels.shape))   # 60000 x 10

  # create RDDs of vectors
  imageRDD = sc.parallelize(images.reshape(shape[0], shape[1] * shape[2]), num_partitions)
  labelRDD = sc.parallelize(labels, num_partitions)

  output_images = output + "/images"
  output_labels = output + "/labels"

  # save RDDs as specific format
  if format == "pickle":
    imageRDD.saveAsPickleFile(output_images)
    labelRDD.saveAsPickleFile(output_labels)
  elif format == "csv":
    imageRDD.map(toCSV).saveAsTextFile(output_images)
    labelRDD.map(toCSV).saveAsTextFile(output_labels)
  elif format == "csv2":
    imageRDD.map(toCSV).zip(labelRDD).map(lambda x: str(x[1]) + "|" + x[0]).saveAsTextFile(output)
  else:  # format == "tfr":
    tfRDD = imageRDD.zip(labelRDD).map(lambda x: (bytearray(toTFExample(x[0], x[1])), None))
    # requires: --jars tensorflow-hadoop-1.0-SNAPSHOT.jar
    tfRDD.saveAsNewAPIHadoopFile(output, "org.tensorflow.hadoop.io.TFRecordFileOutputFormat",
                                 keyClass="org.apache.hadoop.io.BytesWritable",
                                 valueClass="org.apache.hadoop.io.NullWritable")
284685
def _do_merge(orig_files, out_file, config, region):
    
    if not utils.file_exists(out_file):
        with file_transaction(config, out_file) as tx_out_file:
            _check_samples_nodups(orig_files)
            prep_files = run_multicore(p_bgzip_and_index, [[x, config] for x in orig_files], config)
            input_vcf_file = "%s-files.txt" % utils.splitext_plus(out_file)[0]
            with open(input_vcf_file, "w") as out_handle:
                for fname in prep_files:
                    out_handle.write(fname + "\n")
            bcftools = config_utils.get_program("bcftools", config)
            output_type = "z" if out_file.endswith(".gz") else "v"
            region_str = "-r {}".format(region) if region else ""
            cmd = "{bcftools} merge -O {output_type} {region_str} `cat {input_vcf_file}` > {tx_out_file}"
            do.run(cmd.format(**locals()), "Merge variants")
    if out_file.endswith(".gz"):
        bgzip_and_index(out_file, config)
    return out_file
146333
def write_output(self, data, args=None, filename=None, label=None):
        """"""
        if args:
            if not args.outlog:
                return 0
        if not filename: filename=args.outlog
        lastpath = ''
        with open(str(filename), 'w') as output_file:
            for entry in data['entries']:
                if args.label:
                    if entry['source_path'] == lastpath:
                        output_file.write(entry['raw_text'] + '\n')
                    elif args.label == 'fname':
                        output_file.write('======== ' + \
                                entry['source_path'].split('/')[-1] + \
                                ' >>>>\n' + entry['raw_text'] + '\n')
                    elif args.label == 'fpath':
                        output_file.write('======== ' + \
                                entry['source_path']  + \
                                ' >>>>\n' + entry['raw_text'] + '\n')
                else: output_file.write(entry['raw_text'] + '\n')
                lastpath = entry['source_path']
502809
def combine_pyramid_and_save(g_video, orig_video, enlarge_multiple, fps, save_filename='media/output.avi'):
    """"""
    width, height = get_frame_dimensions(orig_video[0])
    fourcc = cv2.VideoWriter_fourcc(*'MJPG')
    print("Outputting to %s" % save_filename)
    writer = cv2.VideoWriter(save_filename, fourcc, fps, (width, height), 1)
    for x in range(0, g_video.shape[0]):
        img = np.ndarray(shape=g_video[x].shape, dtype='float')
        img[:] = g_video[x]
        for i in range(enlarge_multiple):
            img = cv2.pyrUp(img)

        img[:height, :width] = img[:height, :width] + orig_video[x]
        res = cv2.convertScaleAbs(img[:height, :width])
        writer.write(res)