%prog stack fastafile

    Create landscape plots that show the amounts of genic sequences, and repetitive
    sequences along the chromosomes.
557515
def stack(args):
    
    p = OptionParser(stack.__doc__)
    p.add_option("--top", default=10, type="int",
                 help="Draw the first N chromosomes [default: %default]")
    p.add_option("--stacks",
                 default="Exons,Introns,DNA_transposons,Retrotransposons",
                 help="Features to plot in stackplot [default: %default]")
    p.add_option("--switch",
                 help="Change chr names based on two-column file [default: %default]")
    add_window_options(p)
    opts, args, iopts = p.set_image_options(args, figsize="8x8")

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

    fastafile, = args
    top = opts.top
    window, shift, subtract, merge = check_window_options(opts)
    switch = opts.switch
    if switch:
        switch = DictFile(opts.switch)

    stacks = opts.stacks.split(",")
    bedfiles = get_beds(stacks)
    binfiles = get_binfiles(bedfiles, fastafile, shift,
                            subtract=subtract, merge=merge)

    sizes = Sizes(fastafile)
    s = list(sizes.iter_sizes())[:top]
    maxl = max(x[1] for x in s)
    margin = .08
    inner = .02   # y distance between tracks

    pf = fastafile.rsplit(".", 1)[0]
    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])

    # Gauge
    ratio = draw_gauge(root, margin, maxl)

    # Per chromosome
    yinterval = (1 - 2 * margin) / (top + 1)
    xx = margin
    yy = 1 - margin
    for chr, clen in s:
        yy -= yinterval
        xlen = clen / ratio
        cc = chr
        if "_" in chr:
            ca, cb = chr.split("_")
            cc = ca[0].upper() + cb

        if switch and cc in switch:
            cc = "\n".join((cc, "({0})".format(switch[cc])))

        root.add_patch(Rectangle((xx, yy), xlen, yinterval - inner, color=gray))
        ax = fig.add_axes([xx, yy, xlen, yinterval - inner])

        nbins = clen / shift
        if clen % shift:
            nbins += 1

        stackplot(ax, binfiles, nbins, palette, chr, window, shift)
        root.text(xx - .04, yy + .5 * (yinterval - inner), cc, ha="center", va="center")

        ax.set_xlim(0, nbins)
        ax.set_ylim(0, 1)
        ax.set_axis_off()

    # Legends
    yy -= yinterval
    xx = margin
    for b, p in zip(bedfiles, palette):
        b = b.rsplit(".", 1)[0].replace("_", " ")
        b = Registration.get(b, b)

        root.add_patch(Rectangle((xx, yy), inner, inner, color=p, lw=0))
        xx += 2 * inner
        root.text(xx, yy, b, size=13)
        xx += len(b) * .012 + inner

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

    image_name = pf + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
1278208
def write_last_and_beds(pf, GenePositions, ContigStarts):
    
    qbedfile = pf + "tigs.bed"
    sbedfile = pf + "chr.bed"
    lastfile = "{}tigs.{}chr.last".format(pf, pf)
    qbedfw = open(qbedfile, "w")
    sbedfw = open(sbedfile, "w")
    lastfw = open(lastfile, "w")

    GeneContigs = np.searchsorted(ContigStarts, GenePositions) - 1
    for i, (c, gstart) in enumerate(zip(GeneContigs, GenePositions)):
        gene = "gene{:05d}".format(i)
        tig = "tig{:04d}".format(c)
        start = ContigStarts[c]
        cstart = gstart - start
        print("\t".join(str(x) for x in
                                   (tig, cstart, cstart + 1, gene)), file=qbedfw)
        print("\t".join(str(x) for x in
                                   ("chr1", gstart, gstart + 1, gene)), file=sbedfw)
        lastatoms = [gene, gene, 100] + [0] * 8 + [100]
        print("\t".join(str(x) for x in lastatoms), file=lastfw)

    qbedfw.close()
    sbedfw.close()
    lastfw.close()
1244672
def roc(args):
    
    outputfile = args.outfile
    # Default extension for image
    if outputfile and not outputfile.endswith(".png"):
        outputfile += ".png"
    
    motifs = read_motifs(args.pwmfile, fmt="pwm")

    ids = []
    if args.ids:
        ids = args.ids.split(",")
    else:
        ids = [m.id for m in motifs]
    motifs = [m for m in motifs if (m.id in ids)]
    
    stats = [
            "phyper_at_fpr",
            "roc_auc", 
            "pr_auc", 
            "enr_at_fpr",
            "recall_at_fdr", 
            "roc_values",
            "matches_at_fpr",
            ]
    
    plot_x = []
    plot_y = []
    legend = []
    
    f_out = sys.stdout
    if args.outdir:
        if not os.path.exists(args.outdir):
            os.makedirs(args.outdir)
        f_out = open(args.outdir + "/gimme.roc.report.txt", "w")
    
    # Print the metrics
    f_out.write("Motif\t# matches\t# matches background\tP-value\tlog10 P-value\tROC AUC\tPR AUC\tEnr. at 1% FPR\tRecall at 10% FDR\n")
    
    
    for motif_stats in calc_stats_iterator(motifs, args.sample, args.background, 
            genome=args.genome, stats=stats, ncpus=args.ncpus):
    
        for motif in motifs:
            if str(motif) in motif_stats:
                if outputfile:
                    x, y = motif_stats[str(motif)]["roc_values"]
                    plot_x.append(x)
                    plot_y.append(y)
                    legend.append(motif.id)
                log_pvalue = np.inf
                if motif_stats[str(motif)]["phyper_at_fpr"] > 0:
                    log_pvalue = -np.log10(motif_stats[str(motif)]["phyper_at_fpr"])
                f_out.write("{}\t{:d}\t{:d}\t{:.2e}\t{:.3f}\t{:.3f}\t{:.3f}\t{:.2f}\t{:0.4f}\n".format(
                      motif.id, 
                      motif_stats[str(motif)]["matches_at_fpr"][0], 
                      motif_stats[str(motif)]["matches_at_fpr"][1], 
                      motif_stats[str(motif)]["phyper_at_fpr"], 
                      log_pvalue, 
                      motif_stats[str(motif)]["roc_auc"], 
                      motif_stats[str(motif)]["pr_auc"], 
                      motif_stats[str(motif)]["enr_at_fpr"], 
                      motif_stats[str(motif)]["recall_at_fdr"],
                      ))
    f_out.close() 
    
    if args.outdir:
        html_report(
            args.outdir,
            args.outdir + "/gimme.roc.report.txt",
            args.pwmfile,
            0.01,
            )

    # Plot the ROC curve
    if outputfile:
        roc_plot(outputfile, plot_x, plot_y, ids=legend)
416522
def write_fasta(
    init_fasta, info_frags, output=DEFAULT_NEW_GENOME_NAME, junction=False
):

    

    init_genome = {
        record.id: record.seq for record in SeqIO.parse(init_fasta, "fasta")
    }
    my_new_records = []
    with open(info_frags, "r") as info_frags_handle:
        current_seq = ""
        current_id = None
        previous_contig = None
        for line in info_frags_handle:
            if line.startswith(">"):
                previous_contig = None
                if current_id is not None:
                    new_record = SeqRecord(
                        current_seq, id=current_id, description=""
                    )
                    my_new_records.append(new_record)
                current_seq = ""
                current_id = str(line[1:])
            elif line.startswith("init_contig"):
                previous_contig = None
            else:
                (init_contig, _, orientation, pos_start, pos_end) = str(
                    line[:-1]
                ).split("\t")

                start = int(pos_start)
                end = int(pos_end)
                ori = int(orientation)

                assert start < end
                assert ori in {-1, 1}

                seq_to_add = init_genome[init_contig][start:end]
                if ori == 1:
                    current_seq += seq_to_add
                elif ori == -1:
                    current_seq += seq_to_add.reverse_complement()

                if junction and previous_contig not in {init_contig, None}:
                    error_was_raised = False
                    try:
                        extra_seq = Seq(junction, IUPAC.ambiguous_dna)
                        current_seq = extra_seq + current_seq
                    except TypeError:
                        if not error_was_raised:
                            print("Invalid junction sequence")
                            error_was_raised = True
                previous_contig = init_contig

        new_record = SeqRecord(current_seq, id=current_id, description="")
        my_new_records.append(new_record)
    SeqIO.write(my_new_records, output, "fasta")
167965
def needle(args):
    
    from jcvi.formats.fasta import Fasta, SeqIO

    p = OptionParser(needle.__doc__)

    opts, args = p.parse_args(args)

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

    manager = mp.Manager()
    results = manager.list()
    needle_pool = mp.Pool(processes=mp.cpu_count())

    pairsfile, apep, bpep = args
    afasta, bfasta = Fasta(apep), Fasta(bpep)
    fp = must_open(pairsfile)
    for i, row in enumerate(fp):
        a, b = row.split()
        a, b = afasta[a], bfasta[b]
        fa, fb = must_open("{0}_{1}_a.fasta".format(pairsfile, i), "w"), \
            must_open("{0}_{1}_b.fasta".format(pairsfile, i), "w")
        SeqIO.write([a], fa, "fasta")
        SeqIO.write([b], fb, "fasta")
        fa.close()
        fb.close()

        needlefile = "{0}_{1}_ab.needle".format(pairsfile, i)
        needle_pool.apply_async(_needle, \
            (fa.name, fb.name, needlefile, a.id, b.id, results))

    needle_pool.close()
    needle_pool.join()

    fp.close()

    scoresfile = "{0}.scores".format(pairsfile.rsplit(".")[0])
    fw = must_open(scoresfile, "w")
    for result in results:
        print(result, file=fw)
    fw.close()
642252
def fasta2bed(fastafile):
    
    dustfasta = fastafile.rsplit(".", 1)[0] + ".dust.fasta"
    for name, seq in parse_fasta(dustfasta):
        for islower, ss in groupby(enumerate(seq), key=lambda x: x[-1].islower()):
            if not islower:
                continue
            ss = list(ss)
            ms, mn = min(ss)
            xs, xn = max(ss)
            print("\t".join(str(x) for x in (name, ms, xs)))
1011581
def ploidy(args):
    
    p = OptionParser(ploidy.__doc__)
    opts, args, iopts = p.set_image_options(args, figsize="8x7")

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

    seqidsfile, klayout = args

    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])

    Karyotype(fig, root, seqidsfile, klayout)

    fc = "darkslategrey"
    radius = .012
    ot = -.05  # use this to adjust vertical position of the left panel
    TextCircle(root, .1, .9 + ot, r'$\gamma$', radius=radius, fc=fc)
    root.text(.1, .88 + ot, r"$\times3$", ha="center", va="top", color=fc)
    TextCircle(root, .08, .79 + ot, r'$\alpha$', radius=radius, fc=fc)
    TextCircle(root, .12, .79 + ot, r'$\beta$', radius=radius, fc=fc)
    root.text(.1, .77 + ot, r"$\times3\times2\times2$", ha="center", va="top", color=fc)
    root.text(.1, .67 + ot, r"Brassica triplication", ha="center",
                va="top", color=fc, size=11)
    root.text(.1, .65 + ot, r"$\times3\times2\times2\times3$", ha="center", va="top", color=fc)
    root.text(.1, .42 + ot, r"Allo-tetraploidy", ha="center",
                va="top", color=fc, size=11)
    root.text(.1, .4 + ot, r"$\times3\times2\times2\times3\times2$", ha="center", va="top", color=fc)

    bb = dict(boxstyle="round,pad=.5", fc="w", ec="0.5", alpha=0.5)
    root.text(.5, .2 + ot, r"\noindent\textit{Brassica napus}\\"
                "(A$\mathsf{_n}$C$\mathsf{_n}$ genome)", ha="center",
                size=16, color="k", bbox=bb)

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

    pf = "napus"
    image_name = pf + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
841006
def get_chromosome_lengths(fasta_file, fancy_sort=True):
    """"""
    chromlen = []
    with gzip.open(fasta_file, 'rt', encoding='ascii') as fh:
        fasta = SeqIO.parse(fh, 'fasta')
        for i, f in enumerate(fasta):
            chromlen.append((f.id, len(f.seq)))
            _LOGGER.info('Processed chromosome "%s"...', f.id)
            #print(dir(f))
            #if i == 1: break
            
    # convert to pandas Series
    chromlen = pd.Series(OrderedDict(chromlen))
    chromlen.index.name = 'Chromosome'
    chromlen.name = 'Length'

    if fancy_sort:
        # sort using fancy ordering
        chrom_for_sorting = chromlen.index.to_series().apply(_transform_chrom)
        a = chrom_for_sorting.argsort(kind='mergesort')
        chromlen = chromlen.iloc[a]

    return chromlen
574294
def add_options(parser):
    
    partial_action = common.partial_append_action
    file_mods = parser.add_argument_group("Sequence File Modification")
    file_mods.add_argument('--line-wrap', dest='line_wrap', metavar='N',
        type=int, help='Adjust line wrap for sequence strings.  '
        'When N is 0, all line breaks are removed. Only fasta files '
        'are supported for the output format.')
    file_mods.add_argument('--sort', dest='sort',
        choices=['length-asc', 'length-desc', 'name-asc', 'name-desc'],
        help='Perform sorting by length or name, ascending or descending. '
        'ASCII sorting is performed for names')

    parser.epilog = """Filters using regular expressions are case-sensitive
    by default. Append "(?i)" to a pattern to make it case-insensitive."""

    seq_mods = parser.add_argument_group("Sequence Modificaton")
    seq_mods.add_argument('--apply-function', type=module_function,
            metavar='/path/to/module.py:function_name[:parameter]',
            help="""Specify a custom function to apply to the input sequences,
            specified as /path/to/file.py:function_name. Function should accept
            an iterable of Bio.SeqRecord objects, and yield SeqRecords. If the
            parameter is specified, it will be passed as a string as the second
            argument to the function. Specify more than one to chain.""",
            default=[], action='append')
    seq_mods.add_argument('--cut', dest='transforms',
            metavar="start:end[,start2:end2]",
            type=common.sequence_slices,
            action=partial_action(transform.multi_cut_sequences, 'slices'),
            help="""Keep only the residues within the 1-indexed start and end
            positions specified, : separated. Includes last item. Start or end
            can be left unspecified to indicate start/end of sequence. A
            negative start may be provided to indicate an offset from the end
            of the sequence. Note that to prevent negative numbers being
            interpreted as flags, this should be written with an equals
            sign between `--cut` and the argument, e.g.: `--cut=-10:`""")
    seq_mods.add_argument('--relative-to', dest='cut_relative', metavar='ID',
            help="""Apply --cut relative to the indexes of non-gap residues in
            sequence identified by ID""")
    seq_mods.add_argument('--drop', dest='transforms',
            metavar='start:end[,start2:end2]',
            type=common.sequence_slices,
            action=partial_action(transform.drop_columns, 'slices'),
            help="""Remove the residues at the specified indices. Same format as `--cut`.""")
    seq_mods.add_argument('--dash-gap',
            action=partial_action(transform.dashes_cleanup), dest='transforms',
            help="""Replace any of the characters "?.:~" with a "-" for all
            sequences""")
    seq_mods.add_argument('--lower',
            action=partial_action(transform.lower_sequences),
            dest='transforms', help='Translate the sequences to lower case')
    seq_mods.add_argument('--mask', metavar="start1:end1[,start2:end2]",
            action=partial_action(transform.multi_mask_sequences, 'slices'),
            type=common.sequence_slices, dest='transforms', help="""Replace
            residues in 1-indexed slice with gap-characters. If --relative-to
            is also specified, coordinates are relative to the sequence ID
            provided.""")
    seq_mods.add_argument('--reverse',
            action=partial_action(transform.reverse_sequences),
            dest='transforms', help='Reverse the order of sites in sequences')
    seq_mods.add_argument('--reverse-complement', dest='transforms',
            action=partial_action(transform.reverse_complement_sequences),
            help='Convert sequences into reverse complements')
    seq_mods.add_argument('--squeeze', action=partial_action(transform.squeeze),
            dest='transforms',
            help='''Remove any gaps that are present in the same
            position across all sequences in an alignment (equivalent to
            --squeeze-threshold=1.0)''')
    seq_mods.add_argument('--squeeze-threshold', dest='transforms',
            action=partial_action(transform.squeeze, 'gap_threshold'),
            type=common.typed_range(float, 0.0, 1.0),
            metavar='PROP', help="""Trim columns from an alignment which
            have gaps in least the specified proportion of sequences.""")
    seq_mods.add_argument('--transcribe', dest='transforms',
            action=partial_action(transform.transcribe, 'transcribe'),
            choices=('dna2rna', 'rna2dna'), help="""Transcription and back
            transcription for generic DNA and RNA. Source sequences must be the
            correct alphabet or this action will likely produce incorrect
            results.""")
    seq_mods.add_argument('--translate', dest='transforms',
            action=partial_action(transform.translate, 'translate'),
            choices=['dna2protein', 'rna2protein', 'dna2proteinstop',
                'rna2proteinstop'], help="""Translate from generic DNA/RNA to
            proteins. Options with "stop" suffix will NOT translate through
            stop codons .  Source sequences must be the correct alphabet or
            this action will likely produce incorrect results.""")
    seq_mods.add_argument('--ungap',
            action=partial_action(transform.ungap_sequences),
            dest='transforms', help='Remove gaps in the sequence alignment')
    seq_mods.add_argument('--upper',
            action=partial_action(transform.upper_sequences),
            dest='transforms', help='Translate the sequences to upper case')

    seq_select = parser.add_argument_group("Record Selection")

    seq_select.add_argument('--deduplicate-sequences',
        action='store_const', const=None, default=False,
         dest='deduplicate_sequences', help='Remove any duplicate sequences '
         'by sequence content, keep the first instance seen')
    seq_select.add_argument('--deduplicated-sequences-file', action='store',
        metavar='FILE', dest='deduplicate_sequences', default=False,
        type=common.FileType('wt'),
        help='Write all of the deduplicated sequences to a file')
    seq_select.add_argument('--deduplicate-taxa',
            action=partial_action(transform.deduplicate_taxa),
            dest='transforms', help="""Remove any duplicate sequences by ID,
            keep the first instance seen""")
    seq_select.add_argument('--exclude-from-file', metavar='FILE',
            type=common.FileType('rt'), help="""Filter sequences, removing
            those sequence IDs in the specified file""", dest='transforms',
            action=partial_action(transform.exclude_from_file, 'handle'))
    seq_select.add_argument('--include-from-file', metavar='FILE',
            type=common.FileType('rt'), help="""Filter sequences, keeping only
            those sequence IDs in the specified file""", dest='transforms',
            action=partial_action(transform.include_from_file, 'handle'))
    seq_select.add_argument('--head', metavar='N', dest='transforms',
            action=partial_action(transform.head, 'head'), help="""Trim
            down to top N sequences. With the leading `-', print all but the last N sequences.""")
    seq_select.add_argument('--max-length', dest='transforms', metavar='N',
            action=partial_action(transform.max_length_discard, 'max_length'),
            type=int, help="""Discard any sequences beyond the specified
            maximum length.  This operation occurs *before* all length-changing
            options such as cut and squeeze.""")
    seq_select.add_argument('--min-length', dest='transforms', metavar='N',
            action=partial_action(transform.min_length_discard, 'min_length'),
            type=int, help="""Discard any sequences less than the specified
            minimum length.  This operation occurs *before* cut and squeeze.""")
    seq_select.add_argument('--min-ungapped-length', metavar='N',
            action=partial_action(transform.min_ungap_length_discard,
                'min_length'), type=int, help="""Discard any sequences less
                than the specified minimum length, excluding gaps. This
                operation occurs *before* cut and squeeze.""",
                dest='transforms')
    seq_select.add_argument('--pattern-include', metavar='REGEX',
            action=partial_action(transform.name_include, 'filter_regex'),
            dest='transforms', help="""Filter the sequences by regular
            expression in ID or description""")
    seq_select.add_argument('--pattern-exclude', metavar='REGEX',
            action=partial_action(transform.name_exclude, 'filter_regex'),
            dest='transforms', help="""Filter the sequences by regular
            expression in ID or description""")
    seq_select.add_argument('--prune-empty',
            action=partial_action(transform.prune_empty), dest='transforms',
            help="Prune sequences containing only gaps ('-')")
    seq_select.add_argument('--sample', metavar='N', dest='transforms', type=int,
            action=partial_action(transform.sample, 'k'),
            help = """ Select a random sampling of sequences """)
    seq_select.add_argument('--sample-seed', metavar='N', type=int,
            help = """Set random seed for sampling of sequences""")
    seq_select.add_argument('--seq-pattern-include', metavar='REGEX',
            action=partial_action(transform.seq_include, 'filter_regex'),
            dest='transforms', help="""Filter the sequences by regular
            expression in sequence""")
    seq_select.add_argument('--seq-pattern-exclude', metavar='REGEX',
            action=partial_action(transform.seq_exclude, 'filter_regex'),
            dest='transforms', help="""Filter the sequences by regular
            expression in sequence""")
    seq_select.add_argument('--tail', metavar='N', dest='transforms',
            action=partial_action(transform.tail, 'tail'),
        help="""Trim down to bottom N sequences.  Use +N to output sequences starting with the Nth.""")

    id_mods = parser.add_argument_group("Sequence ID Modification")
    id_mods.add_argument('--first-name',
            action=partial_action(transform.first_name_capture),
            dest='transforms', help='''Take only the first whitespace-delimited
            word as the name of the sequence''')
    id_mods.add_argument('--name-suffix', metavar='SUFFIX',
            action=partial_action(transform.name_append_suffix, 'suffix'),
            dest='transforms', help='Append a suffix to all IDs.')
    id_mods.add_argument('--name-prefix', metavar='PREFIX',
            action=partial_action(transform.name_insert_prefix, 'prefix'),
            dest='transforms', help="""Insert a prefix for all
            IDs.""")
    id_mods.add_argument('--pattern-replace', nargs=2,
            metavar=('search_pattern', 'replace_pattern'),
            action=partial_action(transform.name_replace, ('search_regex',
                'replace_pattern')),
            dest='transforms', help="""Replace regex pattern "search_pattern"
            with "replace_pattern" in sequence ID and description""")
    id_mods.add_argument('--strip-range', dest='transforms',
            action=partial_action(transform.strip_range), help="""Strip ranges
            from sequences IDs, matching </x-y>""")

    format_group = parser.add_argument_group('Format Options')
    format_group.add_argument('--input-format', metavar='FORMAT',
            help="Input file format (default: determine from extension)")
    format_group.add_argument('--output-format', metavar='FORMAT',
            help="Output file format (default: determine from extension)")

    parser.add_argument('--alphabet', choices=ALPHABETS,
            help="""Input alphabet. Required for writing NEXUS.""")

    return parser
423875
def fastq(args):
    
    from jcvi.formats.fastq import guessoffset

    p = OptionParser(fastq.__doc__)
    p.add_option("--outtie", dest="outtie", default=False, action="store_true",
            help="Are these outie reads? [default: %default]")
    p.set_phred()
    p.set_size()

    opts, args = p.parse_args(args)

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

    fastqfiles = [get_abs_path(x) for x in args]
    size = opts.size
    outtie = opts.outtie
    if size > 1000 and (not outtie):
        logging.debug("[warn] long insert size {0} but not outtie".format(size))

    mated = (size != 0)
    libname = op.basename(args[0]).split(".")[0]
    libname = libname.replace("_1_sequence", "")

    frgfile = libname + ".frg"
    mean, sv = get_mean_sv(opts.size)

    cmd = "fastqToCA"
    cmd += " -libraryname {0} ".format(libname)
    fastqs = " ".join("-reads {0}".format(x) for x in fastqfiles)
    if mated:
        assert len(args) in (1, 2), "you need one or two fastq files for mated library"
        fastqs = "-mates {0}".format(",".join(fastqfiles))
        cmd += "-insertsize {0} {1} ".format(mean, sv)
    cmd += fastqs

    offset = int(opts.phred) if opts.phred else guessoffset([fastqfiles[0]])
    illumina = (offset == 64)
    if illumina:
        cmd += " -type illumina"
    if outtie:
        cmd += " -outtie"

    sh(cmd, outfile=frgfile)
296192
def clean(args):
    
    p = OptionParser(clean.__doc__)
    p.add_option("--fancy", default=False, action="store_true",
                 help="Pretty print the sequence [default: %default]")
    p.add_option("--canonical", default=False, action="store_true",
                 help="Use only acgtnACGTN [default: %default]")
    p.set_outfile()

    opts, args = p.parse_args(args)

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

    fastafile, = args
    fw = must_open(opts.outfile, "w")
    if opts.fancy:
        for header, seq in iter_clean_fasta(fastafile):
            print(">" + header, file=fw)
            fancyprint(fw, seq)

        return 0

    iterator = iter_canonical_fasta if opts.canonical else iter_clean_fasta

    for header, seq in iterator(fastafile):
        seq = Seq(seq)
        s = SeqRecord(seq, id=header, description="")
        SeqIO.write([s], fw, "fasta")
1273575
def write_sequences_to_fasta(path, seqs):
    
    from Bio import SeqIO
    from Bio.Seq import Seq
    from Bio.SeqRecord import SeqRecord

    path = Path(path)
    records = []

    for id, seq in seqs.items():
        record = SeqRecord(Seq(seq), id=id, description='')
        records.append(record)

    SeqIO.write(records, str(path), 'fasta')
805578
def from23andme(args):
    
    p = OptionParser(from23andme.__doc__)
    p.set_ref()
    opts, args = p.parse_args(args)

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

    txtfile, seqid = args
    ref_dir = opts.ref
    fastafile = op.join(ref_dir, "hs37d5.fa")
    fasta = Fasta(fastafile)

    pf = txtfile.rsplit(".", 1)[0]
    px = CM[seqid]
    chrvcf = pf + ".{0}.vcf".format(px)
    legend = op.join(ref_dir, "1000GP_Phase3/{0}.rsids".format(px))
    register = read_rsid(seqid, legend)

    fw = open(chrvcf, "w")
    print(get_vcfstanza(fastafile, fasta, txtfile), file=fw)

    fp = open(txtfile)
    seen = set()
    duplicates = skipped = missing = 0
    for row in fp:
        if row[0] == '#':
            continue
        rsid, chr, pos, genotype = row.split()
        if chr != seqid:
            continue
        pos = int(pos)
        if (chr, pos) in seen:
            duplicates += 1
            continue
        seen.add((chr, pos))
        genotype = list(genotype)
        if "-" in genotype:  # missing daa
            missing += 1
            continue

        # Y or MT
        if not register:
            assert len(genotype) == 1
            ref = fasta[chr][pos - 1].seq.upper()
            if "D" in genotype or "I" in genotype:
                skipped += 1
                continue
            genotype = genotype[0]
            code = "0/0" if ref == genotype else "1/1"
            alt = "." if ref == genotype else genotype
            print("\t".join(str(x) for x in \
                    (chr, pos, rsid, ref, alt, ".", ".", "PR", "GT", code)), file=fw)
            continue

        # If rsid is seen in the db, use that
        if rsid in register:
            pos, ref, alt = register[rsid]
        elif pos in register:
            pos, ref, alt = register[pos]
        else:
            skipped += 1  # Not in reference panel
            continue

        assert fasta[chr][pos - 1:pos + len(ref) - 1].seq.upper() == ref
        # Keep it bi-allelic
        not_seen = [x for x in alt if x not in genotype]
        while len(alt) > 1 and not_seen:
            alt.remove(not_seen.pop())
        if len(alt) > 1:
            alt = [alt[0]]
        alleles = [ref] + alt

        if len(genotype) == 1:
            genotype = [genotype[0]] * 2

        alt = ",".join(alt) or "."
        if "D" in genotype or "I" in genotype:
            max_allele = max((len(x), x) for x in alleles)[1]
            alleles = [("I" if x == max_allele else "D") for x in alleles]
            assert "I" in alleles and "D" in alleles
        a, b = genotype
        try:
            ia, ib = alleles.index(a), alleles.index(b)
        except ValueError:  # alleles not seen
            logging.error("{0}: alleles={1}, genotype={2}".\
                            format(rsid, alleles, genotype))
            skipped += 1
            continue
        code = "/".join(str(x) for x in sorted((ia, ib)))

        print("\t".join(str(x) for x in \
                (chr, pos, rsid, ref, alt, ".", ".", "PR", "GT", code)), file=fw)

    logging.debug("duplicates={0} skipped={1} missing={2}".\
                    format(duplicates, skipped, missing))
1134260
static public LinkedHashMap<String, ChromosomeSequence> loadFastaAddGeneFeaturesFromGeneIDGFF2(File fastaSequenceFile, File gffFile) throws Exception {
		LinkedHashMap<String, DNASequence> dnaSequenceList = FastaReaderHelper.readFastaDNASequence(fastaSequenceFile);
		LinkedHashMap<String, ChromosomeSequence> chromosomeSequenceList = GeneFeatureHelper.getChromosomeSequenceFromDNASequence(dnaSequenceList);
		FeatureList listGenes = GeneIDGFF2Reader.read(gffFile.getAbsolutePath());
		addGeneIDGFF2GeneFeatures(chromosomeSequenceList, listGenes);
		return chromosomeSequenceList;
	}
511061
def fasta(args):
    
    from jcvi.formats.sizes import Sizes
    from jcvi.formats.agp import OO, build

    p = OptionParser(fasta.__doc__)
    opts, args = p.parse_args(args)

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

    bedfile, scffasta, pmolfasta = args
    pf = bedfile.rsplit(".", 1)[0]
    bed = Bed(bedfile)
    selected = select_bed(bed)
    oo = OO()
    seen = set()
    sizes = Sizes(scffasta).mapping
    agpfile = pf + ".agp"
    agp = open(agpfile, "w")
    for b in selected:
        scf = range_parse(b.accn).seqid
        chr = b.seqid
        cs = (chr, scf)
        if cs not in seen:
            oo.add(chr, scf, sizes[scf], b.strand)
            seen.add(cs)
        else:
            logging.debug("Seen {0}, ignored.".format(cs))

    oo.write_AGP(agp, gaptype="contig")
    agp.close()
    build([agpfile, scffasta, pmolfasta])
1082246
def make_matepairs(fastafile):
    
    assert op.exists(fastafile)

    matefile = fastafile.rsplit(".", 1)[0] + ".mates"
    if op.exists(matefile):
        logging.debug("matepairs file `{0}` found".format(matefile))
    else:
        logging.debug("parsing matepairs from `{0}`".format(fastafile))
        matefw = open(matefile, "w")
        it = SeqIO.parse(fastafile, "fasta")
        for fwd, rev in zip(it, it):
            print("{0}\t{1}".format(fwd.id, rev.id), file=matefw)

        matefw.close()

    return matefile
657568
def insert_sequences_into_tree(seqs, moltype, params={},
                                           write_log=True):
    

    ih = '_input_as_multiline_string'

    raxml_app = Raxml(params=params,
                      InputHandler=ih,
                      WorkingDir=None,
                      SuppressStderr=False,
                      SuppressStdout=False,
                      HALT_EXEC=False)

    raxml_result = raxml_app(seqs)

    # write a log file
    if write_log:
        log_fp = join(params["-w"],'log_raxml_'+split(get_tmp_filename())[-1])
        log_file=open(log_fp,'w')
        log_file.write(raxml_result['StdOut'].read())
        log_file.close()

    '''
    # getting setup since parsimony doesn't output tree..only jplace, however
    # it is currently corrupt

    # use guppy to convert json file into a placement tree
    guppy_params={'tog':None}

    new_tree=build_tree_from_json_using_params(raxml_result['json'].name, \
                                               output_dir=params["-w"], \
                                               params=guppy_params)
    '''

    # get tree from 'Result Names'
    new_tree=raxml_result['Result'].readlines()
    filtered_tree=re.sub('\[I\d+\]','',str(new_tree))
    tree = DndParser(filtered_tree, constructor=PhyloNode)

    raxml_result.cleanUp()

    return tree