%prog compare NA12878_array_hg38.bed *.seg

    Compare cnv output to known ground truths.
497169
def compare(args):
    
    p = OptionParser(compare.__doc__)
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    truths = args[0]
    cnvoutputs = args[1:]
    cpus = min(len(cnvoutputs), opts.cpus)
    p = Pool(processes=cpus)
    results = []
    files = [(x, truths) for x in cnvoutputs]
    r = p.map_async(compare_worker, files, callback=results.append)
    r.wait()

    for res in results:
        print("\n".join(res))
246439
def check_bed_coords(in_file, data):
    
    if dd.get_ref_file(data):
        contig_sizes = {}
        for contig in ref.file_contigs(dd.get_ref_file(data)):
            contig_sizes[contig.name] = contig.size
        with utils.open_gzipsafe(in_file) as in_handle:
            for line in in_handle:
                if not line.startswith(("#", "track", "browser", "@")) and line.strip():
                    parts = line.split()
                    if len(parts) > 3:
                        try:
                            end = int(parts[2])
                        except ValueError:
                            continue
                        contig = parts[0]
                        check_size = contig_sizes.get(contig)
                        if check_size and end > check_size:
                            raise ValueError("Found BED coordinate off the end of the chromosome:\n%s%s\n"
                                             "Is the input BED from the right genome build?" %
                                             (line, in_file))
648438
def validate(args):
    
    p = OptionParser(validate.__doc__)
    opts, args = p.parse_args(args)

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

    imputed, withheld = args
    register = {}
    fp = open(withheld)
    for row in fp:
        if row[0] == "#":
            continue
        v = VcfLine(row)
        register[(v.seqid, v.pos)] = v.genotype

    logging.debug("Imported {0} records from `{1}`".\
                    format(len(register), withheld))

    fp = must_open(imputed)
    hit = concordant = 0
    seen = set()
    for row in fp:
        if row[0] == "#":
            continue
        v = VcfLine(row)
        chr, pos, genotype = v.seqid, v.pos, v.genotype
        if (chr, pos) in seen:
            continue
        seen.add((chr, pos))
        if (chr, pos) not in register:
            continue
        truth = register[(chr, pos)]
        imputed = genotype.split(":")[0]
        if "|" in imputed:
            imputed = "/".join(sorted(genotype.split(":")[0].split("|")))
            #probs = [float(x) for x in genotype.split(":")[-1].split(",")]
            #imputed = max(zip(probs, ["0/0", "0/1", "1/1"]))[-1]
        hit += 1
        if truth == imputed:
            concordant += 1
        else:
            print(row.strip(), "truth={0}".format(truth), file=sys.stderr)

    logging.debug("Total concordant: {0}".\
            format(percentage(concordant, hit)))
1350545
def main (args):
  
  if len(args) != 3:
    print(main.__doc__)
    sys.exit(1)

  d = libsedml.readSedML(args[1]);
  if ( d.getErrorLog().getNumFailsWithSeverity(libsedml.LIBSEDML_SEV_ERROR) > 0):
    print (d.getErrorLog().toString());
  else:
    libsedml.writeSedML(d, args[2]);

  return 0;
916906
def diff(args):
    
    from jcvi.utils.table import banner

    p = OptionParser(diff.__doc__)
    p.add_option("--ignore_case", default=False, action="store_true",
            help="ignore case when comparing sequences [default: %default]")
    p.add_option("--ignore_N", default=False, action="store_true",
            help="ignore N and X's when comparing sequences [default: %default]")
    p.add_option("--ignore_stop", default=False, action="store_true",
            help="ignore stop codon when comparing sequences [default: %default]")
    p.add_option("--rc", default=False, action="store_true",
            help="also consider reverse complement [default: %default]")
    p.add_option("--quiet", default=False, action="store_true",
            help="don't output comparison details [default: %default]")

    opts, args = p.parse_args(args)

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

    afasta, bfasta = args

    afastan = len(Fasta(afasta))
    bfastan = len(Fasta(bfasta))

    if afastan == bfastan:
        print(green("Two sets contain the same number of sequences ({0}, {1})".\
                format(afastan, bfastan)))
    else:
        print(red("Two sets contain different number of sequences ({0}, {1})".\
                format(afastan, bfastan)))

    ah = SeqIO.parse(afasta, "fasta")
    bh = SeqIO.parse(bfasta, "fasta")

    problem_ids = []
    for arec, brec in zip(ah, bh):

        if opts.ignore_stop:
            arec.seq = arec.seq.rstrip("*")
            brec.seq = brec.seq.rstrip("*")

        asize, bsize = len(arec), len(brec)

        if not opts.quiet:
            print(banner(str(arec), [str(brec)]))
            if asize == bsize:
                print(green("Two sequence size match (%d)" % asize))
            else:
                print(red("Two sequence size do not match (%d, %d)" % (asize, bsize)))

        # print out the first place the two sequences diff
        fd = print_first_difference(arec, brec, ignore_case=opts.ignore_case,
                ignore_N=opts.ignore_N, rc=opts.rc, report_match=not opts.quiet)
        if not fd:
            logging.error("Two sets of sequences differ at `{0}`".format(arec.id))
            problem_ids.append("\t".join(str(x) for x in (arec.id, asize, bsize,
                    abs(asize - bsize))))

    if problem_ids:
        print(red("A total of {0} records mismatch.".format(len(problem_ids))))
        fw = must_open("Problems.ids", "w")
        print("\n".join(problem_ids), file=fw)
1233175
def _add_seg_to_output(out, data, enumerate_chroms=False):
    
    out_file = "%s.seg" % 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"), "export",
                   "seg"]
            if enumerate_chroms:
                cmd += ["--enumerate-chroms"]
            cmd += ["-o", tx_out_file, out["cns"]]
            do.run(cmd, "CNVkit export seg")
    out["seg"] = out_file
    return out
1122384
def _check_sig_dims(sig, read_len, n_sig, samps_per_frame):
    
    if isinstance(sig, np.ndarray):
        if sig.shape != (read_len, n_sig):
            raise ValueError('Samples were not loaded correctly')
    else:
        if len(sig) != n_sig:
            raise ValueError('Samples were not loaded correctly')
        for ch in range(n_sig):
            if len(sig[ch]) != samps_per_frame[ch] * read_len:
                raise ValueError('Samples were not loaded correctly')
1046376
def benchmark(args):
    
    p = OptionParser(benchmark.__doc__)
    opts, args = p.parse_args(args)

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

    pf, bedfile = args
    truth = pf + ".truth"
    synfind = pf + ".synfind"
    mcscanx = pf + ".mcscanx"
    iadhore = pf + ".iadhore"
    orthofinder = pf + ".orthofinder"
    pivots = set([x.accn for x in Bed(bedfile)])

    fp = open(truth)
    truth = set()
    for row in fp:
        a, b = row.strip().split("\t")[:2]
        pivots.add(a)
        truth.add(tuple(sorted((a, b))))
    logging.debug("Truth: {0} pairs".format(len(truth)))

    fp = open(synfind)
    benchmarkfile = pf + ".benchmark"
    fw = must_open(benchmarkfile, "w")
    synfind = set()
    for row in fp:
        atoms = row.strip().split("\t")
        query, hit, tag = atoms[:3]
        if tag != "S":
            continue
        synfind.add(tuple(sorted((query, hit))))
    calc_sensitivity_specificity(synfind, truth, "SynFind", fw)

    fp = open(mcscanx)
    mcscanx = set()
    for row in fp:
        if row[0] == '#':
            continue
        atoms = row.strip().split(":")[1].split()
        query, hit = atoms[:2]
        mcscanx.add(tuple(sorted((query, hit))))
    calc_sensitivity_specificity(mcscanx, truth, "MCScanX", fw)

    fp = open(iadhore)
    iadhore = set()
    next(fp)
    for row in fp:
        atoms = row.strip().split("\t")
        query, hit = atoms[3:5]
        iadhore.add(tuple(sorted((query, hit))))
    calc_sensitivity_specificity(iadhore, truth, "iADHoRe", fw)

    fp = open(orthofinder)
    orthofinder = set()
    next(fp)
    for row in fp:
        row = row.replace('"', "")
        atoms = row.replace(",", " ").split()
        genes = [x.strip() for x in atoms if not x.startswith("OG")]
        genes = [gene_name(x) for x in genes]
        pps = [x for x in genes if x in pivots]
        for p in pps:
            for g in genes:
                if p == g:
                    continue
                orthofinder.add(tuple(sorted((p, g))))
    #write_pairs(orthofinder, "orthofinder.pairs")
    calc_sensitivity_specificity(orthofinder, truth, "OrthoFinder", fw)
    fw.close()
1108225
def find_all_output_in_range(self, ifo, currSeg, useSplitLists=False):
        
        if not useSplitLists:
            # Slower, but simpler method
            outFiles = [i for i in self if ifo in i.ifo_list]
            outFiles = [i for i in outFiles \
                                      if i.segment_list.intersects_segment(currSeg)]
        else:
            # Faster, but more complicated
            # Basically only check if a subset of files intersects_segment by
            # using a presorted list. Sorting only happens once.
            if not self._check_split_list_validity():
                # FIXME: DO NOT hard code this.
                self._temporal_split_list(100)
            startIdx = int( (currSeg[0] - self._splitListsStart) / \
                                                          self._splitListsStep )
            # Add some small rounding here
            endIdx = (currSeg[1] - self._splitListsStart) / self._splitListsStep
            endIdx = int(endIdx - 0.000001)

            outFiles = []
            for idx in range(startIdx, endIdx + 1):
                if idx < 0 or idx >= self._splitListsNum:
                    continue
                outFilesTemp = [i for i in self._splitLists[idx] \
                                                            if ifo in i.ifo_list]
                outFiles.extend([i for i in outFilesTemp \
                                      if i.segment_list.intersects_segment(currSeg)])
                # Remove duplicates
                outFiles = list(set(outFiles))

        return self.__class__(outFiles)
650742
def assert_output(output, assert_equal):
    
    sorted_output = sorted(output)
    sorted_assert = sorted(assert_equal)
    if sorted_output != sorted_assert:
        raise ValueError(ASSERT_ERROR.format(sorted_output, sorted_assert))
381158
def checkIsConsistent(self):
        
        if is_an_array(self.mask) and self.mask.shape != self.data.shape:
            raise ConsistencyError("Shape mismatch mask={}, data={}"
                                   .format(self.mask.shape != self.data.shape))
893440
def dump(self):
        
        for seg in self.seglist:
            print("==== %08x-%08x" % (seg.startea, seg.endea))
            if seg.endea - seg.startea < 30:
                for ea in range(seg.startea, seg.endea):
                    print("    %08x: %08x" % (ea, self.getFlags(ea)))
            else:
                for ea in range(seg.startea, seg.startea + 10):
                    print("    %08x: %08x" % (ea, self.getFlags(ea)))
                print("...")
                for ea in range(seg.endea - 10, seg.endea):
                    print("    %08x: %08x" % (ea, self.getFlags(ea)))
1187584
def infer(args):
    
    from jcvi.apps.grid import WriteJobs
    from jcvi.formats.bed import sort

    p = OptionParser(infer.__doc__)
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    scaffoldsf, genomef = args
    inferbed = "infer-components.bed"
    if need_update((scaffoldsf, genomef), inferbed):
        scaffolds = Fasta(scaffoldsf, lazy=True)
        genome = Fasta(genomef)
        genome = genome.tostring()
        args = [(scaffold_name, scaffold, genome) \
                for scaffold_name, scaffold in scaffolds.iteritems_ordered()]

        pool = WriteJobs(map_one_scaffold, args, inferbed, cpus=opts.cpus)
        pool.run()

    sort([inferbed, "-i"])
    bed = Bed(inferbed)
    inferagpbed = "infer.bed"
    fw = open(inferagpbed, "w")
    seen = []
    for b in bed:
        r = (b.seqid, b.start, b.end)
        if check_seen(r, seen):
            continue
        print("\t".join(str(x) for x in \
                    (b.accn, 0, b.span, b.seqid, b.score, b.strand)), file=fw)
        seen.append(r)
    fw.close()

    frombed([inferagpbed])
903804
function _compare (hay, needle) {
    for(var i in set) {
      var x = u.get(hay, set[i]), y = needle[i]
      if(x !== y) return compare(x, y)
    }
    return 0
  }
640086
def check_ordered(self): 
    
    sys.stderr.write("error unimplemented check_ordered\n")
    sys.exit()
    seen_chrs = set()
    curr_chr = None
    prevstart = 0
    for l in self._lines:
      if not l['rng']: continue
      if l['rng'].chr != curr_chr:
        prevstart = 0
        if l['rng'].chr in seen_chrs:
          return False
        curr_chr = l['rng'].chr
        seen_chrs.add(curr_chr)
      if l['rng'].start < prevstart:  return False
      prevstart = l['rng'].start
    return True
845384
def main(argv=None):
    

    if argv is None:
        argv = sys.argv

    # setup command line parser
    parser = U.OptionParser(version="%prog version: $Id$",
                            usage=usage,
                            description=globals()["__doc__"])

    group = U.OptionGroup(parser, "extract-specific options")

    # (Experimental option) Retain the UMI in the sequence read"
    group.add_option("--retain-umi", dest="retain_umi", action="store_true",
                     help=optparse.SUPPRESS_HELP)
    group.add_option("--read2-out", dest="read2_out", type="string",
                     help="file to output processed paired read to")
    group.add_option("--read2-stdout", dest="read2_stdout",
                     action="store_true",
                     help="Paired reads, send read2 to stdout, discarding read1")
    group.add_option("--quality-filter-threshold",
                     dest="quality_filter_threshold", type="int",
                     help=("Remove reads where any UMI base quality score "
                           "falls below this threshold"))
    group.add_option("--quality-filter-mask",
                     dest="quality_filter_mask", type="int",
                     help=("If a UMI base has a quality below this threshold, "
                           "replace the base with 'N'"))
    group.add_option("--quality-encoding",
                     dest="quality_encoding", type="choice",
                     choices=["phred33", "phred64", "solexa"],
                     help=("Quality score encoding. Choose from 'phred33'"
                           "[33-77] 'phred64' [64-106] or 'solexa' [59-106]"))
    group.add_option("--filter-cell-barcode",
                     dest="filter_cell_barcode",
                     action="store_true",
                     help="Filter the cell barcodes")
    group.add_option("--error-correct-cell",
                     dest="error_correct_cell",
                     action="store_true",
                     help=("Correct errors in the cell barcode"))
    group.add_option("--whitelist",
                     dest="whitelist", type="string",
                     help=("A whitelist of accepted cell barcodes"))
    group.add_option("--blacklist",
                     dest="blacklist", type="string",
                     help=("A blacklist of rejected cell barcodes"))
    group.add_option("--subset-reads", "--reads-subset",
                     dest="reads_subset", type="int",
                     help=("Only extract from the first N reads. If N is "
                           "greater than the number of reads, all reads will "
                           "be used"))
    group.add_option("--reconcile-pairs",
                     dest="reconcile", action="store_true",
                     help=("Allow the presences of reads in read2 input that "
                           "are not present in read1 input. This allows cell "
                           "barcode filtering of read1s without "
                           "considering read2s"))

    parser.add_option_group(group)

    group = U.OptionGroup(parser, "[EXPERIMENTAl] barcode extraction options")

    group.add_option("--either-read", dest="either_read", action="store_true",
                     help="UMI may be on either read (see "
                     "--either-read-resolve) for options to resolve cases where"
                     "UMI is on both reads")
    group.add_option("--either-read-resolve",
                     dest="either_read_resolve", type="choice",
                     choices=["discard", "quality"],
                     help=("How to resolve instances where both reads "
                           "contain a UMI but using --either-read."
                           "Choose from 'discard' or 'quality'"
                           "(use highest quality). default=dicard"))

    parser.add_option_group(group)

    parser.set_defaults(extract_method="string",
                        filter_cell_barcodes=False,
                        whitelist=None,
                        blacklist=None,
                        error_correct_cell=False,
                        pattern=None,
                        pattern2=None,
                        read2_in=None,
                        read2_out=False,
                        read2_stdout=False,
                        quality_filter_threshold=None,
                        quality_encoding=None,
                        reconcile=False,
                        either_read=False,
                        either_read_resolve="discard")

    # add common options (-h/--help, ...) and parse command line

    (options, args) = U.Start(parser, argv=argv,
                              add_extract_options=True,
                              add_group_dedup_options=False,
                              add_umi_grouping_options=False,
                              add_sam_options=False)

    if options.retain_umi and not options.extract_method == "regex":
        raise ValueError("option --retain-umi only works with --extract-method=regex")

    if options.quality_filter_threshold or options.quality_filter_mask:
        if not options.quality_encoding:
            U.error("must provide a quality encoding (--quality-"
                    "encoding) to filter UMIs by quality (--quality"
                    "-filter-threshold) or mask low quality bases "
                    "with (--quality-filter-mask)")

    extract_cell, extract_umi = U.validateExtractOptions(options)

    if options.either_read:
        if extract_cell:
            U.error("Option to extract from either read (--either-read) "
                    "is not currently compatible with cell barcode extraction")
        if not options.extract_method == "regex":
            U.error("Option to extract from either read (--either-read)"
                    "requires --extract-method=regex")
        if not options.pattern or not options.pattern2:
            U.error("Option to extract from either read (--either-read)"
                    "requires --bc-pattern=[PATTERN1] and"
                    "--bc-pattern2=[PATTERN2]")

    if options.filter_cell_barcodes:

        if not options.whitelist:
                U.error("must provide a whitelist (--whitelist) if using "
                        "--filter-cell-barcode option")

        if not extract_cell:
            if options.extract_method == "string":
                U.error("barcode pattern(s) do not include any cell bases "
                        "(marked with 'Cs') %s, %s" % (
                            options.pattern, options.pattern2))
            elif options.extract_method == "regex":
                U.error("barcode regex(es) do not include any cell groups "
                        "(starting with 'cell_') %s, %s" (
                            options.pattern, options.pattern2))

    read1s = umi_methods.fastqIterate(options.stdin)

    # set up read extractor
    ReadExtractor = extract_methods.ExtractFilterAndUpdate(
        options.extract_method,
        options.pattern,
        options.pattern2,
        options.prime3,
        extract_cell,
        options.quality_encoding,
        options.quality_filter_threshold,
        options.quality_filter_mask,
        options.filter_cell_barcode,
        options.retain_umi,
        options.either_read,
        options.either_read_resolve)

    if options.filter_cell_barcode:
        cell_whitelist, false_to_true_map = whitelist_methods.getUserDefinedBarcodes(
            options.whitelist, options.error_correct_cell)

        ReadExtractor.cell_whitelist = cell_whitelist
        ReadExtractor.false_to_true_map = false_to_true_map

    if options.blacklist:
        blacklist = set()
        with U.openFile(options.blacklist, "r") as inf:
            for line in inf:
                blacklist.add(line.strip().split("\t")[0])
        ReadExtractor.cell_blacklist = blacklist

    # variables for progress monitor
    progCount = 0
    displayMax = 100000
    U.info("Starting barcode extraction")

    if options.read2_in is None:
        for read in read1s:

            # incrementing count for monitoring progress
            progCount += 1

            # Update display in every 100kth iteration
            if progCount % displayMax == 0:
                U.info("Parsed {} reads".format(progCount))

            new_read = ReadExtractor(read)

            if options.reads_subset:
                if (ReadExtractor.read_counts['Input Reads'] >
                    options.reads_subset):
                    break

            if not new_read:
                continue

            options.stdout.write(str(new_read) + "\n")

    else:
        read2s = umi_methods.fastqIterate(U.openFile(options.read2_in))

        if options.read2_out:
            read2_out = U.openFile(options.read2_out, "w")

        if options.reconcile:
            strict = False
        else:
            strict = True

        for read1, read2 in umi_methods.joinedFastqIterate(
                read1s, read2s, strict):

            # incrementing count for monitoring progress
            progCount += 1

            # Update display in every 100kth iteration
            if progCount % displayMax == 0:
                U.info("Parsed {} reads".format(progCount))
                sys.stdout.flush()

            reads = ReadExtractor(read1, read2)

            if options.reads_subset:
                if (ReadExtractor.read_counts['Input Reads'] >
                    options.reads_subset):
                    break

            if not reads:
                continue
            else:
                new_read1, new_read2 = reads

            if options.read2_stdout:
                options.stdout.write(str(new_read2) + "\n")
            else:
                options.stdout.write(str(new_read1) + "\n")

                if options.read2_out:
                    read2_out.write(str(new_read2) + "\n")

    if options.read2_out:
        read2_out.close()

    for k, v in ReadExtractor.getReadCounts().most_common():
        U.info("%s: %s" % (k, v))

    U.Stop()
1269185
def find_synonymous(input_file, work_dir):
    
    cwd = os.getcwd()
    os.chdir(work_dir)
    # create the .ctl file
    ctl_file = "yn-input.ctl"
    output_file = "nuc-subs.yn"
    ctl_h = open(ctl_file, "w")
    ctl_h.write("seqfile = %s\noutfile = %s\nverbose = 0\n" %
                (op.basename(input_file), output_file))
    ctl_h.write("icode = 0\nweighting = 0\ncommonf3x4 = 0\n")
    ctl_h.close()

    cl = YnCommandline(ctl_file)
    print("\tyn00:", cl, file=sys.stderr)
    r, e = cl.run()
    ds_value_yn = None
    ds_value_ng = None
    dn_value_yn = None
    dn_value_ng = None

    # Nei-Gojobori
    output_h = open(output_file)
    row = output_h.readline()
    while row:
        if row.find("Nei & Gojobori") >=0:
            for x in xrange(5):
                row = next(output_h)
            dn_value_ng, ds_value_ng = row.split('(')[1].split(')')[0].split()
            break
        row = output_h.readline()
    output_h.close()

    # Yang
    output_h = open(output_file)
    for line in output_h:
        if line.find("+-") >= 0 and line.find("dS") == -1:
            parts = line.split(" +-")
            ds_value_yn = extract_subs_value(parts[1])
            dn_value_yn = extract_subs_value(parts[0])

    if ds_value_yn is None or ds_value_ng is None:
        h = open(output_file)
        print("yn00 didn't work: \n%s" % h.read(), file=sys.stderr)

    os.chdir(cwd)
    return ds_value_yn, dn_value_yn, ds_value_ng, dn_value_ng