def sequin(args):
    """
    %prog sequin inputfasta

    Generate a gapped fasta format with known gap sizes embedded. suitable for
    Sequin submission.

    A gapped sequence represents a newer method for describing non-contiguous
    sequences, but only requires a single sequence identifier. A gap is
    represented by a line that starts with >? and is immediately followed by
    either a length (for gaps of known length) or "unk100" for gaps of unknown
    length. For example, ">?200". The next sequence segment continues on the
    next line, with no separate definition line or identifier. The difference
    between a gapped sequence and a segmented sequence is that the gapped
    sequence uses a single identifier and can specify known length gaps.
    Gapped sequences are preferred over segmented sequences. A sample gapped
    sequence file is shown here:

    >m_gagei [organism=Mansonia gagei] Mansonia gagei NADH dehydrogenase ...
    ATGGAGCATACATATCAATATTCATGGATCATACCGTTTGTGCCACTTCCAATTCCTATTTTAATAGGAA
    TTGGACTCCTACTTTTTCCGACGGCAACAAAAAATCTTCGTCGTATGTGGGCTCTTCCCAATATTTTATT
    >?200
    GGTATAATAACAGTATTATTAGGGGCTACTTTAGCTCTTGC
    TCAAAAAGATATTAAGAGGGGTTTAGCCTATTCTACAATGTCCCAACTGGGTTATATGATGTTAGCTCTA
    >?unk100
    TCAATAAAACTATGGGGTAAAGAAGAACAAAAAATAATTAACAGAAATTTTCGTTTATCTCCTTTATTAA
    TATTAACGATGAATAATAATGAGAAGCCATATAGAATTGGTGATAATGTAAAAAAAGGGGCTCTTATTAC
    """
    p = OptionParser(sequin.__doc__)
    p.add_option("--unk", default=100, type="int",
            help="The size for unknown gaps [default: %default]")
    p.add_option("--newid", default=None,
            help="Use this identifier instead [default: %default]")
    p.add_option("--chromosome", default=None,
            help="Add [chromosome= ] to FASTA header [default: %default]")
    p.add_option("--clone", default=None,
            help="Add [clone= ] to FASTA header [default: %default]")
    p.set_mingap(default=100)
    opts, args = p.parse_args(args)

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

    inputfasta, = args
    unk = opts.unk

    outputfasta = inputfasta.rsplit(".", 1)[0] + ".split"
    rec = SeqIO.parse(must_open(inputfasta), "fasta").next()
    seq = ""
    unknowns, knowns = 0, 0
    for gap, gap_group in groupby(rec.seq, lambda x: x.upper() == 'N'):
        subseq = "".join(gap_group)
        if gap:
            gap_length = len(subseq)
            if gap_length == unk:
                subseq = "\n>?unk{0}\n".format(unk)
                unknowns += 1
            elif opts.mingap >= gap_length:
                subseq = "\n>?{0}\n".format(gap_length)
                knowns += 1
        seq += subseq

    fw = must_open(outputfasta, "w")
    id = opts.newid or rec.id
    fastaheader = ">{0}".format(id)
    if opts.chromosome:
        fastaheader += " [chromosome={0}]".format(opts.chromosome)
    if opts.clone:
        fastaheader += " [clone={0}]".format(opts.clone)

    print >> fw, fastaheader
    print >> fw, seq
    fw.close()
    logging.debug("Sequin FASTA written to `{0}` (gaps: {1} unknowns, {2} knowns).".\
            format(outputfasta, unknowns, knowns))

    return outputfasta, unknowns + knowns
dataset/ETHPy150Open tanghaibao/jcvi/formats/fasta.py/sequin/swapped operands of `>=`