Package Bio :: Package CodonAlign :: Module CodonAlignment'
[hide private]
[frames] | no frames]

Source Code for Module Bio.CodonAlign.CodonAlignment'

  1  # Copyright 2013 by Zheng Ruan (zruan1991@gmail.com). 
  2  # All rights reserved. 
  3  # This code is part of the Biopython distribution and governed by its 
  4  # license.  Please see the LICENSE file that should have been included 
  5  # as part of this package. 
  6  """Code for dealing with Codon Alignment. 
  7   
  8  CodonAlignment class is interited from MultipleSeqAlignment class. This is 
  9  the core class to deal with codon alignment in biopython. 
 10   
 11  """ 
 12  from __future__ import division, print_function 
 13   
 14  __docformat__ = "epytext en"  # Don't just use plain text in epydoc API pages! 
 15   
 16  from Bio.Align import MultipleSeqAlignment 
 17  from Bio.SeqRecord import SeqRecord 
 18   
 19  from Bio.CodonAlign.CodonAlphabet import default_codon_table, default_codon_alphabet 
 20  from Bio.CodonAlign.CodonSeq import _get_codon_list, CodonSeq, cal_dn_ds 
 21  from Bio.CodonAlign.chisq import chisqprob 
22 23 -class CodonAlignment(MultipleSeqAlignment):
24 """Codon Alignment class that inherits from MultipleSeqAlignment. 25 26 >>> from Bio.Alphabet import generic_dna 27 >>> from Bio.SeqRecord import SeqRecord 28 >>> from Bio.Alphabet import IUPAC, Gapped 29 >>> a = SeqRecord(CodonSeq("AAAACGTCG", alphabet=default_codon_alphabet), id="Alpha") 30 >>> b = SeqRecord(CodonSeq("AAA---TCG", alphabet=default_codon_alphabet), id="Beta") 31 >>> c = SeqRecord(CodonSeq("AAAAGGTGG", alphabet=default_codon_alphabet), id="Gamma") 32 >>> print(CodonAlignment([a, b, c])) 33 CodonAlphabet(Standard) CodonAlignment with 3 rows and 9 columns (3 codons) 34 AAAACGTCG Alpha 35 AAA---TCG Beta 36 AAAAGGTGG Gamma 37 38 """
39 - def __init__(self, records='', name=None, alphabet=default_codon_alphabet):
40 41 MultipleSeqAlignment.__init__(self, records, alphabet=alphabet) 42 43 # check the type of the alignment to be nucleotide 44 for rec in self: 45 if not isinstance(rec.seq, CodonSeq): 46 raise TypeError("CodonSeq object are expected in each " 47 "SeqRecord in CodonAlignment") 48 49 assert self.get_alignment_length() % 3 == 0, \ 50 "Alignment length is not a triple number"
51
52 - def __str__(self):
53 """Return a multi-line string summary of the alignment. 54 55 This output is indicated to be readable, but large alignment 56 is shown truncated. A maximum of 20 rows (sequences) and 57 60 columns (20 codons) are shown, with the record identifiers. 58 This should fit nicely on a single screen. e.g. 59 60 """ 61 rows = len(self._records) 62 lines = ["%s CodonAlignment with %i rows and %i columns (%i codons)" 63 % (str(self._alphabet), rows, \ 64 self.get_alignment_length(), self.get_aln_length())] 65 66 if rows <= 60: 67 lines.extend([self._str_line(rec, length=60) \ 68 for rec in self._records]) 69 else: 70 lines.extend([self._str_line(rec, length=60) \ 71 for rec in self._records[:18]]) 72 lines.append("...") 73 lines.append(self._str_line(self._records[-1], length=60)) 74 return "\n".join(lines)
75
76 - def __getitem__(self, index, alphabet=None):
77 """Return a CodonAlignment object for single indexing 78 """ 79 if isinstance(index, int): 80 return self._records[index] 81 elif isinstance(index, slice): 82 return CodonAlignment(self._records[index], self._alphabet) 83 elif len(index) != 2: 84 raise TypeError("Invalid index type.") 85 # Handle double indexing 86 row_index, col_index = index 87 if isinstance(row_index, int): 88 return self._records[row_index][col_index] 89 elif isinstance(col_index, int): 90 return "".join(str(rec[col_index]) for rec in \ 91 self._records[row_index]) 92 else: 93 if alphabet is None: 94 from Bio.Alphabet import generic_nucleotide 95 return MultipleSeqAlignment((rec[col_index] for rec in \ 96 self._records[row_index]), 97 generic_nucleotide) 98 else: 99 return MultipleSeqAlignment((rec[col_index] for rec in \ 100 self._records[row_index]), 101 generic_nucleotide)
102
103 - def get_aln_length(self):
104 return self.get_alignment_length() // 3
105
106 - def toMultipleSeqAlignment(self):
107 """Return a MultipleSeqAlignment containing all the 108 SeqRecord in the CodonAlignment using Seq to store 109 sequences 110 """ 111 alignments = [SeqRecord(rec.seq.toSeq(), id=rec.id) for \ 112 rec in self._records] 113 return MultipleSeqAlignment(alignments)
114
115 - def get_dn_ds_matrix(self, method="NG86"):
116 """Available methods include NG86, LWL85, YN00 and ML. 117 """ 118 from Bio.Phylo.TreeConstruction import _DistanceMatrix as DM 119 names = [i.id for i in self._records] 120 size = len(self._records) 121 dn_matrix = [] 122 ds_matrix = [] 123 for i in range(size): 124 dn_matrix.append([]) 125 ds_matrix.append([]) 126 for j in range(i+1): 127 if i != j: 128 dn, ds = cal_dn_ds(self._records[i], self._records[j], 129 method=method) 130 dn_matrix[i].append(dn) 131 ds_matrix[i].append(ds) 132 else: 133 dn_matrix[i].append(0.0) 134 ds_matrix[i].append(0.0) 135 dn_dm = DM(names, matrix=dn_matrix) 136 ds_dm = DM(names, matrix=ds_matrix) 137 return dn_dm, ds_dm
138
139 - def get_dn_ds_tree(self, dn_ds_method="NG86", tree_method="UPGMA"):
140 """Method for constructing dn tree and ds tree. 141 Argument: 142 - dn_ds_method - Available methods include NG86, LWL85, YN00 143 and ML. 144 - tree_method - Available methods include UPGMA and NJ. 145 """ 146 from Bio.Phylo.TreeConstruction import DistanceTreeConstructor 147 dn_dm, ds_dm = self.get_dn_ds_matrix(method=dn_ds_method) 148 dn_constructor = DistanceTreeConstructor() 149 ds_constructor = DistanceTreeConstructor() 150 if tree_method == "UPGMA": 151 dn_tree = dn_constructor.upgma(dn_dm) 152 ds_tree = ds_constructor.upgma(ds_dm) 153 elif tree_method == "NJ": 154 dn_tree = dn_constructor.nj(dn_dm) 155 ds_tree = ds_constructor.nj(ds_dm) 156 else: 157 raise RuntimeError("Unkown tree method ({0}). Only NJ and UPGMA " 158 "are accepted.".format(tree_method)) 159 return dn_tree, ds_tree
160 161 @classmethod
162 - def from_msa(cls, align, alphabet=default_codon_alphabet):
163 """Function to convert a MultipleSeqAlignment to CodonAlignment. 164 It is the user's responsibility to ensure all the requirement 165 needed by CodonAlignment is met. 166 """ 167 rec = [SeqRecord(CodonSeq(str(i.seq), alphabet=alphabet), id=i.id) \ 168 for i in align._records] 169 return cls(rec, alphabet=alphabet)
170
171 172 -def mktest(codon_alns, codon_table=default_codon_table, alpha=0.05):
173 """McDonald-Kreitman test for neutrality (PMID: 1904993) This method 174 counts changes rather than sites (http://mkt.uab.es/mkt/help_mkt.asp). 175 Arguments: 176 - codon_alns - list of CodonAlignment to compare (each 177 CodonAlignment object corresponds to gene 178 sampled from a species) 179 180 Return the p-value of test result 181 """ 182 import copy 183 if not all([isinstance(i, CodonAlignment) for i in codon_alns]): 184 raise TypeError("mktest accept CodonAlignment list.") 185 codon_aln_len = [i.get_alignment_length() for i in codon_alns] 186 if len(set(codon_aln_len)) != 1: 187 raise RuntimeError("CodonAlignment object for mktest should be of" 188 " equal length.") 189 codon_num = codon_aln_len[0]//3 190 # prepare codon_dict (taking stop codon as an extra amino acid) 191 codon_dict = copy.deepcopy(codon_table.forward_table) 192 for stop in codon_table.stop_codons: 193 codon_dict[stop] = 'stop' 194 # prepare codon_lst 195 codon_lst = [] 196 for codon_aln in codon_alns: 197 codon_lst.append([]) 198 for i in codon_aln: 199 codon_lst[-1].append(_get_codon_list(i.seq)) 200 codon_set = [] 201 for i in range(codon_num): 202 uniq_codons = [] 203 for j in codon_lst: 204 uniq_codon = set([k[i] for k in j]) 205 uniq_codons.append(uniq_codon) 206 codon_set.append(uniq_codons) 207 syn_fix, nonsyn_fix, syn_poly, nonsyn_poly = 0, 0, 0, 0 208 G, nonsyn_G = _get_codon2codon_matrix(codon_table=codon_table) 209 for i in codon_set: 210 all_codon = i[0].union(*i[1:]) 211 if '-' in all_codon or len(all_codon) == 1: 212 continue 213 fix_or_not = all([len(k) == 1 for k in i]) 214 if fix_or_not: 215 # fixed 216 nonsyn_subgraph = _get_subgraph(all_codon, nonsyn_G) 217 subgraph = _get_subgraph(all_codon, G) 218 this_non = _count_replacement(all_codon, nonsyn_subgraph) 219 this_syn = _count_replacement(all_codon, subgraph) - this_non 220 nonsyn_fix += this_non 221 syn_fix += this_syn 222 else: 223 # not fixed 224 nonsyn_subgraph = _get_subgraph(all_codon, nonsyn_G) 225 subgraph = _get_subgraph(all_codon, G) 226 this_non = _count_replacement(all_codon, nonsyn_subgraph) 227 this_syn = _count_replacement(all_codon, subgraph) - this_non 228 nonsyn_poly += this_non 229 syn_poly += this_syn 230 return _G_test([syn_fix, nonsyn_fix, syn_poly, nonsyn_poly])
231
232 233 -def _get_codon2codon_matrix(codon_table=default_codon_table):
234 """Function to get codon codon subsitution matrix. Elements 235 in the matrix are number of synonymous and nonsynonymous 236 substitutions required for the substitution (PRIVATE). 237 """ 238 import platform 239 if platform.python_implementation() == 'PyPy': 240 import numpypy as np 241 else: 242 import numpy as np 243 base_tuple = ('A', 'T', 'C', 'G') 244 codons = [i for i in list(codon_table.forward_table.keys()) + \ 245 codon_table.stop_codons if 'U' not in i] 246 # set up codon_dict considering stop codons 247 codon_dict = codon_table.forward_table 248 for stop in codon_table.stop_codons: 249 codon_dict[stop] = 'stop' 250 # count site 251 num = len(codons) 252 G = {} # graph for substitution 253 nonsyn_G = {} # graph for nonsynonymous substitution 254 graph = {} 255 graph_nonsyn = {} 256 for i, codon in enumerate(codons): 257 graph[codon] = {} 258 graph_nonsyn[codon] = {} 259 for p, b in enumerate(codon): 260 for j in base_tuple: 261 tmp_codon = codon[0:p] + j + codon[p+1:] 262 if codon_dict[codon] != codon_dict[tmp_codon]: 263 graph_nonsyn[codon][tmp_codon] = 1 264 graph[codon][tmp_codon] = 1 265 else: 266 if codon != tmp_codon: 267 graph_nonsyn[codon][tmp_codon] = 0.1 268 graph[codon][tmp_codon] = 1 269 for codon1 in codons: 270 nonsyn_G[codon1] = {} 271 G[codon1] = {} 272 for codon2 in codons: 273 if codon1 == codon2: 274 nonsyn_G[codon1][codon2] = 0 275 G[codon1][codon2] = 0 276 else: 277 nonsyn_G[codon1][codon2] = _dijkstra(graph_nonsyn, codon1, 278 codon2) 279 G[codon1][codon2] = _dijkstra(graph, codon1, codon2) 280 return G, nonsyn_G
281
282 283 -def _dijkstra(graph, start, end):
284 """ 285 Dijkstra's algorithm Python implementation. 286 Algorithm adapted from 287 http://thomas.pelletier.im/2010/02/dijkstras-algorithm-python-implementation/. 288 However, an abvious bug in 289 if D[child_node] >(<) D[node] + child_value: 290 is fixed. 291 This function will return the distance between start and end. 292 293 Arguments: 294 graph: Dictionnary of dictionnary (keys are vertices). 295 start: Start vertex. 296 end: End vertex. 297 Output: 298 List of vertices from the beggining to the end. 299 """ 300 D = {} # Final distances dict 301 P = {} # Predecessor dict 302 # Fill the dicts with default values 303 for node in graph.keys(): 304 D[node] = 100 # Vertices are unreachable 305 P[node] = "" # Vertices have no predecessors 306 D[start] = 0 # The start vertex needs no move 307 unseen_nodes = list(graph.keys()) # All nodes are unseen 308 while len(unseen_nodes) > 0: 309 # Select the node with the lowest value in D (final distance) 310 shortest = None 311 node = '' 312 for temp_node in unseen_nodes: 313 if shortest == None: 314 shortest = D[temp_node] 315 node = temp_node 316 elif D[temp_node] < shortest: 317 shortest = D[temp_node] 318 node = temp_node 319 # Remove the selected node from unseen_nodes 320 unseen_nodes.remove(node) 321 # For each child (ie: connected vertex) of the current node 322 for child_node, child_value in graph[node].items(): 323 if D[child_node] > D[node] + child_value: 324 D[child_node] = D[node] + child_value 325 # To go to child_node, you have to go through node 326 P[child_node] = node 327 if node == end: break 328 # Set a clean path 329 path = [] 330 # We begin from the end 331 node = end 332 distance = 0 333 # While we are not arrived at the beginning 334 while not (node == start): 335 if path.count(node) == 0: 336 path.insert(0, node) # Insert the predecessor of the current node 337 node = P[node] # The current node becomes its predecessor 338 else: 339 break 340 path.insert(0, start) # Finally, insert the start vertex 341 for i in range(len(path)-1): 342 distance += graph[path[i]][path[i+1]] 343 return distance
344
345 346 -def _count_replacement(codon_set, G):
347 """Count replacement needed for a given codon_set (PRIVATE). 348 """ 349 from math import floor 350 if len(codon_set) == 1: 351 return 0, 0 352 elif len(codon_set) == 2: 353 codons = list(codon_set) 354 return floor(G[codons[0]][codons[1]]) 355 else: 356 codons = list(codon_set) 357 return _prim(G)
358
359 360 -def _prim(G):
361 """Prim's algorithm to find minimum spanning tree. Code is adapted 362 from 363 http://programmingpraxis.com/2010/04/09/minimum-spanning-tree-prims-algorithm/ 364 (PRIVATE). 365 """ 366 from math import floor 367 from collections import defaultdict 368 from heapq import heapify, heappop, heappush 369 nodes = [] 370 edges = [] 371 for i in G.keys(): 372 nodes.append(i) 373 for j in G[i]: 374 if (i, j, G[i][j]) not in edges and (j, i, G[i][j]) not in edges: 375 edges.append((i, j, G[i][j])) 376 conn = defaultdict(list) 377 for n1, n2, c in edges: 378 conn[n1].append((c, n1, n2)) 379 conn[n2].append((c, n2, n1)) 380 mst = [] # minimum spanning tree 381 used = set(nodes[0]) 382 usable_edges = conn[nodes[0]][:] 383 heapify(usable_edges) 384 while usable_edges: 385 cost, n1, n2 = heappop(usable_edges) 386 if n2 not in used: 387 used.add(n2) 388 mst.append((n1, n2, cost)) 389 for e in conn[n2]: 390 if e[2] not in used: 391 heappush(usable_edges, e) 392 length = 0 393 for p in mst: 394 length += floor(p[2]) 395 return length
396
397 398 -def _get_subgraph(codons, G):
399 """Get the subgraph that contains all codons in list (PRIVATE). 400 """ 401 subgraph = {} 402 for i in codons: 403 subgraph[i] = {} 404 for j in codons: 405 if i != j: subgraph[i][j] = G[i][j] 406 return subgraph
407
408 409 -def _G_test(site_counts):
410 """G test for 2x2 contingency table (PRIVATE). 411 Argument: 412 - site_counts - [syn_fix, nonsyn_fix, syn_poly, nonsyn_poly] 413 414 >>> round(_G_test([17, 7, 42, 2]), 7) 415 0.004924 416 """ 417 # TODO: 418 # Apply continuity correction for Chi-square test. 419 from math import log 420 #from scipy.stats import chi2 421 G = 0 422 tot = sum(site_counts) 423 tot_syn = site_counts[0] + site_counts[2] 424 tot_non = site_counts[1] + site_counts[3] 425 tot_fix = sum(site_counts[:2]) 426 tot_poly = sum(site_counts[2:]) 427 exp = [tot_fix*tot_syn/tot, tot_fix*tot_non/tot, 428 tot_poly*tot_syn/tot, tot_poly*tot_non/tot] 429 for obs, ex in zip(site_counts, exp): 430 G += obs*log(obs/ex) 431 G *= 2 432 #return 1-chi2.cdf(G, 1) # only 1 dof for 2x2 table 433 return chisqprob(G, 1)
434 435 436 if __name__ == "__main__": 437 from Bio._utils import run_doctest 438 run_doctest() 439