Bio.Align.AlignInfo

29 """Calculate summary info about the alignment. 30 31 This class should be used to caclculate information summarizing the 32 results of an alignment. This may either be straight consensus info 33 or more complicated things. 34 """

35 - def __init__(self, alignment):

36 """Initialize with the alignment to calculate information on. 37 ic_vector attribute. A dictionary. Keys: column numbers. Values: 38 """ 39 self.alignment = alignment 40 self.ic_vector = {}

41

42 - def dumb_consensus(self, threshold = .7, ambiguous = "X", 43 consensus_alpha = None, require_multiple = 0):

44 """Output a fast consensus sequence of the alignment. 45 46 This doesn't do anything fancy at all. It will just go through the 47 sequence residue by residue and count up the number of each type 48 of residue (ie. A or G or T or C for DNA) in all sequences in the 49 alignment. If the percentage of the most common residue type is 50 greater then the passed threshold, then we will add that residue type, 51 otherwise an ambiguous character will be added. 52 53 This could be made a lot fancier (ie. to take a substitution matrix 54 into account), but it just meant for a quick and dirty consensus. 55 56 Arguments: 57 o threshold - The threshold value that is required to add a particular 58 atom. 59 o ambiguous - The ambiguous character to be added when the threshold is 60 not reached. 61 o consensus_alpha - The alphabet to return for the consensus sequence. 62 If this is None, then we will try to guess the alphabet. 63 o require_multiple - If set as 1, this will require that more than 64 1 sequence be part of an alignment to put it in the consensus (ie. 65 not just 1 sequence and gaps). 66 """ 67 # Iddo Friedberg, 1-JUL-2004: changed ambiguous default to "X" 68 consensus = '' 69 70 # find the length of the consensus we are creating 71 con_len = self.alignment.get_alignment_length() 72 73 # go through each seq item 74 for n in range(con_len): 75 # keep track of the counts of the different atoms we get 76 atom_dict = {} 77 num_atoms = 0 78 79 for record in self.alignment._records: 80 # make sure we haven't run past the end of any sequences 81 # if they are of different lengths 82 if n < len(record.seq): 83 if record.seq[n] != '-' and record.seq[n] != '.': 84 if record.seq[n] not in atom_dict: 85 atom_dict[record.seq[n]] = 1 86 else: 87 atom_dict[record.seq[n]] += 1 88 89 num_atoms = num_atoms + 1 90 91 max_atoms = [] 92 max_size = 0 93 94 for atom in atom_dict: 95 if atom_dict[atom] > max_size: 96 max_atoms = [atom] 97 max_size = atom_dict[atom] 98 elif atom_dict[atom] == max_size: 99 max_atoms.append(atom) 100 101 if require_multiple and num_atoms == 1: 102 consensus += ambiguous 103 elif (len(max_atoms) == 1) and ((float(max_size)/float(num_atoms)) 104 >= threshold): 105 consensus += max_atoms[0] 106 else: 107 consensus += ambiguous 108 109 # we need to guess a consensus alphabet if one isn't specified 110 if consensus_alpha is None: 111 consensus_alpha = self._guess_consensus_alphabet(ambiguous) 112 113 return Seq(consensus, consensus_alpha)

114

115 - def gap_consensus(self, threshold = .7, ambiguous = "X", 116 consensus_alpha = None, require_multiple = 0):

117 """Same as dumb_consensus(), but allows gap on the output. 118 119 Things to do: Let the user define that with only one gap, the result 120 character in consensus is gap. Let the user select gap character, now 121 it takes the same is input. 122 """ 123 # Iddo Friedberg, 1-JUL-2004: changed ambiguous default to "X" 124 consensus = '' 125 126 # find the length of the consensus we are creating 127 con_len = self.alignment.get_alignment_length() 128 129 # go through each seq item 130 for n in range(con_len): 131 # keep track of the counts of the different atoms we get 132 atom_dict = {} 133 num_atoms = 0 134 135 for record in self.alignment._records: 136 # make sure we haven't run past the end of any sequences 137 # if they are of different lengths 138 if n < len(record.seq): 139 if record.seq[n] not in atom_dict: 140 atom_dict[record.seq[n]] = 1 141 else: 142 atom_dict[record.seq[n]] += 1 143 144 num_atoms += 1 145 146 max_atoms = [] 147 max_size = 0 148 149 for atom in atom_dict: 150 if atom_dict[atom] > max_size: 151 max_atoms = [atom] 152 max_size = atom_dict[atom] 153 elif atom_dict[atom] == max_size: 154 max_atoms.append(atom) 155 156 if require_multiple and num_atoms == 1: 157 consensus += ambiguous 158 elif (len(max_atoms) == 1) and ((float(max_size)/float(num_atoms)) 159 >= threshold): 160 consensus += max_atoms[0] 161 else: 162 consensus += ambiguous 163 164 # we need to guess a consensus alphabet if one isn't specified 165 if consensus_alpha is None: 166 #TODO - Should we make this into a Gapped alphabet? 167 consensus_alpha = self._guess_consensus_alphabet(ambiguous) 168 169 return Seq(consensus, consensus_alpha)

170

171 - def _guess_consensus_alphabet(self, ambiguous):

172 """Pick an (ungapped) alphabet for an alignment consesus sequence. 173 174 This just looks at the sequences we have, checks their type, and 175 returns as appropriate type which seems to make sense with the 176 sequences we've got. 177 """ 178 #Start with the (un-gapped version of) the alignment alphabet 179 a = Alphabet._get_base_alphabet(self.alignment._alphabet) 180 181 #Now check its compatible with all the rest of the sequences 182 for record in self.alignment: 183 #Get the (un-gapped version of) the sequence's alphabet 184 alt = Alphabet._get_base_alphabet(record.seq.alphabet) 185 if not isinstance(alt, a.__class__): 186 raise ValueError("Alignment contains a sequence with \ 187 an incompatible alphabet.") 188 189 #Check the ambiguous character we are going to use in the consensus 190 #is in the alphabet's list of valid letters (if defined). 191 if hasattr(a, "letters") and a.letters is not None \ 192 and ambiguous not in a.letters: 193 #We'll need to pick a more generic alphabet... 194 if isinstance(a, IUPAC.IUPACUnambiguousDNA): 195 if ambiguous in IUPAC.IUPACUnambiguousDNA().letters: 196 a = IUPAC.IUPACUnambiguousDNA() 197 else: 198 a = Alphabet.generic_dna 199 elif isinstance(a, IUPAC.IUPACUnambiguousRNA): 200 if ambiguous in IUPAC.IUPACUnambiguousRNA().letters: 201 a = IUPAC.IUPACUnambiguousRNA() 202 else: 203 a = Alphabet.generic_rna 204 elif isinstance(a, IUPAC.IUPACProtein): 205 if ambiguous in IUPAC.ExtendedIUPACProtein().letters: 206 a = IUPAC.ExtendedIUPACProtein() 207 else: 208 a = Alphabet.generic_protein 209 else: 210 a = Alphabet.single_letter_alphabet 211 return a

212

213 - def replacement_dictionary(self, skip_chars = []):

214 """Generate a replacement dictionary to plug into a substitution matrix 215 216 This should look at an alignment, and be able to generate the number 217 of substitutions of different residues for each other in the 218 aligned object. 219 220 Will then return a dictionary with this information: 221 {('A', 'C') : 10, ('C', 'A') : 12, ('G', 'C') : 15 ....} 222 223 This also treats weighted sequences. The following example shows how 224 we calculate the replacement dictionary. Given the following 225 multiple sequence alignments: 226 227 GTATC 0.5 228 AT--C 0.8 229 CTGTC 1.0 230 231 For the first column we have: 232 ('A', 'G') : 0.5 * 0.8 = 0.4 233 ('C', 'G') : 0.5 * 1.0 = 0.5 234 ('A', 'C') : 0.8 * 1.0 = 0.8 235 236 We then continue this for all of the columns in the alignment, summing 237 the information for each substitution in each column, until we end 238 up with the replacement dictionary. 239 240 Arguments: 241 o skip_chars - A list of characters to skip when creating the dictionary. 242 For instance, you might have Xs (screened stuff) or Ns, and not want 243 to include the ambiguity characters in the dictionary. 244 """ 245 # get a starting dictionary based on the alphabet of the alignment 246 rep_dict, skip_items = self._get_base_replacements(skip_chars) 247 248 # iterate through each record 249 for rec_num1 in range(len(self.alignment._records)): 250 # iterate through each record from one beyond the current record 251 # to the end of the list of records 252 for rec_num2 in range(rec_num1 + 1, len(self.alignment._records)): 253 # for each pair of records, compare the sequences and add 254 # the pertinent info to the dictionary 255 rep_dict = self._pair_replacement( 256 self.alignment._records[rec_num1].seq, 257 self.alignment._records[rec_num2].seq, 258 self.alignment._records[rec_num1].annotations.get('weight',1.0), 259 self.alignment._records[rec_num2].annotations.get('weight',1.0), 260 rep_dict, skip_items) 261 262 return rep_dict

263

264 - def _pair_replacement(self, seq1, seq2, weight1, weight2, 265 start_dict, ignore_chars):

266 """Compare two sequences and generate info on the replacements seen. 267 268 Arguments: 269 o seq1, seq2 - The two sequences to compare. 270 o weight1, weight2 - The relative weights of seq1 and seq2. 271 o start_dict - The dictionary containing the starting replacement 272 info that we will modify. 273 o ignore_chars - A list of characters to ignore when calculating 274 replacements (ie. '-'). 275 276 Returns: 277 o A replacment dictionary which is modified from initial_dict with 278 the information from the sequence comparison. 279 """ 280 # loop through each residue in the sequences 281 for residue_num in range(len(seq1)): 282 residue1 = seq1[residue_num] 283 try: 284 residue2 = seq2[residue_num] 285 # if seq2 is shorter, then we just stop looking at replacements 286 # and return the information 287 except IndexError: 288 return start_dict 289 290 # if the two residues are characters we want to count 291 if (residue1 not in ignore_chars) and (residue2 not in ignore_chars): 292 try: 293 # add info about the replacement to the dictionary, 294 # modified by the sequence weights 295 start_dict[(residue1, residue2)] += weight1 * weight2 296 297 # if we get a key error, then we've got a problem with alphabets 298 except KeyError: 299 raise ValueError("Residues %s, %s not found in alphabet %s" 300 % (residue1, residue2, 301 self.alignment._alphabet)) 302 303 return start_dict

304

305 - def _get_all_letters(self):

306 """Returns a string containing the expected letters in the alignment.""" 307 all_letters = self.alignment._alphabet.letters 308 if all_letters is None \ 309 or (isinstance(self.alignment._alphabet, Alphabet.Gapped) 310 and all_letters == self.alignment._alphabet.gap_char): 311 #We are dealing with a generic alphabet class where the 312 #letters are not defined! We must build a list of the 313 #letters used... 314 set_letters = set() 315 for record in self.alignment: 316 #Note the built in set does not have a union_update 317 #which was provided by the sets module's Set 318 set_letters = set_letters.union(record.seq) 319 list_letters = list(set_letters) 320 list_letters.sort() 321 all_letters = "".join(list_letters) 322 return all_letters

323

324 - def _get_base_replacements(self, skip_items = []):

325 """Get a zeroed dictionary of all possible letter combinations. 326 327 This looks at the type of alphabet and gets the letters for it. 328 It then creates a dictionary with all possible combinations of these 329 letters as keys (ie. ('A', 'G')) and sets the values as zero. 330 331 Returns: 332 o The base dictionary created 333 o A list of alphabet items to skip when filling the dictionary.Right 334 now the only thing I can imagine in this list is gap characters, but 335 maybe X's or something else might be useful later. This will also 336 include any characters that are specified to be skipped. 337 """ 338 base_dictionary = {} 339 all_letters = self._get_all_letters() 340 341 # if we have a gapped alphabet we need to find the gap character 342 # and drop it out 343 if isinstance(self.alignment._alphabet, Alphabet.Gapped): 344 skip_items.append(self.alignment._alphabet.gap_char) 345 all_letters = all_letters.replace(self.alignment._alphabet.gap_char,'') 346 347 # now create the dictionary 348 for first_letter in all_letters: 349 for second_letter in all_letters: 350 if first_letter not in skip_items and \ 351 second_letter not in skip_items: 352 base_dictionary[(first_letter, second_letter)] = 0 353 354 return base_dictionary, skip_items

355

356 - def pos_specific_score_matrix(self, axis_seq = None, 357 chars_to_ignore = []):

358 """Create a position specific score matrix object for the alignment. 359 360 This creates a position specific score matrix (pssm) which is an 361 alternative method to look at a consensus sequence. 362 363 Arguments: 364 o chars_to_ignore - A listing of all characters not to include in 365 the pssm. If the alignment alphabet declares a gap character, 366 then it will be excluded automatically. 367 o axis_seq - An optional argument specifying the sequence to 368 put on the axis of the PSSM. This should be a Seq object. If nothing 369 is specified, the consensus sequence, calculated with default 370 parameters, will be used. 371 372 Returns: 373 o A PSSM (position specific score matrix) object. 374 """ 375 # determine all of the letters we have to deal with 376 all_letters = self._get_all_letters() 377 assert all_letters 378 379 if not isinstance(chars_to_ignore, list): 380 raise TypeError("chars_to_ignore should be a list.") 381 382 # if we have a gap char, add it to stuff to ignore 383 if isinstance(self.alignment._alphabet, Alphabet.Gapped): 384 chars_to_ignore.append(self.alignment._alphabet.gap_char) 385 386 for char in chars_to_ignore: 387 all_letters = all_letters.replace(char, '') 388 389 if axis_seq: 390 left_seq = axis_seq 391 assert len(axis_seq) == self.alignment.get_alignment_length() 392 else: 393 left_seq = self.dumb_consensus() 394 395 pssm_info = [] 396 # now start looping through all of the sequences and getting info 397 for residue_num in range(len(left_seq)): 398 score_dict = self._get_base_letters(all_letters) 399 for record in self.alignment._records: 400 try: 401 this_residue = record.seq[residue_num] 402 # if we hit an index error we've run out of sequence and 403 # should not add new residues 404 except IndexError: 405 this_residue = None 406 407 if this_residue and this_residue not in chars_to_ignore: 408 weight = record.annotations.get('weight', 1.0) 409 try: 410 score_dict[this_residue] += weight 411 # if we get a KeyError then we have an alphabet problem 412 except KeyError: 413 raise ValueError("Residue %s not found in alphabet %s" 414 % (this_residue, 415 self.alignment._alphabet)) 416 417 pssm_info.append((left_seq[residue_num], 418 score_dict)) 419 420 return PSSM(pssm_info)

421

422 - def _get_base_letters(self, letters):

423 """Create a zeroed dictionary with all of the specified letters. 424 """ 425 base_info = {} 426 for letter in letters: 427 base_info[letter] = 0 428 429 return base_info

430

431 - def information_content(self, start = 0, 432 end = None, 433 e_freq_table = None, log_base = 2, 434 chars_to_ignore = []):

435 """Calculate the information content for each residue along an alignment. 436 437 Arguments: 438 o start, end - The starting an ending points to calculate the 439 information content. These points should be relative to the first 440 sequence in the alignment, starting at zero (ie. even if the 'real' 441 first position in the seq is 203 in the initial sequence, for 442 the info content, we need to use zero). This defaults to the entire 443 length of the first sequence. 444 o e_freq_table - A FreqTable object specifying the expected frequencies 445 for each letter in the alphabet we are using (e.g. {'G' : 0.4, 446 'C' : 0.4, 'T' : 0.1, 'A' : 0.1}). Gap characters should not be 447 included, since these should not have expected frequencies. 448 o log_base - The base of the logathrim to use in calculating the 449 information content. This defaults to 2 so the info is in bits. 450 o chars_to_ignore - A listing of characterw which should be ignored 451 in calculating the info content. 452 453 Returns: 454 o A number representing the info content for the specified region. 455 456 Please see the Biopython manual for more information on how information 457 content is calculated. 458 """ 459 # if no end was specified, then we default to the end of the sequence 460 if end is None: 461 end = len(self.alignment._records[0].seq) 462 463 if start < 0 or end > len(self.alignment._records[0].seq): 464 raise ValueError("Start (%s) and end (%s) are not in the \ 465 range %s to %s" 466 % (start, end, 0, len(self.alignment._records[0].seq))) 467 # determine random expected frequencies, if necessary 468 random_expected = None 469 if not e_freq_table: 470 #TODO - What about ambiguous alphabets? 471 base_alpha = Alphabet._get_base_alphabet(self.alignment._alphabet) 472 if isinstance(base_alpha, Alphabet.ProteinAlphabet): 473 random_expected = Protein20Random 474 elif isinstance(base_alpha, Alphabet.NucleotideAlphabet): 475 random_expected = Nucleotide4Random 476 else: 477 errstr = "Error in alphabet: not Nucleotide or Protein, " 478 errstr += "supply expected frequencies" 479 raise ValueError(errstr) 480 del base_alpha 481 elif not isinstance(e_freq_table, FreqTable.FreqTable): 482 raise ValueError("e_freq_table should be a FreqTable object") 483 484 # determine all of the letters we have to deal with 485 all_letters = self._get_all_letters() 486 for char in chars_to_ignore: 487 all_letters = all_letters.replace(char, '') 488 489 info_content = {} 490 for residue_num in range(start, end): 491 freq_dict = self._get_letter_freqs(residue_num, 492 self.alignment._records, 493 all_letters, chars_to_ignore) 494 # print freq_dict, 495 column_score = self._get_column_info_content(freq_dict, 496 e_freq_table, 497 log_base, 498 random_expected) 499 500 info_content[residue_num] = column_score 501 # sum up the score 502 total_info = sum(info_content.itervalues()) 503 # fill in the ic_vector member: holds IC for each column 504 for i in info_content: 505 self.ic_vector[i] = info_content[i] 506 return total_info

507

508 - def _get_letter_freqs(self, residue_num, all_records, letters, to_ignore):

509 """Determine the frequency of specific letters in the alignment. 510 511 Arguments: 512 o residue_num - The number of the column we are getting frequencies 513 from. 514 o all_records - All of the SeqRecords in the alignment. 515 o letters - The letters we are interested in getting the frequency 516 for. 517 o to_ignore - Letters we are specifically supposed to ignore. 518 519 This will calculate the frequencies of each of the specified letters 520 in the alignment at the given frequency, and return this as a 521 dictionary where the keys are the letters and the values are the 522 frequencies. 523 """ 524 freq_info = self._get_base_letters(letters) 525 526 total_count = 0 527 # collect the count info into the dictionary for all the records 528 for record in all_records: 529 try: 530 if record.seq[residue_num] not in to_ignore: 531 weight = record.annotations.get('weight',1.0) 532 freq_info[record.seq[residue_num]] += weight 533 total_count += weight 534 # getting a key error means we've got a problem with the alphabet 535 except KeyError: 536 raise ValueError("Residue %s not found in alphabet %s" 537 % (record.seq[residue_num], 538 self.alignment._alphabet)) 539 540 if total_count == 0: 541 # This column must be entirely ignored characters 542 for letter in freq_info: 543 assert freq_info[letter] == 0 544 #TODO - Map this to NA or NaN? 545 else: 546 # now convert the counts into frequencies 547 for letter in freq_info: 548 freq_info[letter] = freq_info[letter] / total_count 549 550 return freq_info

551

552 - def _get_column_info_content(self, obs_freq, e_freq_table, log_base, 553 random_expected):

554 """Calculate the information content for a column. 555 556 Arguments: 557 o obs_freq - The frequencies observed for each letter in the column. 558 o e_freq_table - An optional argument specifying the expected 559 frequencies for each letter. This is a SubsMat.FreqTable instance. 560 o log_base - The base of the logathrim to use in calculating the 561 info content. 562 """ 563 try: 564 gap_char = self.alignment._alphabet.gap_char 565 except AttributeError: 566 #The alphabet doesn't declare a gap - there could be none 567 #in the sequence... or just a vague alphabet. 568 gap_char = "-" # Safe? 569 570 if e_freq_table: 571 if not isinstance(e_freq_table, FreqTable.FreqTable): 572 raise ValueError("e_freq_table should be a FreqTable object") 573 # check the expected freq information to make sure it is good 574 for key in obs_freq: 575 if (key != gap_char and key not in e_freq_table): 576 raise ValueError("Expected frequency letters %s " 577 "do not match observed %s" 578 % (e_freq_table.keys(), 579 obs_freq.keys() - [gap_char])) 580 581 total_info = 0.0 582 583 for letter in obs_freq: 584 inner_log = 0.0 585 # if we have expected frequencies, modify the log value by them 586 # gap characters do not have expected frequencies, so they 587 # should just be the observed frequency. 588 if letter != gap_char: 589 if e_freq_table: 590 inner_log = obs_freq[letter] / e_freq_table[letter] 591 else: 592 inner_log = obs_freq[letter] / random_expected 593 # if the observed frequency is zero, we don't add any info to the 594 # total information content 595 if inner_log > 0: 596 letter_info = (obs_freq[letter] * 597 math.log(inner_log) / math.log(log_base)) 598 total_info += letter_info 599 return total_info

600

601 - def get_column(self,col):

602 return self.alignment.get_column(col)

Source Code for Module Bio.Align.AlignInfo