def main(argv=None):

    if argv is None:
        argv = sys.argv

    parser = E.OptionParser(
        version="%prog version: $Id: fasta2distances.py 2781 2009-09-10 11:33:14Z andreas $")

    parser.add_option("--filters", dest="filters", type="string",
                      help="Filters to use for filtering sequences [all|codon1|codon2|codon3|d4].")
    parser.add_option("--fields", dest="fields", type="string",
                      help="Fields to output [aligned|nunaligned1|nunaligned2|identical|transitions|transversions|jc69|t92].")

    parser.set_defaults(
        filename_map=None,
        filters="all,codon1,codon2,codon3,d4",
        gap_char="-",
        fields="aligned,unaligned1,unaligned2,identical,transitions,transversions,jc69,t92",
    )

    (options, args) = E.Start(parser, add_pipe_options=True)

    options.filters = options.filters.split(",")
    options.fields = options.fields.split(",")

    iterator = FastaIterator.FastaIterator(options.stdin)

    headers = ["id1", "id2"]
    for f in options.filters:
        headers += list(map(lambda x: "%s_%s" % (f, x), options.fields))

    options.stdout.write("\t".join(headers) + "\n")

    while 1:
        try:
            cur_record = iterator.next()
            if cur_record is None:
                break
            first_record = cur_record
            cur_record = iterator.next()
            if cur_record is None:
                break
            second_record = cur_record

        except StopIteration:
            break

        if len(first_record.sequence) != len(second_record.sequence):
            raise "sequences %s and %s of unequal length" % (
                first_record.title, second_record.title)

        if len(first_record.sequence) % 3 != 0:
            raise "sequence %s not multiple of 3" % first_record.title

        # old: Bio.Alphabet.IUPAC.extended_dna.letters
        alphabet = "ACGT" + options.gap_char

        result = []
        for f in options.filters:

            s1, s2 = FilterAlignedPairForPositions(first_record.sequence,
                                                   second_record.sequence,
                                                   f)

            info = Genomics.CalculatePairIndices(s1, s2, options.gap_char)

            for field in options.fields:

                if field == "aligned":
                    c = "%i" % info.mNAligned
                elif field == "unaligned1":
                    c = "%i" % info.mNUnaligned1
                elif field == "unaligned2":
                    c = "%i" % info.mNUnaligned2
                elif field == "transversions":
                    c = "%i" % info.mNTransversions
                elif field == "transitions":
                    c = "%i" % info.mNTransitions
                elif field == "identical":
                    c = "%i" % info.mNIdentical
                elif field == "jc69":
                    try:
                        c = "%6.4f" % CalculateDistanceJC69(info)[0]
                    except ValueError:
                        c = "nan"
                elif field == "t92":
                    try:
                        c = "%6.4f" % CalculateDistanceT92(info)[0]
                    except ValueError:
                        c = "nan"
                else:
                    raise "Unknown field %s" % field

                result.append(c)

        options.stdout.write("%s\t%s\t%s\n" % (first_record.title,
                                               second_record.title,
                                               "\t".join(result)))

    E.Stop()
script main.
dataset/ETHPy150Open CGATOxford/cgat/scripts/fasta2distances.py/main