Package Bio :: Package SeqUtils :: Module CodonUsage
[hide private]
[frames] | no frames]

Source Code for Module Bio.SeqUtils.CodonUsage

  1  import math 
  2  from CodonUsageIndices import SharpEcoliIndex 
  3  from Bio import SeqIO # To parse a FASTA file 
  4   
  5  CodonsDict = {'TTT':0, 'TTC':0, 'TTA':0, 'TTG':0, 'CTT':0,  
  6  'CTC':0, 'CTA':0, 'CTG':0, 'ATT':0, 'ATC':0,  
  7  'ATA':0, 'ATG':0, 'GTT':0, 'GTC':0, 'GTA':0,  
  8  'GTG':0, 'TAT':0, 'TAC':0, 'TAA':0, 'TAG':0,  
  9  'CAT':0, 'CAC':0, 'CAA':0, 'CAG':0, 'AAT':0,  
 10  'AAC':0, 'AAA':0, 'AAG':0, 'GAT':0, 'GAC':0,  
 11  'GAA':0, 'GAG':0, 'TCT':0, 'TCC':0, 'TCA':0,  
 12  'TCG':0, 'CCT':0, 'CCC':0, 'CCA':0, 'CCG':0,  
 13  'ACT':0, 'ACC':0, 'ACA':0, 'ACG':0, 'GCT':0,  
 14  'GCC':0, 'GCA':0, 'GCG':0, 'TGT':0, 'TGC':0,  
 15  'TGA':0, 'TGG':0, 'CGT':0, 'CGC':0, 'CGA':0,  
 16  'CGG':0, 'AGT':0, 'AGC':0, 'AGA':0, 'AGG':0,  
 17  'GGT':0, 'GGC':0, 'GGA':0, 'GGG':0} 
 18   
 19   
 20  # this dictionary is used to know which codons encode the same AA. 
 21  SynonymousCodons = {'CYS': ['TGT', 'TGC'], 'ASP': ['GAT', 'GAC'], 
 22  'SER': ['TCT', 'TCG', 'TCA', 'TCC', 'AGC', 'AGT'], 
 23  'GLN': ['CAA', 'CAG'], 'MET': ['ATG'], 'ASN': ['AAC', 'AAT'], 
 24  'PRO': ['CCT', 'CCG', 'CCA', 'CCC'], 'LYS': ['AAG', 'AAA'], 
 25  'STOP': ['TAG', 'TGA', 'TAA'], 'THR': ['ACC', 'ACA', 'ACG', 'ACT'], 
 26  'PHE': ['TTT', 'TTC'], 'ALA': ['GCA', 'GCC', 'GCG', 'GCT'], 
 27  'GLY': ['GGT', 'GGG', 'GGA', 'GGC'], 'ILE': ['ATC', 'ATA', 'ATT'], 
 28  'LEU': ['TTA', 'TTG', 'CTC', 'CTT', 'CTG', 'CTA'], 'HIS': ['CAT', 'CAC'], 
 29  'ARG': ['CGA', 'CGC', 'CGG', 'CGT', 'AGG', 'AGA'], 'TRP': ['TGG'], 
 30  'VAL': ['GTA', 'GTC', 'GTG', 'GTT'], 'GLU': ['GAG', 'GAA'], 'TYR': ['TAT', 'TAC']} 
 31   
 32   
33 -class CodonAdaptationIndex(object):
34 """A codon adaptaion index (CAI) implementation. 35 36 This class implements the codon adaptaion index (CAI) described by Sharp and 37 Li (Nucleic Acids Res. 1987 Feb 11;15(3):1281-95). 38 39 methods: 40 41 set_cai_index(Index): 42 43 This method sets-up an index to be used when calculating CAI for a gene. 44 Just pass a dictionary similar to the SharpEcoliIndex in CodonUsageIndices 45 module. 46 47 generate_index(FastaFile): 48 49 This method takes a location of a FastaFile and generates an index. This 50 index can later be used to calculate CAI of a gene. 51 52 cai_for_gene(DNAsequence): 53 54 This method uses the Index (either the one you set or the one you generated) 55 and returns the CAI for the DNA sequence. 56 57 print_index(): 58 This method prints out the index you used. 59 60 NOTE - This implementation does not currently cope with alternative genetic 61 codes, only the synonymous codons in the standard table are considered. 62 """
63 - def __init__(self):
64 self.index = {} 65 self.codon_count={}
66 67 # use this method with predefined CAI index
68 - def set_cai_index(self, Index):
69 self.index = Index
70
71 - def generate_index(self, FastaFile):
72 """Generate a codon usage index from a FASTA file of CDS sequences. 73 74 This method takes a location of a Fasta file containing CDS sequences 75 (which must all have a whole number of codons) and generates a codon 76 usage index. This index can later be used to calculate CAI of a gene. 77 """ 78 # first make sure i am not overwriting an existing index: 79 if self.index != {} or self.codon_count!={}: 80 raise ValueError("an index has already been set or a codon count has been done. cannot overwrite either.") 81 # count codon occurances in the file. 82 self._count_codons(FastaFile) 83 84 # now to calculate the index we first need to sum the number of times 85 # synonymous codons were used all together. 86 for AA in SynonymousCodons: 87 Sum=0.0 88 RCSU=[] # RCSU values are equal to CodonCount/((1/num of synonymous codons) * sum of all synonymous codons) 89 90 for codon in SynonymousCodons[AA]: 91 Sum += self.codon_count[codon] 92 # calculate the RSCU value for each of the codons 93 for codon in SynonymousCodons[AA]: 94 RCSU.append(self.codon_count[codon]/((1.0/len(SynonymousCodons[AA]))*Sum)) 95 # now generate the index W=RCSUi/RCSUmax: 96 RCSUmax = max(RCSU) 97 for i in range(len(SynonymousCodons[AA])): 98 self.index[SynonymousCodons[AA][i]]= RCSU[i]/RCSUmax
99 100
101 - def cai_for_gene(self, DNAsequence):
102 """Calculate the CAI (float) for the provided DNA sequence (string). 103 104 This method uses the Index (either the one you set or the one you generated) 105 and returns the CAI for the DNA sequence. 106 """ 107 caiValue = 0 108 LengthForCai = 0 109 # if no index is set or generated, the default SharpEcoliIndex will be used. 110 if self.index=={}: 111 self.set_cai_index(SharpEcoliIndex) 112 113 if DNAsequence.islower(): 114 DNAsequence = DNAsequence.upper() 115 for i in range (0,len(DNAsequence),3): 116 codon = DNAsequence[i:i+3] 117 if codon in self.index: 118 if codon!='ATG' and codon!= 'TGG': #these two codons are always one, exclude them. 119 caiValue += math.log(self.index[codon]) 120 LengthForCai += 1 121 elif codon not in ['TGA','TAA', 'TAG']: # some indices you will use may not include stop codons. 122 raise TypeError("illegal codon in sequence: %s.\n%s" % (codon, self.index)) 123 return math.exp(caiValue*(1.0/(LengthForCai-1)))
124
125 - def _count_codons(self, FastaFile):
126 handle = open(FastaFile, 'r') 127 128 # make the codon dictionary local 129 self.codon_count = CodonsDict.copy() 130 131 # iterate over sequence and count all the codons in the FastaFile. 132 for cur_record in SeqIO.parse(handle, "fasta"): 133 # make sure the sequence is lower case 134 if str(cur_record.seq).islower(): 135 DNAsequence = str(cur_record.seq).upper() 136 else: 137 DNAsequence = str(cur_record.seq) 138 for i in range(0,len(DNAsequence),3): 139 codon = DNAsequence[i:i+3] 140 if codon in self.codon_count: 141 self.codon_count[codon] += 1 142 else: 143 raise TypeError("illegal codon %s in gene: %s" % (codon, cur_record.id)) 144 handle.close()
145 146 # this just gives the index when the objects is printed.
147 - def print_index (self):
148 """This method prints out the index you used.""" 149 for i in sorted(self.index): 150 print "%s\t%.3f" %(i, self.index[i])
151