- Contribution -

	

	  Usage: julia generate-codon-table.jl FOLDER
	
	
Where FOLDER is a path to a folder containing sequence.fasta (the DNA sequence of the gene) and sequence.gb (GenBank file containing translations).
	
Download
      

	
using GenomicAnnotations
using Dictionaries

folder = ARGS[1]

# Amino acid => [Count, [Associated codons]]
codonMap = Dictionary(Dict([
    'F'=>[0, ["TTT", "TTC"]],
    'L'=>[0, ["TTA", "TTG","CTT", "CTC", "CTA", "CTG"]],
    'S'=>[0, ["TCT", "TCC", "TCA", "TCG", "AGT", "AGC"]],
    'Y'=>[0, ["TAT", "TAC"]],
    '*'=>[0, ["TAA", "TAG", "TGA"]], #Stop codon
    'C'=>[0, ["TGT", "TGC"]],
    'W'=>[0, ["TGG"]],
    
    'P'=>[0, ["CCT", "CCC", "CCA", "CCG"]],
    'H'=>[0, ["CAT", "CAC"]],
    'Q'=>[0, ["CAA", "CAG"]],
    'R'=>[0, ["CGT", "CGC", "CGA", "CGG", "AGA", "AGG"]],

    'I'=>[0, ["ATT", "ATC", "ATA"]],
    'M'=>[0, ["ATG"]],
    'T'=>[0, ["ACT", "ACC", "ACA", "ACG"]],
    'N'=>[0, ["AAT", "AAC"]],
    'K'=>[0, ["AAA", "AAG"]],

    'V'=>[0, ["GTT", "GTC", "GTA", "GTG"]],
    'A'=>[0, ["GCT", "GCC", "GCA", "GCG"]],
    'D'=>[0, ["GAT", "GAC"]],
    'E'=>[0, ["GAA", "GAG"]],
    'G'=>[0, ["GGT", "GGC", "GGA", "GGG"]],
]))

# [Codon, probability]
probabilities = Dictionary(Dict([
    ["TTT", 0.0], ["TTC", 0.0], ["TTA", 0.0], ["TTG", 0.0],
    ["TCT", 0.0], ["TCC", 0.0], ["TCA", 0.0], ["TCG", 0.0],
    ["TAT", 0.0], ["TAC", 0.0], ["TAA", 0.0], ["TAG", 0.0],
    ["TGT", 0.0], ["TGC", 0.0], ["TGA", 0.0], ["TGG", 0.0],
    
    ["CTT", 0.0], ["CTC", 0.0], ["CTA", 0.0], ["CTG", 0.0],
    ["CCT", 0.0], ["CCC", 0.0], ["CCA", 0.0], ["CCG", 0.0],
    ["CAT", 0.0], ["CAC", 0.0], ["CAA", 0.0], ["CAG", 0.0],
    ["CGT", 0.0], ["CGC", 0.0], ["CGA", 0.0], ["CGG", 0.0],

    ["ATT", 0.0], ["ATC", 0.0], ["ATA", 0.0], ["ATG", 0.0],
    ["ACT", 0.0], ["ACC", 0.0], ["ACA", 0.0], ["ACG", 0.0],
    ["AAT", 0.0], ["AAC", 0.0], ["AAA", 0.0], ["AAG", 0.0],
    ["AGT", 0.0], ["AGC", 0.0], ["AGA", 0.0], ["AGG", 0.0],

    ["GTT", 0.0], ["GTC", 0.0], ["GTA", 0.0], ["GTG", 0.0],
    ["GCT", 0.0], ["GCC", 0.0], ["GCA", 0.0], ["GCG", 0.0],
    ["GAT", 0.0], ["GAC", 0.0], ["GAA", 0.0], ["GAG", 0.0],
    ["GGT", 0.0], ["GGC", 0.0], ["GGA", 0.0], ["GGG", 0.0],
]))

# Counts occurence of each codon in the genome file
function countCodons!()
    total = 0
    genomeFile = open("./$(folder)/sequence.fasta", "r")
    gbdata = readgbk("./$(folder)/sequence.gb")

    # Import sequence.fasta as a string
    genomeStr = ""
    for line in readlines(genomeFile)
        if (first(line, 1) != ">")              # Ignore > lines
            genomeStr = string(genomeStr, line) # Append line to string
        end
    end
    
    # For each chromosome in gbdata
    for chr in gbdata
        # For each gene
        for gene in @genes(chr, !ismissing(:translation))
            # Split DNA into codons, and assign relevant amino acid
            dna = String.(Iterators.partition(genomeStr[locus(gene).position], 3)) # Split into 3s
            for i in 1:length(dna)  # For each codon
                total += 1          # Increment codon counter
                codon = dna[i]      # Get codon from DNA sequence

                # Check codon valid
                if (i < length(gene.translation) || codon in keys(probabilities))
                    # Lookup amino acid from codon table or translation sequence
                    aminoacid = i > length(gene.translation) ?
                        getAminoAcid(codon) : gene.translation[i] # Use translation for efficiency if available

                    # If mismatched amino acid, log warning
                    if (i < length(gene.translation))
                        if (aminoacid != gene.translation[i])
                            println("WARNING: $(aminoacid) != $(gene.translation[i]) for codon $(codon) in $(chr.name)")
                        end
                    end

                    if (aminoacid != 'U')           # Ignore artificial amino acids
                        probabilities[codon] += 1   # Increment codon count
                        # codonMap[aminoacid][1] += 1 # Increment amino acid count
                    end
                end
            end
        end
    end
    
    close(genomeFile)
    println("Total number of codons: ", total)
end

function getAminoAcidValues(codon)
    filter(((aminoacid,aminoacidvalues),) -> codon in aminoacidvalues[2],
           pairs(codonMap))
end

function getAminoAcid(codon)
    getAminoAcidValues(codon).values[1][1]
end

function getAminoAcidCount(codon)
    getAminoAcidValues(codon).values[1][2][1]
end

function calculateAminoAcidCounts!()
    # For each codon associated with each amino acid
    for aminoacid in keys(codonMap)
        for codon in codonMap[aminoacid][2]
            # Add the count to the total
            codonMap[aminoacid][1] += probabilities[codon]
        end
    end
end

function calculateProbabilities!()
    calculateAminoAcidCounts!()
    
    # For each codon, divide the codon count by the amino acid count
    map(((codon,p),) -> p / getAminoAcidCount(codon),
         pairs(probabilities))
end

countCodons!()

println("Codon | Frequency")
display(calculateProbabilities!())