--- /dev/null
+import math
+from CodonUsageIndices import SharpEcoliIndex
+from Bio import SeqIO # To parse a FASTA file
+
+CodonsDict = {'TTT':0, 'TTC':0, 'TTA':0, 'TTG':0, 'CTT':0,
+'CTC':0, 'CTA':0, 'CTG':0, 'ATT':0, 'ATC':0,
+'ATA':0, 'ATG':0, 'GTT':0, 'GTC':0, 'GTA':0,
+'GTG':0, 'TAT':0, 'TAC':0, 'TAA':0, 'TAG':0,
+'CAT':0, 'CAC':0, 'CAA':0, 'CAG':0, 'AAT':0,
+'AAC':0, 'AAA':0, 'AAG':0, 'GAT':0, 'GAC':0,
+'GAA':0, 'GAG':0, 'TCT':0, 'TCC':0, 'TCA':0,
+'TCG':0, 'CCT':0, 'CCC':0, 'CCA':0, 'CCG':0,
+'ACT':0, 'ACC':0, 'ACA':0, 'ACG':0, 'GCT':0,
+'GCC':0, 'GCA':0, 'GCG':0, 'TGT':0, 'TGC':0,
+'TGA':0, 'TGG':0, 'CGT':0, 'CGC':0, 'CGA':0,
+'CGG':0, 'AGT':0, 'AGC':0, 'AGA':0, 'AGG':0,
+'GGT':0, 'GGC':0, 'GGA':0, 'GGG':0}
+
+
+# this dictionary is used to know which codons encode the same AA.
+SynonymousCodons = {'CYS': ['TGT', 'TGC'], 'ASP': ['GAT', 'GAC'],
+'SER': ['TCT', 'TCG', 'TCA', 'TCC', 'AGC', 'AGT'],
+'GLN': ['CAA', 'CAG'], 'MET': ['ATG'], 'ASN': ['AAC', 'AAT'],
+'PRO': ['CCT', 'CCG', 'CCA', 'CCC'], 'LYS': ['AAG', 'AAA'],
+'STOP': ['TAG', 'TGA', 'TAA'], 'THR': ['ACC', 'ACA', 'ACG', 'ACT'],
+'PHE': ['TTT', 'TTC'], 'ALA': ['GCA', 'GCC', 'GCG', 'GCT'],
+'GLY': ['GGT', 'GGG', 'GGA', 'GGC'], 'ILE': ['ATC', 'ATA', 'ATT'],
+'LEU': ['TTA', 'TTG', 'CTC', 'CTT', 'CTG', 'CTA'], 'HIS': ['CAT', 'CAC'],
+'ARG': ['CGA', 'CGC', 'CGG', 'CGT', 'AGG', 'AGA'], 'TRP': ['TGG'],
+'VAL': ['GTA', 'GTC', 'GTG', 'GTT'], 'GLU': ['GAG', 'GAA'], 'TYR': ['TAT', 'TAC']}
+
+
+class CodonAdaptationIndex:
+ """A codon adaptaion index (CAI) implementation.
+
+ This class implements the codon adaptaion index (CAI) described by Sharp and
+ Li (Nucleic Acids Res. 1987 Feb 11;15(3):1281-95).
+
+ methods:
+
+ set_cai_index(Index):
+
+ This method sets-up an index to be used when calculating CAI for a gene.
+ Just pass a dictionary similar to the SharpEcoliIndex in CodonUsageIndices
+ module.
+
+ generate_index(FastaFile):
+
+ This method takes a location of a FastaFile and generates an index. This
+ index can later be used to calculate CAI of a gene.
+
+ cai_for_gene(DNAsequence):
+
+ This method uses the Index (either the one you set or the one you generated)
+ and returns the CAI for the DNA sequence.
+
+ print_index():
+ This method prints out the index you used.
+
+ NOTE - This implementation does not currently cope with alternative genetic
+ codes, only the synonymous codons in the standard table are considered.
+ """
+ def __init__(self):
+ self.index = {}
+ self.codon_count={}
+
+ # use this method with predefined CAI index
+ def set_cai_index(self, Index):
+ self.index = Index
+
+ def generate_index(self, FastaFile):
+ """Generate a codon usage index from a FASTA file of CDS sequences.
+
+ This method takes a location of a Fasta file containing CDS sequences
+ (which must all have a whole number of codons) and generates a codon
+ usage index. This index can later be used to calculate CAI of a gene.
+ """
+ # first make sure i am not overwriting an existing index:
+ if self.index != {} or self.codon_count!={}:
+ raise ValueError("an index has already been set or a codon count has been done. cannot overwrite either.")
+ # count codon occurances in the file.
+ self._count_codons(FastaFile)
+
+ # now to calculate the index we first need to sum the number of times
+ # synonymous codons were used all together.
+ for AA in SynonymousCodons.keys():
+ Sum=0.0
+ RCSU=[] # RCSU values are equal to CodonCount/((1/num of synonymous codons) * sum of all synonymous codons)
+
+ for codon in SynonymousCodons[AA]:
+ Sum += self.codon_count[codon]
+ # calculate the RSCU value for each of the codons
+ for codon in SynonymousCodons[AA]:
+ RCSU.append(self.codon_count[codon]/((1.0/len(SynonymousCodons[AA]))*Sum))
+ # now generate the index W=RCSUi/RCSUmax:
+ RCSUmax = max(RCSU)
+ for i in range(len(SynonymousCodons[AA])):
+ self.index[SynonymousCodons[AA][i]]= RCSU[i]/RCSUmax
+
+
+ def cai_for_gene(self, DNAsequence):
+ """Calculate the CAI (float) for the provided DNA sequence (string).
+
+ This method uses the Index (either the one you set or the one you generated)
+ and returns the CAI for the DNA sequence.
+ """
+ caiValue = 0
+ LengthForCai = 0
+ # if no index is set or generated, the default SharpEcoliIndex will be used.
+ if self.index=={}:
+ self.set_cai_index(SharpEcoliIndex)
+
+ if DNAsequence.islower():
+ DNAsequence = DNAsequence.upper()
+ for i in range (0,len(DNAsequence),3):
+ codon = DNAsequence[i:i+3]
+ if codon in self.index:
+ if codon!='ATG' and codon!= 'TGG': #these two codons are always one, exclude them.
+ caiValue += math.log(self.index[codon])
+ LengthForCai += 1
+ elif codon not in ['TGA','TAA', 'TAG']: # some indices you will use may not include stop codons.
+ raise TypeError("illegal codon in sequence: %s.\n%s" % (codon, self.index))
+ return math.exp(caiValue*(1.0/(LengthForCai-1)))
+
+ def _count_codons(self, FastaFile):
+ handle = open(FastaFile, 'r')
+
+ # make the codon dictionary local
+ self.codon_count = CodonsDict.copy()
+
+ # iterate over sequence and count all the codons in the FastaFile.
+ for cur_record in SeqIO.parse(handle, "fasta") :
+ # make sure the sequence is lower case
+ if str(cur_record.seq).islower():
+ DNAsequence = str(cur_record.seq).upper()
+ else:
+ DNAsequence = str(cur_record.seq)
+ for i in range(0,len(DNAsequence),3):
+ codon = DNAsequence[i:i+3]
+ if codon in self.codon_count:
+ self.codon_count[codon] += 1
+ else:
+ raise TypeError("illegal codon %s in gene: %s" % (codon, cur_record.id))
+ handle.close()
+
+ # this just gives the index when the objects is printed.
+ def print_index (self):
+ """This method prints out the index you used."""
+ X=self.index.keys()
+ X.sort()
+ for i in X:
+ print "%s\t%.3f" %(i, self.index[i])
+
git://source.jalview.org / jabaws.git / blobdiff