<commit_before>module TestSeq

using Compat
using FactCheck
using Bio
using Bio.Seq

# 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(String, 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(String, 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("DNA 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("Sequences") do
        function dna_complement(seq::String)
            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(String, seqc)
        end

        function rna_complement(seq::String)
            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(String, seqc)
        end

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

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

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

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

        function check_rna_revcomp(T, seq)
            return reverse(rna_complement(seq)) ==
                convert(String, 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::String)
                        return convert(String, NucleotideSequence{T}(seq)) == uppercase(seq)
                    end

                    for len in [0, 1, 10, 32, 1000, 10000, 100000]
                        @fact all([check_string_construction(DNANucleotide, random_dna(len)) for _ in 1:reps]) => true
                        @fact all([check_string_construction(RNANucleotide, random_rna(len)) for _ in 1:reps]) => true
                        @fact all([check_string_construction(DNANucleotide, lowercase(random_dna(len))) for _ in 1:reps]) => true
                        @fact all([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
            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([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
            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(String, copy(NucleotideSequence{T}(seq))) == seq
                end

                for len in [1, 10, 32, 1000, 10000, 100000]
                    @fact all([check_copy(DNANucleotide, random_dna(len)) for _ in 1:reps]) => true
                    @fact all([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") => 0
                    @fact start(dna"")      => 0

                    @fact next(dna"ACTGN", 1) => (DNA_C, 2)
                    @fact next(dna"ACTGN", 4) => (DNA_N, 5)
                    @fact_throws next(dna"ACTGN", 5)
                    @fact_throws next(dna"ACTGN", -1)

                    @fact done(dna"", 1)       => true
                    @fact done(dna"ACTGN", 1)  => false
                    @fact done(dna"ACTGN", 4)  => false
                    @fact done(dna"ACTGN", 5)  => true
                    @fact done(dna"ACTGN", -1) => false


                    dna_vector = [DNA_A, DNA_C, DNA_T, DNA_G]
                    @fact all([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") => 0
                    @fact start(rna"")      => 0

                    @fact next(rna"ACUGN", 1) => (RNA_C, 2)
                    @fact next(rna"ACUGN", 4) => (RNA_N, 5)
                    @fact_throws next(rna"ACUGN", 5)
                    @fact_throws next(rna"ACUGN", -1)

                    @fact done(rna"", 1)       => true
                    @fact done(rna"ACUGN", 1)  => false
                    @fact done(rna"ACUGN", 4)  => false
                    @fact done(rna"ACUGN", 5)  => true
                    @fact done(rna"ACUGN", -1) => false

                    # Iteration through RNA Sequence

                    rna_vector = [RNA_A, RNA_C, RNA_U, RNA_G]
                    @fact all([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(String, 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(String, 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
                end
            end

            context("Transformations") do
                context("Reversal") do
                    for len in [0, 1, 10, 32, 1000, 10000, 100000]
                        @fact all([check_reversal(DNASequence, random_dna(len)) for _ in 1:reps]) => true
                        @fact all([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([check_dna_complement(DNASequence, random_dna(len)) for _ in 1:reps]) => true
                        @fact all([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([check_dna_revcomp(DNASequence, random_dna(len)) for _ in 1:reps]) => true
                        @fact all([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([check_mismatches(DNASequence, random_dna(len), random_dna(len))
                               for _ in 1:reps]) => true
                    @fact all([check_mismatches(RNASequence, random_rna(len), random_rna(len))
                               for _ in 1:reps]) => true
                end
            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([check_ns(DNASequence, random_dna(len)) for _ in 1:reps]) => true
                @fact all([check_ns(RNASequence, random_rna(len)) for _ in 1:reps]) => true
            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::String)
                    return convert(String, 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(String, [convert(Char, c) for c in seq]) ==
                           convert(String, kmer(seq...))
                end

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



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

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

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

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

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

                    # Roundabout conversions
                    @fact all([check_roundabout_construction(DNANucleotide, random_dna_kmer(len))
                               for _ in 1:reps]) => true
                    @fact all([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([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
                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
                    # 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([nucleotide == dna_kmer_vector[i] for (i, nucleotide) in enumerate(dna_kmer)]) =>  true
                end

                context("Access RNA Kmer") do
                    # 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
<commit_msg>More tests for Kmers<commit_after>module TestSeq

using Compat
using FactCheck
using Bio
using Bio.Seq

# 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(String, 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(String, 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("DNA 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("Sequences") do
        function dna_complement(seq::String)
            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(String, seqc)
        end

        function rna_complement(seq::String)
            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(String, seqc)
        end

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

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

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

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

        function check_rna_revcomp(T, seq)
            return reverse(rna_complement(seq)) ==
                convert(String, 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::String)
                        return convert(String, NucleotideSequence{T}(seq)) == uppercase(seq)
                    end

                    for len in [0, 1, 10, 32, 1000, 10000, 100000]
                        @fact all([check_string_construction(DNANucleotide, random_dna(len)) for _ in 1:reps]) => true
                        @fact all([check_string_construction(RNANucleotide, random_rna(len)) for _ in 1:reps]) => true
                        @fact all([check_string_construction(DNANucleotide, lowercase(random_dna(len))) for _ in 1:reps]) => true
                        @fact all([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
            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([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
            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(String, copy(NucleotideSequence{T}(seq))) == seq
                end

                for len in [1, 10, 32, 1000, 10000, 100000]
                    @fact all([check_copy(DNANucleotide, random_dna(len)) for _ in 1:reps]) => true
                    @fact all([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") => 0
                    @fact start(dna"")      => 0

                    @fact next(dna"ACTGN", 1) => (DNA_C, 2)
                    @fact next(dna"ACTGN", 4) => (DNA_N, 5)
                    @fact_throws next(dna"ACTGN", 5)
                    @fact_throws next(dna"ACTGN", -1)

                    @fact done(dna"", 1)       => true
                    @fact done(dna"ACTGN", 1)  => false
                    @fact done(dna"ACTGN", 4)  => false
                    @fact done(dna"ACTGN", 5)  => true
                    @fact done(dna"ACTGN", -1) => false


                    dna_vector = [DNA_A, DNA_C, DNA_T, DNA_G]
                    @fact all([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") => 0
                    @fact start(rna"")      => 0

                    @fact next(rna"ACUGN", 1) => (RNA_C, 2)
                    @fact next(rna"ACUGN", 4) => (RNA_N, 5)
                    @fact_throws next(rna"ACUGN", 5)
                    @fact_throws next(rna"ACUGN", -1)

                    @fact done(rna"", 1)       => true
                    @fact done(rna"ACUGN", 1)  => false
                    @fact done(rna"ACUGN", 4)  => false
                    @fact done(rna"ACUGN", 5)  => true
                    @fact done(rna"ACUGN", -1) => false

                    # Iteration through RNA Sequence

                    rna_vector = [RNA_A, RNA_C, RNA_U, RNA_G]
                    @fact all([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(String, 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(String, 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
                end
            end

            context("Transformations") do
                context("Reversal") do
                    for len in [0, 1, 10, 32, 1000, 10000, 100000]
                        @fact all([check_reversal(DNASequence, random_dna(len)) for _ in 1:reps]) => true
                        @fact all([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([check_dna_complement(DNASequence, random_dna(len)) for _ in 1:reps]) => true
                        @fact all([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([check_dna_revcomp(DNASequence, random_dna(len)) for _ in 1:reps]) => true
                        @fact all([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([check_mismatches(DNASequence, random_dna(len), random_dna(len))
                               for _ in 1:reps]) => true
                    @fact all([check_mismatches(RNASequence, random_rna(len), random_rna(len))
                               for _ in 1:reps]) => true
                end
            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([check_ns(DNASequence, random_dna(len)) for _ in 1:reps]) => true
                @fact all([check_ns(RNASequence, random_rna(len)) for _ in 1:reps]) => true
            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::String)
                    return convert(String, 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(String, [convert(Char, c) for c in seq]) ==
                           convert(String, kmer(seq...))
                end

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



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

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

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

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

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

                    # Roundabout conversions
                    @fact all([check_roundabout_construction(DNANucleotide, random_dna_kmer(len))
                               for _ in 1:reps]) => true
                    @fact all([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([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([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

{'fraction_non_alphanumeric': 0.11278207279399081, 'fraction_numerical': 0.021860151484263513, 'mean_word_length': 1.4321128909229597, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 363, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 64, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6160655', 'n_tokens_mistral': 21409, 'n_tokens_neox': 18212, 'n_words': 4489}