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

Source Code for Module Bio.SeqUtils.CodonUsage

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