module TestSeq

using FactCheck,
    Bio.Seq,
    YAML

import ..get_bio_fmt_specimens

# Return a random DNA/RNA sequence of the given length
function random_seq(n::Integer, nts, probs)
    cumprobs = cumsum(probs)
    x = Array(Char, n)
    for i in 1:n
        x[i] = nts[searchsorted(cumprobs, rand()).start]
    end
    return convert(AbstractString, x)
end


function random_array(n::Integer, elements, probs)
    cumprobs = cumsum(probs)
    x = Array(eltype(elements), n)
    for i in 1:n
        x[i] = elements[searchsorted(cumprobs, rand()).start]
    end
    return x
end


function random_dna(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'T', 'N'], probs)
end


function random_rna(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'U', 'N'], probs)
end

const codons = [
        "AAA", "AAC", "AAG", "AAU",
        "ACA", "ACC", "ACG", "ACU",
        "AGA", "AGC", "AGG", "AGU",
        "AUA", "AUC", "AUG", "AUU",
        "CAA", "CAC", "CAG", "CAU",
        "CCA", "CCC", "CCG", "CCU",
        "CGA", "CGC", "CGG", "CGU",
        "CUA", "CUC", "CUG", "CUU",
        "GAA", "GAC", "GAG", "GAU",
        "GCA", "GCC", "GCG", "GCU",
        "GGA", "GGC", "GGG", "GGU",
        "GUA", "GUC", "GUG", "GUU",
        "UAC", "UAU", "UCA", "UCC",
        "UCG", "UCU", "UGC", "UGG",
        "UGU", "UUA", "UUC", "UUG",
        "UUU",
        # translatable ambiguities in the standard code
        "CUN", "CCN", "CGN", "ACN",
        "GUN", "GCN", "GGN", "UCN"
        ]


function random_translatable_rna(n)
    probs = fill(1.0 / length(codons), length(codons))
    cumprobs = cumsum(probs)
    r = rand()
    x = Array(AbstractString, n)
    for i in 1:n
        x[i] = codons[searchsorted(cumprobs, rand()).start]
    end

    return string(x...)
end


function random_dna_kmer(len)
    return random_dna(len, [0.25, 0.25, 0.25, 0.25])
end


function random_rna_kmer(len)
    return random_rna(len, [0.25, 0.25, 0.25, 0.25])
end


function random_dna_kmer_nucleotides(len)
    return random_array(len, [DNA_A, DNA_C, DNA_G, DNA_T],
                        [0.25, 0.25, 0.25, 0.25])
end


function random_rna_kmer_nucleotides(len)
    return random_array(len, [RNA_A, RNA_C, RNA_G, RNA_U],
                        [0.25, 0.25, 0.25, 0.25])
end


function random_aa(len)
    return random_seq(len,
        ['A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I',
         'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V', 'X' ],
        push!(fill(0.049, 20), 0.02))
end


function random_interval(minstart, maxstop)
    start = rand(minstart:maxstop)
    return start:rand(start:maxstop)
end

facts("Nucleotides") do
    context("Conversions") do
        context("UInt8") do
            context("DNA conversions from UInt8") do
                @fact convert(DNANucleotide, UInt8(0)) --> DNA_A
                @fact convert(DNANucleotide, UInt8(1)) --> DNA_C
                @fact convert(DNANucleotide, UInt8(2)) --> DNA_G
                @fact convert(DNANucleotide, UInt8(3)) --> DNA_T
                @fact convert(DNANucleotide, UInt8(4)) --> DNA_N
            end

            context("RNA conversions from UInt8") do
                @fact convert(RNANucleotide, UInt8(0)) --> RNA_A
                @fact convert(RNANucleotide, UInt8(1)) --> RNA_C
                @fact convert(RNANucleotide, UInt8(2)) --> RNA_G
                @fact convert(RNANucleotide, UInt8(3)) --> RNA_U
                @fact convert(RNANucleotide, UInt8(4)) --> RNA_N
            end

            context("DNA conversions to UInt8") do
                @fact convert(UInt8, DNA_A) --> UInt8(0)
                @fact convert(UInt8, DNA_C) --> UInt8(1)
                @fact convert(UInt8, DNA_G) --> UInt8(2)
                @fact convert(UInt8, DNA_T) --> UInt8(3)
                @fact convert(UInt8, DNA_N) --> UInt8(4)
            end

            context("RNA conversions to UInt8") do
                @fact convert(UInt8, RNA_A) --> UInt8(0)
                @fact convert(UInt8, RNA_C) --> UInt8(1)
                @fact convert(UInt8, RNA_G) --> UInt8(2)
                @fact convert(UInt8, RNA_U) --> UInt8(3)
                @fact convert(UInt8, RNA_N) --> UInt8(4)
            end
        end

        context("UInt64") do
            context("DNA conversions from UInt64") do
                @fact convert(DNANucleotide, UInt64(0)) --> DNA_A
                @fact convert(DNANucleotide, UInt64(1)) --> DNA_C
                @fact convert(DNANucleotide, UInt64(2)) --> DNA_G
                @fact convert(DNANucleotide, UInt64(3)) --> DNA_T
                @fact convert(DNANucleotide, UInt64(4)) --> DNA_N
            end

            context("RNA conversions from UInt64") do
                @fact convert(RNANucleotide, UInt64(0)) --> RNA_A
                @fact convert(RNANucleotide, UInt64(1)) --> RNA_C
                @fact convert(RNANucleotide, UInt64(2)) --> RNA_G
                @fact convert(RNANucleotide, UInt64(3)) --> RNA_U
                @fact convert(RNANucleotide, UInt64(4)) --> RNA_N
            end

            context("DNA conversions to UInt64") do
                @fact convert(UInt64, DNA_A) --> UInt64(0)
                @fact convert(UInt64, DNA_C) --> UInt64(1)
                @fact convert(UInt64, DNA_G) --> UInt64(2)
                @fact convert(UInt64, DNA_T) --> UInt64(3)
                @fact convert(UInt64, DNA_N) --> UInt64(4)
            end

            context("RNA conversions to UInt64") do
                @fact convert(UInt64, RNA_A) --> UInt64(0)
                @fact convert(UInt64, RNA_C) --> UInt64(1)
                @fact convert(UInt64, RNA_G) --> UInt64(2)
                @fact convert(UInt64, RNA_U) --> UInt64(3)
                @fact convert(UInt64, RNA_N) --> UInt64(4)
            end
        end

        context("Char") do
            context("DNA conversions from Char") do
                @fact convert(DNANucleotide, 'A') --> DNA_A
                @fact convert(DNANucleotide, 'C') --> DNA_C
                @fact convert(DNANucleotide, 'G') --> DNA_G
                @fact convert(DNANucleotide, 'T') --> DNA_T
                @fact convert(DNANucleotide, 'N') --> DNA_N
            end

            context("RNA conversions from Char") do
                @fact convert(RNANucleotide, 'A') --> RNA_A
                @fact convert(RNANucleotide, 'C') --> RNA_C
                @fact convert(RNANucleotide, 'G') --> RNA_G
                @fact convert(RNANucleotide, 'U') --> RNA_U
                @fact convert(RNANucleotide, 'N') --> RNA_N
            end

            context("DNA conversions to Char") do
                @fact convert(Char, DNA_A) --> 'A'
                @fact convert(Char, DNA_C) --> 'C'
                @fact convert(Char, DNA_G) --> 'G'
                @fact convert(Char, DNA_T) --> 'T'
                @fact convert(Char, DNA_N) --> 'N'
            end

            context("RNA conversions to Char") do
                @fact convert(Char, RNA_A) --> 'A'
                @fact convert(Char, RNA_C) --> 'C'
                @fact convert(Char, RNA_G) --> 'G'
                @fact convert(Char, RNA_U) --> 'U'
                @fact convert(Char, RNA_N) --> 'N'
            end
        end
    end

    context("Show DNA") do
        buf = IOBuffer()
        for nt in [DNA_A, DNA_C, DNA_G, DNA_T, DNA_N]
            show(buf, nt)
        end
        @fact takebuf_string(buf) --> "ACGTN"
    end
    
    context("Show RNA") do
        buf = IOBuffer()
        for nt in [RNA_A, RNA_C, RNA_G, RNA_U, RNA_N]
            show(buf, nt)
        end
        @fact takebuf_string(buf) --> "ACGUN"
    end

    context("Sequences") do
        function dna_complement(seq::AbstractString)
            seqc = Array(Char, length(seq))
            for (i, c) in enumerate(seq)
                if c     ==   'A'
                    seqc[i] = 'T'
                elseif c ==   'C'
                    seqc[i] = 'G'
                elseif c ==   'G'
                    seqc[i] = 'C'
                elseif c ==   'T'
                    seqc[i] = 'A'
                else
                    seqc[i] = 'N'
                end
            end

            return convert(AbstractString, seqc)
        end

        function rna_complement(seq::AbstractString)
            seqc = Array(Char, length(seq))
            for (i, c) in enumerate(seq)
                if c == 'A'
                    seqc[i] = 'U'
                elseif c == 'C'
                    seqc[i] = 'G'
                elseif c == 'G'
                    seqc[i] = 'C'
                elseif c == 'U'
                    seqc[i] = 'A'
                else
                    seqc[i] = 'N'
                end
            end

            return convert(AbstractString, seqc)
        end

        function check_reversal(T, seq)
            return reverse(seq) == convert(AbstractString, reverse(convert(T, seq)))
        end

        function check_dna_complement(T, seq)
            return dna_complement(seq) ==
                convert(AbstractString, complement(convert(T, seq)))
        end

        function check_rna_complement(T, seq)
            return rna_complement(seq) ==
                convert(AbstractString, complement(convert(T, seq)))
        end

        function check_dna_revcomp(T, seq)
            return reverse(dna_complement(seq)) ==
                convert(AbstractString, reverse_complement(convert(T, seq)))
        end

        function check_rna_revcomp(T, seq)
            return reverse(rna_complement(seq)) ==
                convert(AbstractString, reverse_complement(convert(T, seq)))
        end

        function check_mismatches(T, a, b)
            count = 0
            for (ca, cb) in zip(a, b)
                if ca != cb
                    count += 1
                end
            end
            return mismatches(convert(T, a), convert(T, b)) == count
        end

        context("Nucleotide Sequences") do
            reps = 10
            context("Construction and Conversions") do
                context("Constructing empty sequences") do
                    # Check construction of empty nucleotide sequences
                    #  using RNASequence and DNASequence functions
                    @fact RNASequence() --> NucleotideSequence(RNANucleotide)
                    @fact DNASequence() --> NucleotideSequence(DNANucleotide)
                end

                context("Conversion from/to Strings") do
                    # Check that sequences in strings survive round trip conversion:
                    #   String → NucleotideSequence → String
                    function check_string_construction(T::Type, seq::AbstractString)
                        return convert(AbstractString, NucleotideSequence{T}(seq)) == uppercase(seq)
                    end

                    for len in [0, 1, 10, 32, 1000, 10000, 100000]
                        @fact all(Bool[check_string_construction(DNANucleotide, random_dna(len)) for _ in 1:reps]) --> true
                        @fact all(Bool[check_string_construction(RNANucleotide, random_rna(len)) for _ in 1:reps]) --> true
                        @fact all(Bool[check_string_construction(DNANucleotide, lowercase(random_dna(len))) for _ in 1:reps]) --> true
                        @fact all(Bool[check_string_construction(RNANucleotide, lowercase(random_rna(len))) for _ in 1:reps]) --> true
                    end

                    # Non-nucleotide characters should throw
                    @fact_throws DNASequence("ACCNNCATTTTTTAGATXATAG")
                    @fact_throws RNASequence("ACCNNCATTTTTTAGATXATAG")
                end

                context("Conversion between RNA and DNA") do
                    @fact convert(RNASequence, DNASequence("ACGTN")) --> rna"ACGUN"
                    @fact convert(DNASequence, RNASequence("ACGUN")) --> dna"ACGTN"
                end

                context("Construction from nucleotide vectors") do
                    function check_vector_construction(T::Type, seq::AbstractString)
                        xs = T[convert(T, c) for c in seq]
                        return NucleotideSequence{T}(xs) == NucleotideSequence{T}(seq)
                    end

                    for len in [0, 1, 10, 32, 1000, 10000, 100000]
                        @fact all(Bool[check_vector_construction(DNANucleotide, random_dna(len)) for _ in 1:reps]) --> true
                        @fact all(Bool[check_vector_construction(RNANucleotide, random_rna(len)) for _ in 1:reps]) --> true
                    end
                end

                context("Concatenation") do
                    function check_concatenation(::Type{DNANucleotide}, n)
                        chunks = [random_dna(rand(100:300)) for i in 1:n]
                        parts = Any[]
                        for i in 1:n
                            start = rand(1:length(chunks[i]))
                            stop = rand(start:length(chunks[i]))
                            push!(parts, start:stop)
                        end

                        str = string([chunk[parts[i]]
                                      for (i, chunk) in enumerate(chunks)]...)

                        seq = *([DNASequence(chunk)[parts[i]]
                                 for (i, chunk) in enumerate(chunks)]...)

                        return convert(AbstractString, seq) == uppercase(str)
                    end

                    @fact all(Bool[check_concatenation(DNANucleotide, rand(1:10)) for _ in 1:100]) --> true
                end

                context("Repetition") do
                    function check_repetition(::Type{DNANucleotide}, n)
                        chunk = random_dna(rand(100:300))
                        start = rand(1:length(chunk))
                        stop = rand(start:length(chunk))

                        str = chunk[start:stop] ^ n
                        seq = DNASequence(chunk)[start:stop] ^ n
                        return convert(AbstractString, seq) == uppercase(str)
                    end

                    @fact all(Bool[check_repetition(DNANucleotide, rand(1:10)) for _ in 1:100]) --> true
                end
            end

            context("Equality") do
                reps = 10
                function check_seq_equality(len)
                    a = random_dna(len)
                    return a == copy(a)
                end

                for len in [1, 10, 32, 1000]
                    @fact all(Bool[check_seq_equality(len) for _ in 1:reps]) --> true
                end

                a = b = dna"ACTGN"
                @fact ==(a, b) --> true
                @fact ==(dna"ACTGN", dna"ACTGN") --> true
                @fact ==(dna"ACTGN", dna"ACTGA") --> false
                @fact ==(dna"ACTGN", dna"ACTG") --> false
                @fact ==(dna"ACTG", dna"ACTGN") --> false

                c = d = rna"ACUGN"
                @fact ==(c, d) --> true
                @fact ==(rna"ACUGN", rna"ACUGN") --> true
                @fact ==(rna"ACUGN", rna"ACUGA") --> false
                @fact ==(rna"ACUGN", rna"ACUG") --> false
                @fact ==(rna"ACUG", rna"ACUGN") --> false

                a = dna"ACGTNACGTN"
                b = dna"""
                ACGTN
                ACGTN
                """
                @fact ==(a, b) --> true

                c = rna"ACUGNACUGN"
                d = rna"""
                ACUGN
                ACUGN
                """
                @fact ==(c, d) --> true
            end

            context("Length") do
                for len in [0, 1, 10, 32, 1000, 10000, 100000]
                    @fact length(DNASequence(random_dna(len))) --> len
                    @fact endof(DNASequence(random_dna(len))) --> len
                    @fact length(DNASequence(random_dna(len))) --> len
                    @fact endof(DNASequence(random_dna(len))) --> len
                end
            end

            context("Copy") do
                function check_copy(T, seq)
                    return convert(AbstractString, copy(NucleotideSequence{T}(seq))) == seq
                end

                for len in [1, 10, 32, 1000, 10000, 100000]
                    @fact all(Bool[check_copy(DNANucleotide, random_dna(len)) for _ in 1:reps]) --> true
                    @fact all(Bool[check_copy(RNANucleotide, random_rna(len)) for _ in 1:reps]) --> true
                end
            end

            context("Access and Iterations") do
                dna_seq = dna"ACTG"
                rna_seq = rna"ACUG"

                context("Access DNA Sequence") do
                    # Access indexes out of bounds
                    @fact_throws dna_seq[-1]
                    @fact_throws dna_seq[0]
                    @fact_throws dna_seq[5]
                    @fact_throws getindex(dna_seq,-1)
                    @fact_throws getindex(dna_seq, 0)
                    @fact_throws getindex(dna_seq, 5)
                end

                context("Iteration through DNA Sequence") do
                    @fact start(dna"ACNTG") --> (1,3)
                    @fact start(dna"")      --> (1,1)

                    @fact next(dna"ACTGN", (2,5)) --> (DNA_C, (3,5))
                    @fact next(dna"ACTGN", (5,5)) --> (DNA_N, (6,6))

                    @fact done(dna"", (1,1))       --> true
                    @fact done(dna"ACTGN", (2,5))  --> false
                    @fact done(dna"ACTGN", (5,5))  --> false
                    @fact done(dna"ACTGN", (6,5))  --> true
                    @fact done(dna"ACTGN", (0,5))  --> false

                    dna_vector = [DNA_A, DNA_C, DNA_T, DNA_G]
                    @fact all(Bool[nucleotide == dna_vector[i] for (i, nucleotide) in enumerate(dna_seq)]) -->  true
                end

                context("Access RNA Sequence") do
                    # Access indexes out of bounds
                    @fact_throws rna_seq[-1]
                    @fact_throws rna_seq[0]
                    @fact_throws rna_seq[5]
                    @fact_throws getindex(rna_seq, -1)
                    @fact_throws getindex(rna_seq, 0)
                    @fact_throws getindex(rna_seq, 5)
                end

                context("Iteration through RNA Sequence") do
                    @fact start(rna"ACNUG") --> (1,3)
                    @fact start(rna"")      --> (1,1)

                    @fact next(rna"ACUGN", (2,5)) --> (RNA_C, (3,5))
                    @fact next(rna"ACUGN", (5,5)) --> (RNA_N, (6,6))

                    @fact done(rna"", (1,1))       --> true
                    @fact done(rna"ACUGN", (2,5))  --> false
                    @fact done(rna"ACUGN", (5,5))  --> false
                    @fact done(rna"ACUGN", (6,5))  --> true
                    @fact done(rna"ACUGN", (0,5))  --> false

                    # Iteration through RNA Sequence

                    rna_vector = [RNA_A, RNA_C, RNA_U, RNA_G]
                    @fact all(Bool[nucleotide == rna_vector[i] for (i, nucleotide) in enumerate(rna_seq)]) -->  true
                end

                context("Indexing with Ranges") do
                    @fact getindex(dna"ACTGNACTGN", 1:5) --> dna"ACTGN"
                    @fact getindex(rna"ACUGNACUGN", 1:5) --> rna"ACUGN"
                    @fact getindex(dna"ACTGNACTGN", 5:1) --> dna""
                    @fact getindex(rna"ACUGNACUGN", 5:1) --> rna""
                end
            end

            context("Subsequence Construction") do
                for len in [1, 10, 32, 1000, 10000, 100000]
                    seq = random_dna(len)
                    dnaseq = DNASequence(seq)

                    results = Bool[]
                    for _ in 1:reps
                        part = random_interval(1, length(seq))
                        push!(results, seq[part] == convert(AbstractString, dnaseq[part]))
                    end
                    @fact all(results) --> true
                end

                for len in [1, 10, 32, 1000, 10000, 100000]
                    seq = random_rna(len)
                    rnaseq = RNASequence(seq)

                    results = Bool[]
                    for _ in 1:reps
                        part = random_interval(1, length(seq))

                        push!(results, seq[part] == convert(AbstractString, rnaseq[part]))
                    end
                    @fact all(results) --> true
                end

                context("Subsequence Construction from Ranges") do
                    # Subsequence from range
                    @fact (RNASequence(rna"AUCGAUCG", 5:8) == RNASequence("AUCG")) --> true
                    @fact (DNASequence(dna"ATCGATCG", 5:8) == DNASequence("ATCG")) --> true

                    # Invalid ranges
                    @fact_throws RNASequence(rna"AUCGAUCG", 5:10)
                    @fact_throws DNASequence(dna"ATCGATCG", 5:10)

                    # Empty ranges
                    @fact (RNASequence(rna"AUCGAUCG", 5:4) == RNASequence()) --> true
                    @fact (DNASequence(dna"ATCGATCG", 5:4) == DNASequence()) --> true

                    # Subsequence of subsequence
                    @fact dna"ACGTAG"[4:end][1:2] == dna"TA" --> true
                    @fact dna"ACGTAG"[4:end][2:3] == dna"AG" --> true
                    @fact_throws dna"ACGTAG"[4:end][0:1]
                    @fact_throws dna"ACGTAG"[4:end][3:4]
                end
            end

            context("Transformations") do
                context("Reversal") do
                    for len in [0, 1, 10, 32, 1000, 10000, 100000]
                        @fact all(Bool[check_reversal(DNASequence, random_dna(len)) for _ in 1:reps]) --> true
                        @fact all(Bool[check_reversal(RNASequence, random_rna(len)) for _ in 1:reps]) --> true
                    end
                end

                context("Complement") do
                    for len in [1, 10, 32, 1000, 10000, 100000]
                        @fact all(Bool[check_dna_complement(DNASequence, random_dna(len)) for _ in 1:reps]) --> true
                        @fact all(Bool[check_rna_complement(RNASequence, random_rna(len)) for _ in 1:reps]) --> true
                    end
                end

                context("Reverse Complement") do
                    for len in [1, 10, 32, 1000, 10000, 100000]
                        @fact all(Bool[check_dna_revcomp(DNASequence, random_dna(len)) for _ in 1:reps]) --> true
                        @fact all(Bool[check_rna_revcomp(RNASequence, random_rna(len)) for _ in 1:reps]) --> true
                    end
                end
            end

            context("Mismatches") do
                for len in [1, 10, 32, 1000, 10000, 100000]
                    @fact all(Bool[check_mismatches(DNASequence, random_dna(len), random_dna(len))
                                   for _ in 1:reps]) --> true
                    @fact all(Bool[check_mismatches(RNASequence, random_rna(len), random_rna(len))
                                   for _ in 1:reps]) --> true
                end
            end

            context("Mutability") do
                seq = dna"ACGTACGT"
                @fact ismutable(seq) --> false
                @fact_throws seq[1] = DNA_C

                seq2 = seq[1:4]
                @fact ismutable(seq2) --> false

                mutable!(seq)
                @fact ismutable(seq) --> true

                seq[1] = DNA_C
                @fact seq[1] --> DNA_C
                @fact seq2[1] --> DNA_A

                seq[2] = DNA_N
                @fact seq[2] --> DNA_N
                @fact seq2[2] --> DNA_C

                seq[2] = DNA_G
                @fact seq[2] --> DNA_G

                immutable!(seq)
                @fact_throws seq[1] = DNA_A

                seq = dna"ACGTACGT"
                mutable!(seq)
                rnaseq = convert(RNASequence, seq)
                @fact ismutable(rnaseq) --> true
                @fact_throws rnaseq[1] = DNA_C
                rnaseq[1] = RNA_C
                @fact seq --> dna"ACGTACGT"
                @fact rnaseq --> rna"CCGUACGU"
            end
        end

        context("SequenceNIterator") do
            function check_ns(T, seq)
                expected = Int[]
                for i in 1:length(seq)
                    if seq[i] == 'N'
                        push!(expected, i)
                    end
                end

                collect(npositions(T(seq))) == expected
            end

            reps = 10
            for len in [1, 10, 32, 1000, 10000, 100000]
                @fact all(Bool[check_ns(DNASequence, random_dna(len)) for _ in 1:reps]) --> true
                @fact all(Bool[check_ns(RNASequence, random_rna(len)) for _ in 1:reps]) --> true
            end

            dna_seq   = dna"ANANANA"
            dna_niter = npositions(dna_seq)
            @fact start(dna_niter) --> 2
            @fact next(dna_niter, 2) --> (2,4)
            @fact done(dna_niter, 6) --> false
            @fact done(dna_niter, 8) --> true


            ns = [2,4,6]
            for (i, n) in enumerate(dna_niter)
                @fact n --> ns[i]
            end

            rna_seq   = rna"ANANANA"
            rna_niter = npositions(rna_seq)
            @fact start(rna_niter) --> 2
            @fact next(rna_niter, 2) --> (2,4)
            @fact done(rna_niter, 6) --> false
            @fact done(rna_niter, 8) --> true

            ns = [2,4,6]
            for (i, n) in enumerate(rna_niter)
                @fact n --> ns[i]
            end

        end

        context("Kmer") do
            reps = 100
            context("Construction and Conversions") do
                # Check that kmers in strings survive round trip conversion:
                #   UInt64 → Kmer → UInt64
                function check_uint64_convertion(T::Type, n::UInt64, len::Int)
                    return convert(UInt64, convert(Kmer{T, len}, n)) === n
                end

                # Check that kmers in strings survive round trip conversion:
                #   String → Kmer → String
                function check_string_construction(T::Type, seq::AbstractString)
                    return convert(AbstractString, convert(Kmer{T}, seq)) == uppercase(seq)
                end

                # Check that dnakmers can be constructed from a DNASequence
                #   DNASequence → Kmer → DNASequence
                function check_dnasequence_construction(seq::DNASequence)
                    return convert(DNASequence, convert(DNAKmer, seq)) == seq
                end

                # Check that rnakmers can be constructed from a RNASequence
                #   RNASequence → Kmer → RNASequence
                function check_rnasequence_construction(seq::RNASequence)
                    return convert(RNASequence, convert(RNAKmer, seq)) == seq
                end

                # Check that kmers can be constructed from a NucleotideSequence
                #   NucleotideSequence → Kmer → NucleotideSequence
                function check_nucsequence_construction(seq::NucleotideSequence)
                    return convert(NucleotideSequence, convert(Kmer, seq)) == seq
                end

                # Check that kmers can be constructed from an array of nucleotides
                #   Vector{Nucleotide} → Kmer → Vector{Nucleotide}
                function check_nucarray_kmer{T <: Nucleotide}(seq::Vector{T})
                    return convert(AbstractString, [convert(Char, c) for c in seq]) ==
                           convert(AbstractString, kmer(seq...))
                end

                # Check that kmers in strings survive round trip conversion:
                #   String → NucleotideSequence → Kmer → NucleotideSequence → String
                function check_roundabout_construction(T::Type, seq::AbstractString)
                    return convert(AbstractString,
                               convert(NucleotideSequence{T},
                                   convert(Kmer,
                                       convert(NucleotideSequence{T}, seq)))) == uppercase(seq)
                end



                for len in [0, 1, 16, 32]
                    # UInt64 conversions
                    @fact all(Bool[check_uint64_convertion(DNANucleotide, rand(UInt64), len)
                              for _ in 1:reps]) --> true
                    @fact all(Bool[check_uint64_convertion(RNANucleotide, rand(UInt64), len)
                              for _ in 1:reps]) --> true

                    # String construction
                    @fact all(Bool[check_string_construction(DNANucleotide, random_dna_kmer(len))
                               for _ in 1:reps]) --> true
                    @fact all(Bool[check_string_construction(RNANucleotide, random_rna_kmer(len))
                               for _ in 1:reps]) --> true

                    # DNA/RNASequence Constructions
                    @fact all(Bool[check_dnasequence_construction(DNASequence(random_dna_kmer(len))) for _ in 1:reps]) --> true
                    @fact all(Bool[check_rnasequence_construction(RNASequence(random_rna_kmer(len))) for _ in 1:reps]) --> true

                    # NucleotideSequence Construction
                    @fact all(Bool[check_nucsequence_construction(DNASequence(random_dna_kmer(len))) for _ in 1:reps]) --> true
                    @fact all(Bool[check_nucsequence_construction(RNASequence(random_rna_kmer(len))) for _ in 1:reps]) --> true

                    # Construction from nucleotide arrays
                    if len > 0
                        @fact all(Bool[check_nucarray_kmer(random_dna_kmer_nucleotides(len))
                                   for _ in 1:reps]) --> true
                        @fact all(Bool[check_nucarray_kmer(random_rna_kmer_nucleotides(len))
                                   for _ in 1:reps]) --> true
                    end

                    # Roundabout conversions
                    @fact all(Bool[check_roundabout_construction(DNANucleotide, random_dna_kmer(len))
                               for _ in 1:reps]) --> true
                    @fact all(Bool[check_roundabout_construction(RNANucleotide, random_rna_kmer(len))
                               for _ in 1:reps]) --> true
                end

                @fact_throws kmer() # can't construct 0-mer using `kmer()`
                @fact_throws kmer(RNA_A, RNA_C, RNA_G, RNA_N, RNA_U) # no Ns in kmers
                @fact_throws kmer(DNA_A, DNA_C, DNA_G, DNA_N, DNA_T) # no Ns in kmers
                @fact_throws kmer(rna"ACGNU")# no Ns in kmers
                @fact_throws rnmakmer(rna"ACGNU")# no Ns in kmers
                @fact_throws kmer(dna"ACGNT") # no Ns in kmers
                @fact_throws dnmakmer(dna"ACGNT") # no Ns in kmers
                @fact_throws kmer(RNA_A, DNA_A) # no mixing of RNA and DNA
                @fact_throws kmer(random_rna(33)) # no kmer larger than 32nt
                @fact_throws kmer(random_dna(33)) # no kmer larger than 32nt
                @fact_throws kmer(RNA_A, RNA_C, RNA_G, RNA_U, # no kmer larger than 32nt
                                  RNA_A, RNA_C, RNA_G, RNA_U,
                                  RNA_A, RNA_C, RNA_G, RNA_U,
                                  RNA_A, RNA_C, RNA_G, RNA_U,
                                  RNA_A, RNA_C, RNA_G, RNA_U,
                                  RNA_A, RNA_C, RNA_G, RNA_U,
                                  RNA_A, RNA_C, RNA_G, RNA_U,
                                  RNA_A, RNA_C, RNA_G, RNA_U,
                                  RNA_A, RNA_C, RNA_G, RNA_U)
                @fact_throws kmer(DNA_A, DNA_C, DNA_G, DNA_T, # no kmer larger than 32nt
                                  DNA_A, DNA_C, DNA_G, DNA_T,
                                  DNA_A, DNA_C, DNA_G, DNA_T,
                                  DNA_A, DNA_C, DNA_G, DNA_T,
                                  DNA_A, DNA_C, DNA_G, DNA_T,
                                  DNA_A, DNA_C, DNA_G, DNA_T,
                                  DNA_A, DNA_C, DNA_G, DNA_T,
                                  DNA_A, DNA_C, DNA_G, DNA_T,
                                  DNA_A, DNA_C, DNA_G, DNA_T)

                context("From strings") do
                    @fact dnakmer("ACTG") --> convert(Kmer, DNASequence("ACTG"))
                    @fact rnakmer("ACUG") --> convert(Kmer, RNASequence("ACUG"))

                    # N is not allowed in Kmers
                    @fact_throws dnakmer("ACGTNACGT")
                    @fact_throws rnakmer("ACGUNACGU")
                end
            end

            context("Comparisons") do
                context("Equality") do
                    function check_seq_kmer_equality(len)
                        a = dnakmer(random_dna_kmer(len))
                        b = convert(DNASequence, a)
                        return a == b && b == a
                    end

                    for len in [1, 10, 32]
                        @fact all(Bool[check_seq_kmer_equality(len) for _ in 1:reps]) --> true
                    end

                    # True negatives
                    @fact (dnakmer("ACG") == rnakmer("ACG")) --> false
                    @fact (dnakmer("T")   == rnakmer("U"))   --> false
                    @fact (dnakmer("AC")  == dnakmer("AG"))  --> false
                    @fact (rnakmer("AC")  == rnakmer("AG"))  --> false

                    @fact (dnakmer("ACG") == rna"ACG") --> false
                    @fact (dnakmer("T")   == rna"U")   --> false
                    @fact (dnakmer("AC")  == dna"AG")  --> false
                    @fact (rnakmer("AC")  == rna"AG")  --> false

                    @fact (rna"ACG" == dnakmer("ACG")) --> false
                    @fact (rna"U"   == dnakmer("T"))   --> false
                    @fact (dna"AG"  == dnakmer("AC"))  --> false
                    @fact (rna"AG"  == rnakmer("AC"))  --> false
                end

                context("Inequality") do
                    for len in [1, 10, 32]
                        @fact isless(convert(DNAKmer{1}, UInt64(0)), convert(DNAKmer{1}, UInt64(1))) --> true
                        @fact isless(convert(DNAKmer{1}, UInt64(0)), convert(DNAKmer{1}, UInt64(0))) --> false
                        @fact isless(convert(DNAKmer{1}, UInt64(1)), convert(DNAKmer{1}, UInt64(0))) --> false

                        @fact isless(convert(RNAKmer{1}, UInt64(0)), convert(RNAKmer{1}, UInt64(1))) --> true
                        @fact isless(convert(RNAKmer{1}, UInt64(0)), convert(RNAKmer{1}, UInt64(0))) --> false
                        @fact isless(convert(RNAKmer{1}, UInt64(1)), convert(RNAKmer{1}, UInt64(0))) --> false
                    end
                end
            end

            context("Length") do
                for len in [0, 1, 16, 32]
                    @fact length(dnakmer(random_dna_kmer(len))) --> len
                    @fact length(rnakmer(random_rna_kmer(len))) --> len
                end
            end

            context("Access and Iterations") do
                dna_kmer = dnakmer("ACTG")
                rna_kmer = rnakmer("ACUG")

                context("Access DNA Kmer") do
                    @fact dna_kmer[1] --> DNA_A
                    @fact dna_kmer[2] --> DNA_C
                    @fact dna_kmer[3] --> DNA_T
                    @fact dna_kmer[4] --> DNA_G

                    # Access indexes out of bounds
                    @fact_throws dna_kmer[-1]
                    @fact_throws dna_kmer[0]
                    @fact_throws dna_kmer[5]
                    @fact_throws getindex(dna_kmer,-1)
                    @fact_throws getindex(dna_kmer, 0)
                    @fact_throws getindex(dna_kmer, 5)
                end

                context("Iteration through DNA Kmer") do
                    @fact start(dnakmer("ACTG"))  --> 1
                    @fact start(dnakmer(""))      --> 1

                    @fact next(dnakmer("ACTG"), 1) --> (DNA_A, 2)
                    @fact next(dnakmer("ACTG"), 4) --> (DNA_G, 5)

                    @fact done(dnakmer(""), 1)      --> true
                    @fact done(dnakmer("ACTG"), 1)  --> false
                    @fact done(dnakmer("ACTG"), 4)  --> false
                    @fact done(dnakmer("ACTG"), 5)  --> true
                    @fact done(dnakmer("ACTG"), -1) --> false


                    dna_kmer_vector = [DNA_A, DNA_C, DNA_T, DNA_G]
                    @fact all(Bool[nucleotide == dna_kmer_vector[i] for (i, nucleotide) in enumerate(dna_kmer)]) -->  true
                end

                context("Access RNA Kmer") do
                    @fact rna_kmer[1] --> RNA_A
                    @fact rna_kmer[2] --> RNA_C
                    @fact rna_kmer[3] --> RNA_U
                    @fact rna_kmer[4] --> RNA_G

                    # Access indexes out of bounds
                    @fact_throws rna_kmer[-1]
                    @fact_throws rna_kmer[0]
                    @fact_throws rna_kmer[5]
                    @fact_throws getindex(rna_kmer, -1)
                    @fact_throws getindex(rna_kmer, 0)
                    @fact_throws getindex(rna_kmer, 5)
                end

                context("Iteration through RNA Kmer") do
                    @fact start(rnakmer("ACUG"))  --> 1
                    @fact start(rnakmer(""))      --> 1

                    @fact next(rnakmer("ACUG"), 1) --> (RNA_A, 2)
                    @fact next(rnakmer("ACUG"), 4) --> (RNA_G, 5)

                    @fact done(rnakmer(""), 1)      --> true
                    @fact done(rnakmer("ACUG"), 1)  --> false
                    @fact done(rnakmer("ACUG"), 4)  --> false
                    @fact done(rnakmer("ACUG"), 5)  --> true
                    @fact done(rnakmer("ACUG"), -1) --> false

                    rna_kmer_vector = [RNA_A, RNA_C, RNA_U, RNA_G]
                    @fact all(Bool[nucleotide == rna_kmer_vector[i] for (i, nucleotide) in enumerate(rna_kmer)]) -->  true
                end
            end

            reps = 10
            context("Transformations") do
                context("Reversal") do
                    for len in [0, 1, 16, 32]
                        @fact all(Bool[check_reversal(DNAKmer, random_dna_kmer(len)) for _ in 1:reps]) --> true
                        @fact all(Bool[check_reversal(RNAKmer, random_rna_kmer(len)) for _ in 1:reps]) --> true
                    end
                end

                context("Complement") do
                    for len in [0, 1, 16, 32]
                        @fact all(Bool[check_dna_complement(DNAKmer, random_dna_kmer(len)) for _ in 1:reps]) --> true
                        @fact all(Bool[check_rna_complement(RNAKmer, random_rna_kmer(len)) for _ in 1:reps]) --> true
                    end
                end

                context("Reverse Complement") do
                    for len in [0, 1, 16, 32]
                        @fact all(Bool[check_dna_revcomp(DNAKmer, random_dna_kmer(len)) for _ in 1:reps]) --> true
                        @fact all(Bool[check_rna_revcomp(RNAKmer, random_rna_kmer(len)) for _ in 1:reps]) --> true
                    end
                end
            end

            context("Mismatches") do
                for len in [0, 1, 16, 32]
                    @fact all(Bool[check_mismatches(DNAKmer, random_dna_kmer(len), random_dna_kmer(len))
                              for _ in 1:reps]) --> true
                    @fact all(Bool[check_mismatches(RNAKmer, random_rna_kmer(len), random_rna_kmer(len))
                              for _ in 1:reps]) --> true
                end
            end
        end

        context("EachKmer") do
            function string_eachkmer(seq::AbstractString, k, step=1)
                kmers = AbstractString[]
                i = 1
                for i in 1:step:length(seq) - k + 1
                    subseq = seq[i:i + k - 1]
                    if !in('N', subseq)
                        push!(kmers, subseq)
                    end
                end
                return kmers
            end

            function check_eachkmer(T, seq::AbstractString, k, step=1)
                xs = [convert(AbstractString, x) for (i, x) in collect(each(Kmer{T, k}, NucleotideSequence{T}(seq), step))]
                ys = string_eachkmer(seq, k, step)
                return xs == ys
            end

            reps = 10
            len = 10000

            for k in [1, 3, 16, 32]
                @fact all(Bool[check_eachkmer(DNANucleotide, random_dna(len), k) for _ in 1:reps]) --> true
                @fact all(Bool[check_eachkmer(RNANucleotide, random_rna(len), k) for _ in 1:reps]) --> true
            end

            for k in [1, 3, 16, 32]
                @fact all(Bool[check_eachkmer(DNANucleotide, random_dna(len), k, 3) for _ in 1:reps]) --> true
                @fact all(Bool[check_eachkmer(RNANucleotide, random_rna(len), k, 3) for _ in 1:reps]) --> true
            end

            @fact isempty(collect(each(DNAKmer{1}, dna""))) --> true
            @fact isempty(collect(each(DNAKmer{1}, dna"NNNNNNNNNN"))) --> true
            @fact isempty(collect(each(DNAKmer{0}, dna"ACGT"))) --> true
            @fact_throws each(DNAKmer{-1}, dna"ACGT")
            @fact_throws each(DNAKmer{33}, dna"ACGT")
        end

        context("Nucleotide Counting") do
            function string_nucleotide_count(::Type{DNANucleotide}, seq::AbstractString)
                counts = Dict{DNANucleotide, Int}(
                    DNA_A => 0,
                    DNA_C => 0,
                    DNA_G => 0,
                    DNA_T => 0,
                    DNA_N => 0 )
                for c in seq
                    counts[convert(DNANucleotide, c)] += 1
                end

                return counts
            end

            function string_nucleotide_count(::Type{RNANucleotide}, seq::AbstractString)
                counts = Dict{RNANucleotide, Int}(
                    RNA_A => 0,
                    RNA_C => 0,
                    RNA_G => 0,
                    RNA_U => 0,
                    RNA_N => 0 )
                for c in seq
                    counts[convert(RNANucleotide, c)] += 1
                end

                return counts
            end

            function check_nucleotide_count(::Type{DNANucleotide}, seq::AbstractString)
                string_counts = string_nucleotide_count(DNANucleotide, seq)
                seq_counts = NucleotideCounts(DNASequence(seq))
                return string_counts[DNA_A] == seq_counts[DNA_A] &&
                       string_counts[DNA_C] == seq_counts[DNA_C] &&
                       string_counts[DNA_G] == seq_counts[DNA_G] &&
                       string_counts[DNA_T] == seq_counts[DNA_T] &&
                       string_counts[DNA_N] == seq_counts[DNA_N]
            end

            function check_nucleotide_count(::Type{RNANucleotide}, seq::AbstractString)
                string_counts = string_nucleotide_count(RNANucleotide, seq)
                seq_counts = NucleotideCounts(RNASequence(seq))
                return string_counts[RNA_A] == seq_counts[RNA_A] &&
                       string_counts[RNA_C] == seq_counts[RNA_C] &&
                       string_counts[RNA_G] == seq_counts[RNA_G] &&
                       string_counts[RNA_U] == seq_counts[RNA_U] &&
                       string_counts[RNA_N] == seq_counts[RNA_N]
            end

            function check_kmer_nucleotide_count(::Type{DNANucleotide}, seq::AbstractString)
                string_counts = string_nucleotide_count(DNANucleotide, seq)
                kmer_counts = NucleotideCounts(dnakmer(seq))
                return string_counts[DNA_A] == kmer_counts[DNA_A] &&
                       string_counts[DNA_C] == kmer_counts[DNA_C] &&
                       string_counts[DNA_G] == kmer_counts[DNA_G] &&
                       string_counts[DNA_T] == kmer_counts[DNA_T] &&
                       string_counts[DNA_N] == kmer_counts[DNA_N]
            end

            function check_kmer_nucleotide_count(::Type{RNANucleotide}, seq::AbstractString)
                string_counts = string_nucleotide_count(RNANucleotide, seq)
                kmer_counts = NucleotideCounts(rnakmer(seq))
                return string_counts[RNA_A] == kmer_counts[RNA_A] &&
                       string_counts[RNA_C] == kmer_counts[RNA_C] &&
                       string_counts[RNA_G] == kmer_counts[RNA_G] &&
                       string_counts[RNA_U] == kmer_counts[RNA_U] &&
                       string_counts[RNA_N] == kmer_counts[RNA_N]
            end

            reps = 10
            for len in [1, 10, 32, 1000, 10000, 100000]
                @fact all(Bool[check_nucleotide_count(DNANucleotide, random_dna(len))
                           for _ in 1:reps]) --> true
                @fact all(Bool[check_nucleotide_count(RNANucleotide, random_rna(len))
                           for _ in 1:reps]) --> true
            end

            for len in [1, 10, 32]
                @fact all(Bool[check_kmer_nucleotide_count(DNANucleotide, random_dna_kmer(len))
                           for _ in 1:reps]) --> true
                @fact all(Bool[check_kmer_nucleotide_count(RNANucleotide, random_rna_kmer(len))
                           for _ in 1:reps]) --> true
            end
        end

        context("Kmer Counting") do
            function string_kmer_count{T <: Nucleotide}(::Type{T}, seq::AbstractString, k, step)
                counts = Dict{Kmer{T, k}, Int}()
                for x in UInt64(0):UInt64(4^k-1)
                    counts[convert(Kmer{T, k}, x)] = 0
                end

                for i in 1:step:(length(seq)-k+1)
                    s = seq[i:i+k-1]
                    if 'N' in s
                        continue
                    end
                    counts[convert(Kmer{T, k}, s)] += 1
                end

                return counts
            end

            function check_kmer_count{T <: Nucleotide}(::Type{T}, seq::AbstractString, k, step)
                string_counts = string_kmer_count(T, seq, k, step)
                kmer_counts = KmerCounts{T, k}(convert(NucleotideSequence{T}, seq), step)
                for y in UInt64(0):UInt64(4^k-1)
                    x = convert(Kmer{T, k}, y)
                    if string_counts[x] != kmer_counts[x]
                        return false
                    end
                end
                return true
            end

            reps = 10
            for len in [1, 10, 32, 1000, 10000]
                for k in [1, 2, 5]
                    for step in [1, 3]
                        @fact all(Bool[check_kmer_count(DNANucleotide, random_dna(len), k, step)
                                   for _ in 1:reps]) --> true
                        @fact all(Bool[check_kmer_count(RNANucleotide, random_rna(len), k, step)
                                   for _ in 1:reps]) --> true
                    end
                end
            end
        end
    end
end

facts("Aminoacids") do
    context("AminoAcid Sequences") do
        reps = 10
        context("Construction") do
            # Non-aa characters should throw
            @fact_throws AminoAcidSequence("ATGHLMYZZACAGNM")

            # Check that sequences in strings survive round trip conversion:
            #   String → AminoAcidSequence → String
            function check_string_construction(seq::AbstractString)
                return convert(AbstractString, AminoAcidSequence(seq)) == uppercase(seq)
            end

            for len in [0, 1, 10, 32, 1000, 10000, 100000]
                @fact all(Bool[check_string_construction(random_aa(len)) for _ in 1:reps]) --> true
                @fact all(Bool[check_string_construction(lowercase(random_aa(len))) for _ in 1:reps]) --> true
            end

            # Check creation of empty
        end

        context("Conversion") do
            seq = aa"ARNDCQEGHILKMFPSTWYVX"
            @fact convert(AminoAcidSequence, [aa for aa in seq]) --> seq
            @fact convert(Vector{AminoAcid}, seq) --> [aa for aa in seq]

            @fact_throws convert(AminoAcidSequence, [convert(AminoAcid, UInt8(30))])
        end

        context("Concatenation") do
            function check_concatenation(n)
                chunks = [random_aa(rand(100:300)) for i in 1:n]
                parts = Any[]
                for i in 1:n
                    start = rand(1:length(chunks[i]))
                    stop = rand(start:length(chunks[i]))
                    push!(parts, start:stop)
                end

                str = string([chunk[parts[i]]
                              for (i, chunk) in enumerate(chunks)]...)

                seq = *([AminoAcidSequence(chunk)[parts[i]]
                         for (i, chunk) in enumerate(chunks)]...)

                return convert(AbstractString, seq) == uppercase(str)
            end

            @fact all(Bool[check_concatenation(rand(1:10)) for _ in 1:100]) --> true
        end

        context("Equality") do
            seq = aa"ARNDCQEGHILKMFPSTWYVX"
            @fact ==(seq, aa"ARNDCQEGHILKMFPSTWYVX") --> true
            @fact ==(seq, aa"ARNDCQEGHILKMFPSTWYXV") --> false
            @fact ==(seq, aa"ARNDCQEGHLKMFPSTWYVX")  --> false

            seq′ = aa"""
            ARNDCQEGHI
            LKMFPSTWYV
            X
            """
            @fact ==(seq, seq′) --> true
        end

        context("Repetition") do
            function check_repetition(n)
                chunk = random_aa(rand(100:300))
                start = rand(1:length(chunk))
                stop = rand(start:length(chunk))

                str = chunk[start:stop] ^ n
                seq = AminoAcidSequence(chunk)[start:stop] ^ n

                return convert(AbstractString, seq) == uppercase(str)
            end

            @fact all(Bool[check_repetition(rand(1:10)) for _ in 1:100]) --> true
        end

        context("Copy") do
            function check_copy(seq)
                return convert(AbstractString, copy(AminoAcidSequence(seq))) == uppercase(seq)
            end

            for len in [1, 10, 32, 1000, 10000, 100000]
                @fact all(Bool[check_copy(random_aa(len)) for _ in 1:reps]) --> true
            end
        end

        context("Subsequence Construction") do
            for len in [1, 10, 32, 1000, 10000, 100000]
                seq = random_aa(len)
                aaseq = AminoAcidSequence(seq)

                results = Bool[]
                for _ in 1:reps
                    part = random_interval(1, length(seq))
                    push!(results, seq[part] == convert(AbstractString, aaseq[part]))
                end
                @fact all(results) --> true
            end
        end

        context("Mutability") do
            seq = aa"ARNDCQEGHILKMFPSTWYVX"
            @fact ismutable(seq) --> false
            @fact_throws seq[1] = AA_C

            seq2 = seq[1:4]
            @fact ismutable(seq2) --> false

            mutable!(seq)
            @fact ismutable(seq) --> true
            seq[1] = AA_C
            @fact seq[1] --> AA_C
            @fact seq2[1] --> AA_A

            immutable!(seq)
            @fact_throws seq[1] = AA_A
        end
    end

    context("Parsers") do
        context("Valid Cases") do
            # case-insensitive and ignores spaces
            @fact parse(AminoAcid, "a") --> AA_A
            @fact parse(AminoAcid, "Ala") --> AA_A
            @fact parse(AminoAcid, "aLa ") --> AA_A
            @fact parse(AminoAcid, " alA ") --> AA_A
            @fact parse(AminoAcid, "\tAlA\n") --> AA_A
            @fact parse(AminoAcid, "x") --> AA_X
            @fact parse(AminoAcid, "X") --> AA_X
            aas = [
                ("A", "ALA", AA_A),
                ("R", "ARG", AA_R),
                ("N", "ASN", AA_N),
                ("D", "ASP", AA_D),
                ("C", "CYS", AA_C),
                ("E", "GLU", AA_E),
                ("Q", "GLN", AA_Q),
                ("G", "GLY", AA_G),
                ("H", "HIS", AA_H),
                ("I", "ILE", AA_I),
                ("L", "LEU", AA_L),
                ("K", "LYS", AA_K),
                ("M", "MET", AA_M),
                ("F", "PHE", AA_F),
                ("P", "PRO", AA_P),
                ("S", "SER", AA_S),
                ("T", "THR", AA_T),
                ("W", "TRP", AA_W),
                ("Y", "TYR", AA_Y),
                ("V", "VAL", AA_V),
            ]
            @fact length(aas) --> 20
            for (one, three, aa) in aas
                @fact parse(AminoAcid, one) --> aa
                @fact parse(AminoAcid, three) --> aa
            end
        end

        context("Invalid Cases") do
            @fact_throws parse(AminoAcid, "")
            @fact_throws parse(AminoAcid, "AL")
            @fact_throws parse(AminoAcid, "LA")
            @fact_throws parse(AminoAcid, "ALAA")
        end
    end
end

facts("Translation") do
    # crummy string translation to test against
    standard_genetic_code_dict = Dict{AbstractString, Char}(
        "AAA" => 'K', "AAC" => 'N', "AAG" => 'K', "AAU" => 'N',
        "ACA" => 'T', "ACC" => 'T', "ACG" => 'T', "ACU" => 'T',
        "AGA" => 'R', "AGC" => 'S', "AGG" => 'R', "AGU" => 'S',
        "AUA" => 'I', "AUC" => 'I', "AUG" => 'M', "AUU" => 'I',
        "CAA" => 'Q', "CAC" => 'H', "CAG" => 'Q', "CAU" => 'H',
        "CCA" => 'P', "CCC" => 'P', "CCG" => 'P', "CCU" => 'P',
        "CGA" => 'R', "CGC" => 'R', "CGG" => 'R', "CGU" => 'R',
        "CUA" => 'L', "CUC" => 'L', "CUG" => 'L', "CUU" => 'L',
        "GAA" => 'E', "GAC" => 'D', "GAG" => 'E', "GAU" => 'D',
        "GCA" => 'A', "GCC" => 'A', "GCG" => 'A', "GCU" => 'A',
        "GGA" => 'G', "GGC" => 'G', "GGG" => 'G', "GGU" => 'G',
        "GUA" => 'V', "GUC" => 'V', "GUG" => 'V', "GUU" => 'V',
        "UAA" => '*', "UAC" => 'Y', "UAG" => '*', "UAU" => 'Y',
        "UCA" => 'S', "UCC" => 'S', "UCG" => 'S', "UCU" => 'S',
        "UGA" => '*', "UGC" => 'C', "UGG" => 'W', "UGU" => 'C',
        "UUA" => 'L', "UUC" => 'F', "UUG" => 'L', "UUU" => 'F',

        # translatable ambiguities in the standard code
        "CUN" => 'L', "CCN" => 'P', "CGN" => 'R', "ACN" => 'T',
        "GUN" => 'V', "GCN" => 'A', "GGN" => 'G', "UCN" => 'S'

    )

    function string_translate(seq::AbstractString)
        @assert length(seq) % 3 == 0
        aaseq = Array(Char, div(length(seq), 3))
        for i in 1:3:length(seq) - 3 + 1
            aaseq[div(i, 3) + 1] = standard_genetic_code_dict[seq[i:i+2]]
        end
        return convert(AbstractString, aaseq)
    end

    function check_translate(seq::AbstractString)
        return string_translate(seq) == convert(AbstractString, translate(RNASequence(seq)))
    end

    reps = 10
    #for len in [1, 10, 32, 1000, 10000, 100000]
    for len in [1, 10, 32]
        @fact all(Bool[check_translate(random_translatable_rna(len)) for _ in 1:reps]) --> true
    end

    @fact_throws translate(dna"ACGTACGTA") # can't translate DNA
    @fact_throws translate(rna"ACGUACGU")  # can't translate non-multiples of three
    @fact_throws translate(rna"ACGUACGNU") # can't translate N
end


facts("Sequence Parsing") do
    context("FASTA Parsing") do
        get_bio_fmt_specimens()

        function check_fasta_parse(filename)
            # Reading from a stream
            for seqrec in open(filename, FASTA)
            end

            # Reading from a memory mapped file
            for seqrec in open(filename, FASTA, memory_map=true)
            end

            # Check round trip
            output = IOBuffer()
            expected_entries = Any[]
            for seqrec in open(filename, FASTA)
                write(output, seqrec)
                push!(expected_entries, seqrec)
            end

            read_entries = Any[]
            for seqrec in open(takebuf_array(output), FASTA)
                push!(read_entries, seqrec)
            end

            return expected_entries == read_entries
        end

        path = Pkg.dir("Bio", "test", "BioFmtSpecimens", "FASTA")
        for specimen in YAML.load_file(joinpath(path, "index.yml"))
            tags = specimen["tags"]
            valid = get(specimen, "valid", true)
            # currently unsupported features
            if contains(tags, "gaps") || contains(tags, "comments") || contains(tags, "ambiguity")
                continue
            end
            if valid
                @fact check_fasta_parse(joinpath(path, specimen["filename"])) --> true
            else
                @fact_throws check_fasta_parse(joinpath(path, specimen["filename"]))
            end
        end
    end

    context("FASTQ Parsing") do
        get_bio_fmt_specimens()

        function check_fastq_parse(filename)
            # Reading from a stream
            for seqrec in open(filename, FASTQ)
            end

            # Reading from a memory mapped file
            for seqrec in open(filename, FASTQ, memory_map=true)
            end

            # Check round trip
            output = IOBuffer()
            expected_entries = Any[]
            for seqrec in open(filename, FASTQ)
                write(output, seqrec)
                push!(expected_entries, seqrec)
            end

            read_entries = Any[]
            for seqrec in open(takebuf_array(output), FASTQ)
                push!(read_entries, seqrec)
            end

            return expected_entries == read_entries
        end

        path = Pkg.dir("Bio", "test", "BioFmtSpecimens", "FASTQ")
        for specimen in YAML.load_file(joinpath(path, "index.yml"))
            tags = get(specimen, "tags", "")
            valid = get(specimen, "valid", true)
            # currently unsupported features
            if contains(tags, "rna") || contains(tags, "gaps") ||
               contains(tags, "comments") || contains(tags, "ambiguity")
                continue
            end
            if valid
                @fact check_fastq_parse(joinpath(path, specimen["filename"])) --> true
            else
                @fact_throws check_fastq_parse(joinpath(path, specimen["filename"]))
            end
        end
    end
end

facts("Quality scores") do
    using Bio.Seq: encode_quality_string!,
                   encode_quality_string,
                   decode_quality_string!,
                   decode_quality_string,
                   ILLUMINA13_QUAL_ENCODING,
                   ILLUMINA15_QUAL_ENCODING,
                   ILLUMINA18_QUAL_ENCODING,
                   SANGER_QUAL_ENCODING,
                   SOLEXA_QUAL_ENCODING

    context("Decoding PHRED scores") do

        function test_decode(encoding, values, expected)
            result = Array(Int8, length(expected))
            decode_quality_string!(encoding, values, result, 1,
                                   length(result))
            @fact result --> expected

            # Without start and end does the entire string
            decode_quality_string!(encoding, values, result)
            @fact result --> expected

            # Non in-place version of the above
            result = decode_quality_string(encoding, values, 1, length(result))
            @fact result --> expected
            resutlt = decode_quality_string(encoding, values)
            @fact result --> expected
        end

        test_decode(SANGER_QUAL_ENCODING,
                    UInt8['!', '#', '$', '%', '&', 'I', '~'],
                    Int8[0, 2, 3, 4, 5, 40, 93])

        test_decode(SOLEXA_QUAL_ENCODING,
                    UInt8[';', 'B', 'C', 'D', 'E', 'h', '~'],
                    Int8[-5, 2, 3, 4, 5, 40, 62])

        test_decode(ILLUMINA13_QUAL_ENCODING,
                    UInt8['@', 'B', 'C', 'D', 'E', 'h', '~'],
                    Int8[0, 2, 3, 4, 5, 40, 62])

        test_decode(ILLUMINA15_QUAL_ENCODING,
                    UInt8['C', 'D', 'E', 'h', '~'],
                    Int8[3, 4, 5, 40, 62])

        test_decode(ILLUMINA18_QUAL_ENCODING,
                    UInt8['!', '#', '$', '%', '&', 'I', '~'],
                    Int8[0, 2, 3, 4, 5, 40, 93])
    end

    context("Encoding PHRED scores") do

        function test_encode(encoding, values, expected)
            # With start & end
            result = Array(UInt8, length(expected))
            encode_quality_string!(encoding, values, result, 1,
                                   length(result))
            @fact result --> expected
            # Without start & end means the entire length
            encode_quality_string!(encoding, values, result)
            @fact result --> expected

            result = encode_quality_string(encoding, values, 1, length(result))
            @fact result --> expected
            result = encode_quality_string(encoding, values)
            @fact result --> expected
        end

        test_encode(SANGER_QUAL_ENCODING,
                    Int8[0, 2, 3, 4, 5, 40, 93],
                    UInt8['!', '#', '$', '%', '&', 'I', '~'])

        test_encode(SOLEXA_QUAL_ENCODING,
                    Int8[-5, 2, 3, 4, 5, 40, 62],
                    UInt8[';', 'B', 'C', 'D', 'E', 'h', '~'])

        test_encode(ILLUMINA13_QUAL_ENCODING,
                    Int8[0, 2, 3, 4, 5, 40, 62],
                    UInt8['@', 'B', 'C', 'D', 'E', 'h', '~'])

        test_encode(ILLUMINA15_QUAL_ENCODING,
                    Int8[3, 4, 5, 40, 62],
                    UInt8['C', 'D', 'E', 'h', '~'])

        test_encode(ILLUMINA18_QUAL_ENCODING,
                    Int8[0, 2, 3, 4, 5, 40, 93],
                    UInt8['!', '#', '$', '%', '&', 'I', '~'])
    end

    context("Sequence Writing") do
        context("FASTA writing") do
            dna_seq1 = random_dna(79 * 3) # full lines
            dna_seq2 = random_dna(50)     # short line
            seq_name1 = "Sequence 1"
            seq_name2 = "Sequence 2"
            seq_description1 = "Description 1"
            seq_description2 = "Description 2"

            wrapped_seq1 = join([dna_seq1[1:79], dna_seq1[80:158], dna_seq1[159:end]], "\n")

            expected_seq1 = string(">", seq_name1, " ", seq_description1, "\n", wrapped_seq1, "\n")
            expected_seq2 = string(">", seq_name2, " ", seq_description2, "\n", dna_seq2, "\n")
            expected = string(expected_seq1, expected_seq2)

            fasta_seq1 = Seq.FASTADNASeqRecord(
                seq_name1, DNASequence(dna_seq1), Seq.FASTAMetadata(seq_description1))
            fasta_seq2 = Seq.FASTADNASeqRecord(
                seq_name2, DNASequence(dna_seq2), Seq.FASTAMetadata(seq_description2))
            sequences = [fasta_seq1, fasta_seq2]

            output = IOBuffer()
            for seq in sequences
                write(output, seq)
            end
            @fact takebuf_string(output) --> expected
        end
    end
end

end # TestSeq
