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