%prog phase genbankfiles

    Input has to be gb file. Search the `KEYWORDS` section to look for PHASE.
    Also look for "chromosome" and "clone" in the definition line.
354243
def phase(args):
    
    p = OptionParser(phase.__doc__)
    p.set_outfile()

    opts, args = p.parse_args(args)

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

    fw = must_open(opts.outfile, "w")
    for gbfile in args:
        for rec in SeqIO.parse(gbfile, "gb"):
            bac_phase, keywords = get_phase(rec)
            chr, clone = get_clone(rec)
            keyword_field = ";".join(keywords)
            print("\t".join((rec.id, str(bac_phase), keyword_field,
                    chr, clone)), file=fw)
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)])
723000
def bed(args):
    
    p = OptionParser(bed.__doc__)
    opts, args = p.parse_args(args)

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

    idsfile, = args
    ids = set(x.strip() for x in open(idsfile))
    data = get_bed_from_phytozome(list(ids))

    pf = idsfile.rsplit(".", 1)[0]
    bedfile = pf + ".bed"
    fw = open(bedfile, "w")
    for i, row in enumerate(data):
        row = row.strip()
        if row == "":
            continue

        print(row, file=fw)

    logging.debug("A total of {0} records written to `{1}`.".format(i + 1, bedfile))
1199108
def _stage_input_files(self, file_mapping, dry_run=True):
        """"""
        # print ("Staging input ", file_mapping)
        if self._file_stage is None:
            return
        self._file_stage.copy_to_scratch(file_mapping, dry_run)
1328645
def run_genotyper(bam_file, ref_file, dbsnp_file, config_file):
    
    cl = ["gatk_genotyper.py", config_file, ref_file, bam_file]
    if dbsnp_file:
        cl.append(dbsnp_file)
    subprocess.check_call(cl)
380996
def check_denovo_input(inputfile, params):
    
    background = params["background"]
    
    input_type = determine_file_type(inputfile)
    
    if input_type == "fasta":
        valid_bg = FA_VALID_BGS    
    elif input_type in ["bed", "narrowpeak"]:
        genome = params["genome"]
        valid_bg = BED_VALID_BGS    
        if "genomic" in background or "gc" in background:
            Genome(genome)
        # is it a valid bed-file etc.
        check_bed_file(inputfile)    # bed-specific, will also work for narrowPeak
    else:
        sys.stderr.write("Format of inputfile {} not recognized.\n".format(inputfile))
        sys.stderr.write("Input should be FASTA, BED or narrowPeak.\n")
        sys.stderr.write("See https://genome.ucsc.edu/FAQ/FAQformat.html for specifications.\n")
        sys.exit(1)

    for bg in background:
        if not bg in valid_bg:
            logger.info("Input type is %s, ignoring background type '%s'", 
                            input_type, bg)
        background = [bg for bg in background if bg in valid_bg]

    if len(background) == 0:
        logger.error("No valid backgrounds specified!")
        sys.exit(1)

    return input_type, background
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))
698550
private Input processGA4GHInput(String input) throws IOException, GeneralSecurityException, URISyntaxException {
    GA4GHUrl url = new GA4GHUrl(input);
    SAMFilePump pump;
    if (usingGrpc) {
      factoryGrpc.configure(url.getRootUrl(), 
        new Settings(clientSecretsFilename, apiKey, noLocalServer));
      pump = new ReadIteratorToSAMFilePump<
          com.google.genomics.v1.Read,
          com.google.genomics.v1.ReadGroupSet, 
          com.google.genomics.v1.Reference>(
              factoryGrpc
                .get(url.getRootUrl())
                .getReads(url));
    } else {
      factoryRest.configure(url.getRootUrl(), 
          new Settings(clientSecretsFilename, apiKey, noLocalServer));
        pump = new ReadIteratorToSAMFilePump<
            com.google.api.services.genomics.model.Read,
            com.google.api.services.genomics.model.ReadGroupSet, 
            com.google.api.services.genomics.model.Reference>(
                factoryRest
                  .get(url.getRootUrl())
                  .getReads(url)); 
    }
    return new Input(input, STDIN_FILE_NAME, pump);
  }
636235
def parse_cmdline():
    """"""
    parser = ArgumentParser(prog="genbank_get_genomes_by_taxon.py")
    parser.add_argument(
        "-o",
        "--outdir",
        dest="outdirname",
        required=True,
        action="store",
        default=None,
        help="Output directory (required)")
    parser.add_argument(
        "-t",
        "--taxon",
        dest="taxon",
        action="store",
        default=None,
        help="NCBI taxonomy ID")
    parser.add_argument(
        "-v",
        "--verbose",
        dest="verbose",
        action="store_true",
        default=False,
        help="Give verbose output")
    parser.add_argument(
        "-f",
        "--force",
        dest="force",
        action="store_true",
        default=False,
        help="Force file overwriting")
    parser.add_argument(
        "--noclobber",
        dest="noclobber",
        action="store_true",
        default=False,
        help="Don't nuke existing files")
    parser.add_argument(
        "-l",
        "--logfile",
        dest="logfile",
        action="store",
        default=None,
        help="Logfile location")
    parser.add_argument(
        "--format",
        dest="format",
        action="store",
        default="fasta",
        help="Output file format [gbk|fasta]")
    parser.add_argument(
        "--email",
        dest="email",
        required=True,
        action="store",
        default=None,
        help="Email associated with NCBI queries (required)")
    parser.add_argument(
        "--retries",
        dest="retries",
        action="store",
        default=20,
        type=int,
        help="Number of Entrez retry attempts per request.")
    parser.add_argument(
        "--batchsize",
        dest="batchsize",
        action="store",
        default=10000,
        type=int,
        help="Entrez record return batch size")
    parser.add_argument(
        "--timeout",
        dest="timeout",
        action="store",
        default=10,
        type=int,
        help="Timeout for URL connection (s)")
    return parser.parse_args()
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])
865809
def pasa(args):
    
    from jcvi.formats.base import SetFile
    from jcvi.formats.gff import Gff

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

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

    fastafile, gffile = args
    transcodergff = fastafile + ".transdecoder.gff3"
    transcodergenomegff = fastafile + ".transdecoder.genome.gff3"
    if need_update((fastafile, gffile), (transcodergff, transcodergenomegff)):
        cmd = "{0}/scripts/pasa_asmbls_to_training_set.dbi".format(opts.pasa_home)
        cmd += " --pasa_transcripts_fasta {0} --pasa_transcripts_gff3 {1}".\
                format(fastafile, gffile)
        sh(cmd)

    completeids = fastafile.rsplit(".", 1)[0] + ".complete.ids"
    if need_update(transcodergff, completeids):
        cmd = "grep complete {0} | cut -f1 | sort -u".format(transcodergff)
        sh(cmd, outfile=completeids)

    complete = SetFile(completeids)
    seen = set()
    completegff = transcodergenomegff.rsplit(".", 1)[0] + ".complete.gff3"
    fw = open(completegff, "w")
    gff = Gff(transcodergenomegff)
    for g in gff:
        a = g.attributes
        if "Parent" in a:
            id = a["Parent"][0]
        else:
            id = a["ID"][0]
        asmbl_id = id.split("|")[0]
        if asmbl_id not in complete:
            continue
        print(g, file=fw)
        if g.type == "gene":
            seen.add(id)

    fw.close()
    logging.debug("A total of {0} complete models extracted to `{1}`.".\
                    format(len(seen), completegff))
748198
def instantiate(args):
    
    p = OptionParser(instantiate.__doc__)
    p.set_annot_reformat_opts()
    p.add_option("--extended_stride", default=False, action="store_true",
                 help="Toggle extended strides for gene numbering")
    opts, args = p.parse_args(args)

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

    taggedbed, blacklist, gapsbed = args
    r = NameRegister(prefix=opts.prefix, pad0=opts.pad0, uc=opts.uc)
    r.get_blacklist(blacklist)
    r.get_gaps(gapsbed)

    # Run through the bed, identify stretch of NEW ids to instantiate,
    # identify the flanking FRAMEs, interpolate!
    bed = Bed(taggedbed)
    outputbed = taggedbed.rsplit(".", 1)[0] + ".new.bed"
    fw = open(outputbed, "w")

    tagkey = lambda x: x.rsplit("|", 1)[-1]
    for chr, sbed in bed.sub_beds():
        current_chr = chr_number(chr)
        if not current_chr:
            continue

        sbed = list(sbed)

        ranks = []
        for i, s in enumerate(sbed):
            nametag = s.extra[0]
            tag = tagkey(nametag)

            if tag in (NEW, FRAME):
                ranks.append((i, nametag))

        blocks = []
        for tag, names in groupby(ranks, key=lambda x: tagkey(x[-1])):
            names = list(names)
            if tag == NEW:
                blocks.append((tag, [sbed[x[0]] for x in names]))
            else:
                start, end = names[0][-1], names[-1][-1]
                start, end = atg_name(start, retval="rank"), atg_name(end, retval="rank")
                blocks.append((tag, [start, end]))

        id_table = {}  # old to new name conversion
        for i, (tag, info) in enumerate(blocks):
            if tag != NEW:
                continue

            start_id = 0 if i == 0 else blocks[i - 1][1][-1]
            end_id = start_id + 10000 if i == len(blocks) -1 \
                        else blocks[i + 1][1][0]

            r.allocate(info, chr, start_id, end_id, id_table, extended_stride=opts.extended_stride)

        # Output new names
        for i, s in enumerate(sbed):
            nametag = s.extra[0]
            name, tag = nametag.split("|")

            if tag == NEW:
                assert name == '.'
                name = id_table[s.accn]
            elif tag == OVERLAP:
                if name in id_table:
                    name = id_table[name]

            s.extra[0] = "|".join((name, tag))
            print(s, file=fw)

    fw.close()
1272204
static long getGenerationStampFromSeperateChecksumFile(String[] listdir, String blockName) {
    for (int j = 0; j < listdir.length; j++) {
      String path = listdir[j];
      if (!path.startsWith(blockName)) {
        continue;
      }
      String[] vals = StringUtils.split(path, '_');
      if (vals.length != 3) {     // blk, blkid, genstamp.meta
        continue;
      }
      String[] str = StringUtils.split(vals[2], '.');
      if (str.length != 2) {
        continue;
      }
      return Long.parseLong(str[0]);
    }
    DataNode.LOG.warn("Block " + blockName +
                      " does not have a metafile!");
    return Block.GRANDFATHER_GENERATION_STAMP;
  }
147056
def main(argv=sys.argv[1:]):
    
    parser = setup_argparse()
    args = parser.parse_args(argv)
    if args.gendochelp:
        sys.argv[0] = 'gendoc.py'
        genparser = gendoc.setup_parser()
        genparser.print_help()
        sys.exit(0)
    print 'Preparing output directories'
    print '='*80
    for odir in args.output:
        prepare_dir(odir, not args.nodelete)
    print '\nRunning gendoc'
    print '='*80
    for i, idir in enumerate(args.input):
        if i >= len(args.output):
            odir = args.output[-1]
        else:
            odir = args.output[i]
        run_gendoc(idir, odir, args.gendocargs)
885499
def demo(args):
    
    p = OptionParser(demo.__doc__)
    opts, args = p.parse_args(args)

    fig = plt.figure(1, (8, 5))
    root = fig.add_axes([0, 0, 1, 1])

    panel_space = .23
    dup_space = .025
    # Draw a gene and two regulatory elements at these arbitrary locations
    locs = [(.5, .9), # ancestral gene
            (.5, .9 - panel_space + dup_space), # identical copies
            (.5, .9 - panel_space - dup_space),
            (.5, .9 - 2 * panel_space + dup_space), # degenerate copies
            (.5, .9 - 2 * panel_space - dup_space),
            (.2, .9 - 3 * panel_space + dup_space), # sub-functionalization
            (.2, .9 - 3 * panel_space - dup_space),
            (.5, .9 - 3 * panel_space + dup_space), # neo-functionalization
            (.5, .9 - 3 * panel_space - dup_space),
            (.8, .9 - 3 * panel_space + dup_space), # non-functionalization
            (.8, .9 - 3 * panel_space - dup_space),
            ]

    default_regulator = "gm"
    regulators = [default_regulator,
            default_regulator, default_regulator,
            "wm", default_regulator,
            "wm", "gw",
            "wb", default_regulator,
            "ww", default_regulator,
            ]

    width = .24
    for i, (xx, yy) in enumerate(locs):
        regulator = regulators[i]
        x1, x2 = xx - .5 * width, xx + .5 * width
        Glyph(root, x1, x2, yy)
        if i == 9:  # upper copy for non-functionalization
            continue

        # coding region
        x1, x2 = xx - .16 * width, xx + .45 * width
        Glyph(root, x1, x2, yy, fc="k")

        # two regulatory elements
        x1, x2 = xx - .4 * width, xx - .28 * width
        for xx, fc in zip((x1, x2), regulator):
            if fc == 'w':
                continue

            DoubleCircle(root, xx, yy, fc=fc)

        rotation = 30
        tip = .02
        if i == 0:
            ya = yy + tip
            root.text(x1, ya, "Flower", rotation=rotation, va="bottom")
            root.text(x2, ya, "Root", rotation=rotation, va="bottom")
        elif i == 7:
            ya = yy + tip
            root.text(x2, ya, "Leaf", rotation=rotation, va="bottom")

    # Draw arrows between panels (center)
    arrow_dist = .08
    ar_xpos = .5
    for ar_ypos in (.3, .53, .76):
        root.annotate(" ", (ar_xpos, ar_ypos),
                (ar_xpos, ar_ypos + arrow_dist),
                arrowprops=arrowprops)

    ar_ypos = .3
    for ar_xpos in (.2, .8):
        root.annotate(" ", (ar_xpos, ar_ypos),
                (.5, ar_ypos + arrow_dist),
                arrowprops=arrowprops)

    # Duplication, Degeneration
    xx = .6
    ys = (.76, .53)
    processes = ("Duplication", "Degeneration")
    for yy, process in zip(ys, processes):
        root.text(xx, yy + .02, process, fontweight="bold")

    # Label of fates
    xs = (.2, .5, .8)
    fates = ("Subfunctionalization", "Neofunctionalization",
            "Nonfunctionalization")
    yy = .05
    for xx, fate in zip(xs, fates):
        RoundLabel(root, xx, yy, fate)

    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()

    figname = "demo.pdf"
    savefig(figname, dpi=300)
1242986
def _add_genes_to_bed(in_file, gene_file, fai_file, out_file, data, max_distance=10000):
    
    try:
        input_rec = next(iter(pybedtools.BedTool(in_file)))
    except StopIteration:  # empty file
        utils.copy_plus(in_file, out_file)
        return
    # keep everything after standard chrom/start/end, 1-based
    extra_fields = list(range(4, len(input_rec.fields) + 1))
    # keep the new gene annotation
    gene_index = len(input_rec.fields) + 4
    extra_fields.append(gene_index)
    columns = ",".join([str(x) for x in extra_fields])
    max_column = max(extra_fields) + 1
    ops = ",".join(["distinct"] * len(extra_fields))
    # swap over gene name to '.' if beyond maximum distance
    # cut removes the last distance column which can cause issues
    # with bedtools merge: 'ERROR: illegal character '.' found in integer conversion of string'
    distance_filter = (r"""awk -F$'\t' -v OFS='\t' '{if ($NF > %s || $NF < -%s) $%s = "."} {print}'""" %
                       (max_distance, max_distance, gene_index))
    sort_cmd = bedutils.get_sort_cmd(os.path.dirname(out_file))
    cat_cmd = "zcat" if in_file.endswith(".gz") else "cat"
    # Ensure gene transcripts match reference genome
    ready_gene_file = os.path.join(os.path.dirname(out_file), "%s-genomeonly.bed" %
                                   (utils.splitext_plus(os.path.basename(gene_file))[0]))
    ready_gene_file = bedutils.subset_to_genome(gene_file, ready_gene_file, data)
    exports = "export TMPDIR=%s && %s" % (os.path.dirname(out_file), utils.local_path_export())
    bcbio_py = sys.executable
    gsort = config_utils.get_program("gsort", data)
    cmd = ("{exports}{cat_cmd} {in_file} | grep -v ^track | grep -v ^browser | grep -v ^# | "
           "{bcbio_py} -c 'from bcbio.variation import bedutils; bedutils.remove_bad()' | "
           "{gsort} - {fai_file} | "
            "bedtools closest -g {fai_file} "
            "-D ref -t first -a - -b <({gsort} {ready_gene_file} {fai_file}) | "
            "{distance_filter} | cut -f 1-{max_column} | "
            "bedtools merge -i - -c {columns} -o {ops} -delim ',' -d -10 > {out_file}")
    do.run(cmd.format(**locals()), "Annotate BED file with gene info")
364255
def getUI(prog_name, args):
  """"""
  longDescription = "Given a set of BED intervals, compute a profile of " +\
                    "conservation by averaging over all intervals using a " +\
                    "whole genome alignment to a set of relevent species." +\
                    "\n\n" +\
                    "Usage: " + prog_name + " [options] regions.bed " +\
                    "genome-alig species" +\
                    "\n\n" +\
                    "genome-alig can be either a single MAF file, or a " +\
                    "directory of MAF files. In the latter case, the " +\
                    "directory may also optionally contain index files for " +\
                    "the alignment files."
  shortDescription = longDescription

  ui = CLI(prog_name, shortDescription, longDescription)
  # gotta have two args -- MAF dir/file and BED regions.
  # Input by stdin not allowed
  ui.minArgs = 3
  ui.maxArgs = 4
  ui.addOption(Option(short="o", long="output", argName="filename",
                      description="output to given file, else stdout",
                      required=False, type=str))
  ui.addOption(Option(short="w", long="window", argName="size",
                      description="size of window to compute around each " +
                                  "interval; " +
                                  str(DEFAULT_WINDOW_SIZE) +
                                  " to use whole interval. " +
                                  "Default " + str(DEFAULT_WINDOW_SIZE),
                      required=False, type=int))
  ui.addOption(Option(short="e", long="extensions", argName="extension",
                      description="if genome-alig specifies a directory, " +
                                  "treat files with this extension as " +
                                  "alignment files.", required=False,
                                  type=str))
  ui.addOption(Option(short="i", long="index-extensions", argName="extension",
                      description="if genome-alig specifies a directory, " +
                                  "treat files with this extension as " +
                                  "index files for alignments.",
                                  required=False, type=str))
  ui.addOption(Option(short="f", long="fail-no-index",
                      description="fail if an alignment file without an " +
                                  "index is found; otherwise index-less " +
                                  "alignment files are loaded whole (which " +
                                  "might be slow if they're large, and " +
                                  "might require a lot of memory)",
                      default=False, required=False))
  ui.addOption(Option(short="m", long="missing", argName="strategy",
                      description="how to treat missing sequences in " +
                                  "blocks. Options are " +
                                  ", ".join([str(x.name) for x in
                                             MissingSequenceHandler]),
                                  required=False, type=str))
  ui.addOption(Option(short="s", long="species", argName="species",
                      description="consider only these species. Default is " +
                                  "all.", required=False, type=str))
  ui.addOption(Option(short="v", long="verbose",
                      description="output additional messages to stderr " +
                                  "about run (default: " +
                                  str(DEFAULT_VERBOSITY) + ")",
                      default=DEFAULT_VERBOSITY, required=False))
  ui.addOption(Option(short="h", long="help",
                      description="show this help message ", special=True))
  ui.addOption(Option(short="u", long="test",
                      description="run unit tests ", special=True))

  ui.parseCommandLine(args)
  return ui