Package Bio :: Module pairwise2
[hide private]
[frames] | no frames]

Source Code for Module Bio.pairwise2

  1  # Copyright 2002 by Jeffrey Chang.  All rights reserved. 
  2  # This code is part of the Biopython distribution and governed by its 
  3  # license.  Please see the LICENSE file that should have been included 
  4  # as part of this package. 
  5   
  6  """This package implements pairwise sequence alignment using a dynamic 
  7  programming algorithm. 
  8   
  9  This provides functions to get global and local alignments between two 
 10  sequences.  A global alignment finds the best concordance between all 
 11  characters in two sequences.  A local alignment finds just the 
 12  subsequences that align the best. 
 13   
 14  When doing alignments, you can specify the match score and gap 
 15  penalties.  The match score indicates the compatibility between an 
 16  alignment of two characters in the sequences.  Highly compatible 
 17  characters should be given positive scores, and incompatible ones 
 18  should be given negative scores or 0.  The gap penalties should be 
 19  negative. 
 20   
 21  The names of the alignment functions in this module follow the 
 22  convention 
 23  <alignment type>XX 
 24  where <alignment type> is either "global" or "local" and XX is a 2 
 25  character code indicating the parameters it takes.  The first 
 26  character indicates the parameters for matches (and mismatches), and 
 27  the second indicates the parameters for gap penalties. 
 28   
 29  The match parameters are:: 
 30   
 31      CODE  DESCRIPTION 
 32      x     No parameters.  Identical characters have score of 1, otherwise 0. 
 33      m     A match score is the score of identical chars, otherwise mismatch score. 
 34      d     A dictionary returns the score of any pair of characters. 
 35      c     A callback function returns scores. 
 36   
 37  The gap penalty parameters are:: 
 38   
 39      CODE  DESCRIPTION 
 40      x     No gap penalties. 
 41      s     Same open and extend gap penalties for both sequences. 
 42      d     The sequences have different open and extend gap penalties. 
 43      c     A callback function returns the gap penalties. 
 44   
 45  All the different alignment functions are contained in an object 
 46  "align".  For example: 
 47   
 48      >>> from Bio import pairwise2 
 49      >>> alignments = pairwise2.align.globalxx("ACCGT", "ACG") 
 50   
 51  will return a list of the alignments between the two strings.  The 
 52  parameters of the alignment function depends on the function called. 
 53  Some examples:: 
 54   
 55      # Find the best global alignment between the two sequences. 
 56      # Identical characters are given 1 point.  No points are deducted 
 57      # for mismatches or gaps. 
 58      >>> from Bio.pairwise2 import format_alignment 
 59      >>> for a in pairwise2.align.globalxx("ACCGT", "ACG"): 
 60      ...     print(format_alignment(*a)) 
 61      ACCGT 
 62      ||||| 
 63      AC-G- 
 64        Score=3 
 65      <BLANKLINE> 
 66      ACCGT 
 67      ||||| 
 68      A-CG- 
 69        Score=3 
 70      <BLANKLINE> 
 71   
 72      # Same thing as before, but with a local alignment. 
 73      >>> for a in pairwise2.align.localxx("ACCGT", "ACG"): 
 74      ...     print(format_alignment(*a)) 
 75      ACCGT 
 76      |||| 
 77      AC-G- 
 78        Score=3 
 79      <BLANKLINE> 
 80      ACCGT 
 81      |||| 
 82      A-CG- 
 83        Score=3 
 84      <BLANKLINE> 
 85   
 86      # Do a global alignment.  Identical characters are given 2 points, 
 87      # 1 point is deducted for each non-identical character. 
 88      >>> for a in pairwise2.align.globalmx("ACCGT", "ACG", 2, -1): 
 89      ...     print(format_alignment(*a)) 
 90      ACCGT 
 91      ||||| 
 92      AC-G- 
 93        Score=6 
 94      <BLANKLINE> 
 95      ACCGT 
 96      ||||| 
 97      A-CG- 
 98        Score=6 
 99      <BLANKLINE> 
100   
101      # Same as above, except now 0.5 points are deducted when opening a 
102      # gap, and 0.1 points are deducted when extending it. 
103      >>> for a in pairwise2.align.globalms("ACCGT", "ACG", 2, -1, -.5, -.1): 
104      ...     print(format_alignment(*a)) 
105      ACCGT 
106      ||||| 
107      AC-G- 
108        Score=5 
109      <BLANKLINE> 
110      ACCGT 
111      ||||| 
112      A-CG- 
113        Score=5 
114      <BLANKLINE> 
115   
116  The alignment function can also use known matrices already included in 
117  Biopython ( Bio.SubsMat -> MatrixInfo ):: 
118   
119      >>> from Bio.SubsMat import MatrixInfo as matlist 
120      >>> matrix = matlist.blosum62 
121      >>> for a in pairwise2.align.globaldx("KEVLA", "EVL", matrix): 
122      ...     print(format_alignment(*a)) 
123      KEVLA 
124      ||||| 
125      -EVL- 
126        Score=13 
127      <BLANKLINE> 
128   
129  To see a description of the parameters for a function, please look at 
130  the docstring for the function via the help function, e.g. 
131  type help(pairwise2.align.localds) at the Python prompt. 
132  """ 
133  # The alignment functions take some undocumented keyword parameters: 
134  # - penalize_extend_when_opening: boolean 
135  #   Whether to count an extension penalty when opening a gap.  If 
136  #   false, a gap of 1 is only penalize an "open" penalty, otherwise it 
137  #   is penalized "open+extend". 
138  # - penalize_end_gaps: boolean 
139  #   Whether to count the gaps at the ends of an alignment.  By 
140  #   default, they are counted for global alignments but not for local 
141  #   ones. Setting penalize_end_gaps to (boolean, boolean) allows you to 
142  #   specify for the two sequences separately whether gaps at the end of 
143  #   the alignment should be counted. 
144  # - gap_char: string 
145  #   Which character to use as a gap character in the alignment 
146  #   returned.  By default, uses '-'. 
147  # - force_generic: boolean 
148  #   Always use the generic, non-cached, dynamic programming function. 
149  #   For debugging. 
150  # - score_only: boolean 
151  #   Only get the best score, don't recover any alignments.  The return 
152  #   value of the function is the score. 
153  # - one_alignment_only: boolean 
154  #   Only recover one alignment. 
155   
156  from __future__ import print_function 
157   
158  __docformat__ = "restructuredtext en" 
159   
160  MAX_ALIGNMENTS = 1000   # maximum alignments recovered in traceback 
161   
162   
163 -class align(object):
164 """This class provides functions that do alignments.""" 165
166 - class alignment_function:
167 """This class is callable impersonates an alignment function. 168 The constructor takes the name of the function. This class 169 will decode the name of the function to figure out how to 170 interpret the parameters. 171 172 """ 173 # match code -> tuple of (parameters, docstring) 174 match2args = { 175 'x': ([], ''), 176 'm': (['match', 'mismatch'], 177 """match is the score to given to identical characters. mismatch is 178 the score given to non-identical ones."""), 179 'd': (['match_dict'], 180 """match_dict is a dictionary where the keys are tuples of pairs of 181 characters and the values are the scores, e.g. ("A", "C") : 2.5."""), 182 'c': (['match_fn'], 183 """match_fn is a callback function that takes two characters and 184 returns the score between them."""), 185 } 186 # penalty code -> tuple of (parameters, docstring) 187 penalty2args = { 188 'x': ([], ''), 189 's': (['open', 'extend'], 190 """open and extend are the gap penalties when a gap is opened and 191 extended. They should be negative."""), 192 'd': (['openA', 'extendA', 'openB', 'extendB'], 193 """openA and extendA are the gap penalties for sequenceA, and openB 194 and extendB for sequeneB. The penalties should be negative."""), 195 'c': (['gap_A_fn', 'gap_B_fn'], 196 """gap_A_fn and gap_B_fn are callback functions that takes 1) the 197 index where the gap is opened, and 2) the length of the gap. They 198 should return a gap penalty."""), 199 } 200
201 - def __init__(self, name):
202 # Check to make sure the name of the function is 203 # reasonable. 204 if name.startswith("global"): 205 if len(name) != 8: 206 raise AttributeError("function should be globalXX") 207 elif name.startswith("local"): 208 if len(name) != 7: 209 raise AttributeError("function should be localXX") 210 else: 211 raise AttributeError(name) 212 align_type, match_type, penalty_type = \ 213 name[:-2], name[-2], name[-1] 214 try: 215 match_args, match_doc = self.match2args[match_type] 216 except KeyError as x: 217 raise AttributeError("unknown match type %r" % match_type) 218 try: 219 penalty_args, penalty_doc = self.penalty2args[penalty_type] 220 except KeyError as x: 221 raise AttributeError("unknown penalty type %r" % penalty_type) 222 223 # Now get the names of the parameters to this function. 224 param_names = ['sequenceA', 'sequenceB'] 225 param_names.extend(match_args) 226 param_names.extend(penalty_args) 227 self.function_name = name 228 self.align_type = align_type 229 self.param_names = param_names 230 231 self.__name__ = self.function_name 232 # Set the doc string. 233 doc = "%s(%s) -> alignments\n" % ( 234 self.__name__, ', '.join(self.param_names)) 235 if match_doc: 236 doc += "\n%s\n" % match_doc 237 if penalty_doc: 238 doc += "\n%s\n" % penalty_doc 239 doc += ( 240 """\nalignments is a list of tuples (seqA, seqB, score, begin, end). 241 seqA and seqB are strings showing the alignment between the 242 sequences. score is the score of the alignment. begin and end 243 are indexes into seqA and seqB that indicate the where the 244 alignment occurs. 245 """) 246 self.__doc__ = doc
247
248 - def decode(self, *args, **keywds):
249 # Decode the arguments for the _align function. keywds 250 # will get passed to it, so translate the arguments to 251 # this function into forms appropriate for _align. 252 keywds = keywds.copy() 253 if len(args) != len(self.param_names): 254 raise TypeError("%s takes exactly %d argument (%d given)" 255 % (self.function_name, len(self.param_names), len(args))) 256 i = 0 257 while i < len(self.param_names): 258 if self.param_names[i] in [ 259 'sequenceA', 'sequenceB', 260 'gap_A_fn', 'gap_B_fn', 'match_fn']: 261 keywds[self.param_names[i]] = args[i] 262 i += 1 263 elif self.param_names[i] == 'match': 264 assert self.param_names[i + 1] == 'mismatch' 265 match, mismatch = args[i], args[i + 1] 266 keywds['match_fn'] = identity_match(match, mismatch) 267 i += 2 268 elif self.param_names[i] == 'match_dict': 269 keywds['match_fn'] = dictionary_match(args[i]) 270 i += 1 271 elif self.param_names[i] == 'open': 272 assert self.param_names[i + 1] == 'extend' 273 open, extend = args[i], args[i + 1] 274 pe = keywds.get('penalize_extend_when_opening', 0) 275 keywds['gap_A_fn'] = affine_penalty(open, extend, pe) 276 keywds['gap_B_fn'] = affine_penalty(open, extend, pe) 277 i += 2 278 elif self.param_names[i] == 'openA': 279 assert self.param_names[i + 3] == 'extendB' 280 openA, extendA, openB, extendB = args[i:i + 4] 281 pe = keywds.get('penalize_extend_when_opening', 0) 282 keywds['gap_A_fn'] = affine_penalty(openA, extendA, pe) 283 keywds['gap_B_fn'] = affine_penalty(openB, extendB, pe) 284 i += 4 285 else: 286 raise ValueError("unknown parameter %r" 287 % self.param_names[i]) 288 289 # Here are the default parameters for _align. Assign 290 # these to keywds, unless already specified. 291 pe = keywds.get('penalize_extend_when_opening', 0) 292 default_params = [ 293 ('match_fn', identity_match(1, 0)), 294 ('gap_A_fn', affine_penalty(0, 0, pe)), 295 ('gap_B_fn', affine_penalty(0, 0, pe)), 296 ('penalize_extend_when_opening', 0), 297 ('penalize_end_gaps', self.align_type == 'global'), 298 ('align_globally', self.align_type == 'global'), 299 ('gap_char', '-'), 300 ('force_generic', 0), 301 ('score_only', 0), 302 ('one_alignment_only', 0) 303 ] 304 for name, default in default_params: 305 keywds[name] = keywds.get(name, default) 306 value = keywds['penalize_end_gaps'] 307 try: 308 n = len(value) 309 except TypeError: 310 keywds['penalize_end_gaps'] = tuple([value] * 2) 311 else: 312 assert n == 2 313 return keywds
314
315 - def __call__(self, *args, **keywds):
316 keywds = self.decode(*args, **keywds) 317 return _align(**keywds)
318
319 - def __getattr__(self, attr):
320 return self.alignment_function(attr)
321 align = align() 322 323
324 -def _align(sequenceA, sequenceB, match_fn, gap_A_fn, gap_B_fn, 325 penalize_extend_when_opening, penalize_end_gaps, 326 align_globally, gap_char, force_generic, score_only, 327 one_alignment_only):
328 if not sequenceA or not sequenceB: 329 return [] 330 331 if (not force_generic) and isinstance(gap_A_fn, affine_penalty) \ 332 and isinstance(gap_B_fn, affine_penalty): 333 open_A, extend_A = gap_A_fn.open, gap_A_fn.extend 334 open_B, extend_B = gap_B_fn.open, gap_B_fn.extend 335 x = _make_score_matrix_fast( 336 sequenceA, sequenceB, match_fn, open_A, extend_A, open_B, extend_B, 337 penalize_extend_when_opening, penalize_end_gaps, align_globally, 338 score_only) 339 else: 340 x = _make_score_matrix_generic( 341 sequenceA, sequenceB, match_fn, gap_A_fn, gap_B_fn, 342 penalize_extend_when_opening, penalize_end_gaps, align_globally, 343 score_only) 344 score_matrix, trace_matrix = x 345 346 # print("SCORE %s" % print_matrix(score_matrix)) 347 # print("TRACEBACK %s" % print_matrix(trace_matrix)) 348 349 # Look for the proper starting point. Get a list of all possible 350 # starting points. 351 starts = _find_start( 352 score_matrix, sequenceA, sequenceB, 353 gap_A_fn, gap_B_fn, penalize_end_gaps, align_globally) 354 # Find the highest score. 355 best_score = max([x[0] for x in starts]) 356 357 # If they only want the score, then return it. 358 if score_only: 359 return best_score 360 361 tolerance = 0 # XXX do anything with this? 362 # Now find all the positions within some tolerance of the best 363 # score. 364 starts = [ 365 (score, pos) for score, pos in starts 366 if rint(abs(score - best_score)) <= rint(tolerance) 367 ] 368 369 # Recover the alignments and return them. 370 x = _recover_alignments( 371 sequenceA, sequenceB, starts, score_matrix, trace_matrix, 372 align_globally, gap_char, one_alignment_only) 373 return x
374 375
376 -def _make_score_matrix_generic( 377 sequenceA, sequenceB, match_fn, gap_A_fn, gap_B_fn, 378 penalize_extend_when_opening, penalize_end_gaps, align_globally, 379 score_only):
380 # This is an implementation of the Needleman-Wunsch dynamic 381 # programming algorithm for aligning sequences. 382 383 # Create the score and traceback matrices. These should be in the 384 # shape: 385 # sequenceA (down) x sequenceB (across) 386 lenA, lenB = len(sequenceA), len(sequenceB) 387 score_matrix, trace_matrix = [], [] 388 for i in range(lenA): 389 score_matrix.append([None] * lenB) 390 trace_matrix.append([[None]] * lenB) 391 392 # The top and left borders of the matrices are special cases 393 # because there are no previously aligned characters. To simplify 394 # the main loop, handle these separately. 395 for i in range(lenA): 396 # Align the first residue in sequenceB to the ith residue in 397 # sequence A. This is like opening up i gaps at the beginning 398 # of sequence B. 399 score = match_fn(sequenceA[i], sequenceB[0]) 400 if penalize_end_gaps[1]: 401 score += gap_B_fn(0, i) 402 score_matrix[i][0] = score 403 for i in range(1, lenB): 404 score = match_fn(sequenceA[0], sequenceB[i]) 405 if penalize_end_gaps[0]: 406 score += gap_A_fn(0, i) 407 score_matrix[0][i] = score 408 409 # Fill in the score matrix. Each position in the matrix 410 # represents an alignment between a character from sequenceA to 411 # one in sequence B. As I iterate through the matrix, find the 412 # alignment by choose the best of: 413 # 1) extending a previous alignment without gaps 414 # 2) adding a gap in sequenceA 415 # 3) adding a gap in sequenceB 416 for row in range(1, lenA): 417 for col in range(1, lenB): 418 # First, calculate the score that would occur by extending 419 # the alignment without gaps. 420 best_score = score_matrix[row - 1][col - 1] 421 best_score_rint = rint(best_score) 422 best_indexes = [(row - 1, col - 1)] 423 424 # Try to find a better score by opening gaps in sequenceA. 425 # Do this by checking alignments from each column in the 426 # previous row. Each column represents a different 427 # character to align from, and thus a different length 428 # gap. 429 for i in range(0, col - 1): 430 score = score_matrix[row - 1][i] + gap_A_fn(row, col - 1 - i) 431 score_rint = rint(score) 432 if score_rint == best_score_rint: 433 best_score, best_score_rint = score, score_rint 434 best_indexes.append((row - 1, i)) 435 elif score_rint > best_score_rint: 436 best_score, best_score_rint = score, score_rint 437 best_indexes = [(row - 1, i)] 438 439 # Try to find a better score by opening gaps in sequenceB. 440 for i in range(0, row - 1): 441 score = score_matrix[i][col - 1] + gap_B_fn(col, row - 1 - i) 442 score_rint = rint(score) 443 if score_rint == best_score_rint: 444 best_score, best_score_rint = score, score_rint 445 best_indexes.append((i, col - 1)) 446 elif score_rint > best_score_rint: 447 best_score, best_score_rint = score, score_rint 448 best_indexes = [(i, col - 1)] 449 450 score_matrix[row][col] = best_score + \ 451 match_fn(sequenceA[row], sequenceB[col]) 452 if not align_globally and score_matrix[row][col] < 0: 453 score_matrix[row][col] = 0 454 trace_matrix[row][col] = best_indexes 455 return score_matrix, trace_matrix
456 457
458 -def _make_score_matrix_fast( 459 sequenceA, sequenceB, match_fn, open_A, extend_A, open_B, extend_B, 460 penalize_extend_when_opening, penalize_end_gaps, 461 align_globally, score_only):
462 first_A_gap = calc_affine_penalty(1, open_A, extend_A, 463 penalize_extend_when_opening) 464 first_B_gap = calc_affine_penalty(1, open_B, extend_B, 465 penalize_extend_when_opening) 466 467 # Create the score and traceback matrices. These should be in the 468 # shape: 469 # sequenceA (down) x sequenceB (across) 470 lenA, lenB = len(sequenceA), len(sequenceB) 471 score_matrix, trace_matrix = [], [] 472 for i in range(lenA): 473 score_matrix.append([None] * lenB) 474 trace_matrix.append([[None]] * lenB) 475 476 # The top and left borders of the matrices are special cases 477 # because there are no previously aligned characters. To simplify 478 # the main loop, handle these separately. 479 for i in range(lenA): 480 # Align the first residue in sequenceB to the ith residue in 481 # sequence A. This is like opening up i gaps at the beginning 482 # of sequence B. 483 score = match_fn(sequenceA[i], sequenceB[0]) 484 if penalize_end_gaps[1]: 485 score += calc_affine_penalty( 486 i, open_B, extend_B, penalize_extend_when_opening) 487 score_matrix[i][0] = score 488 for i in range(1, lenB): 489 score = match_fn(sequenceA[0], sequenceB[i]) 490 if penalize_end_gaps[0]: 491 score += calc_affine_penalty( 492 i, open_A, extend_A, penalize_extend_when_opening) 493 score_matrix[0][i] = score 494 495 # In the generic algorithm, at each row and column in the score 496 # matrix, we had to scan all previous rows and columns to see 497 # whether opening a gap might yield a higher score. Here, since 498 # we know the penalties are affine, we can cache just the best 499 # score in the previous rows and columns. Instead of scanning 500 # through all the previous rows and cols, we can just look at the 501 # cache for the best one. Whenever the row or col increments, the 502 # best cached score just decreases by extending the gap longer. 503 504 # The best score and indexes for each row (goes down all columns). 505 # I don't need to store the last row because it's the end of the 506 # sequence. 507 row_cache_score, row_cache_index = [None] * (lenA - 1), [None] * (lenA - 1) 508 # The best score and indexes for each column (goes across rows). 509 col_cache_score, col_cache_index = [None] * (lenB - 1), [None] * (lenB - 1) 510 511 for i in range(lenA - 1): 512 # Initialize each row to be the alignment of sequenceA[i] to 513 # sequenceB[0], plus opening a gap in sequenceA. 514 row_cache_score[i] = score_matrix[i][0] + first_A_gap 515 row_cache_index[i] = [(i, 0)] 516 for i in range(lenB - 1): 517 col_cache_score[i] = score_matrix[0][i] + first_B_gap 518 col_cache_index[i] = [(0, i)] 519 520 # Fill in the score_matrix. 521 for row in range(1, lenA): 522 for col in range(1, lenB): 523 # Calculate the score that would occur by extending the 524 # alignment without gaps. 525 nogap_score = score_matrix[row - 1][col - 1] 526 527 # Check the score that would occur if there were a gap in 528 # sequence A. 529 if col > 1: 530 row_score = row_cache_score[row - 1] 531 else: 532 row_score = nogap_score - 1 # Make sure it's not the best. 533 # Check the score that would occur if there were a gap in 534 # sequence B. 535 if row > 1: 536 col_score = col_cache_score[col - 1] 537 else: 538 col_score = nogap_score - 1 539 540 best_score = max(nogap_score, row_score, col_score) 541 best_score_rint = rint(best_score) 542 best_index = [] 543 if best_score_rint == rint(nogap_score): 544 best_index.append((row - 1, col - 1)) 545 if best_score_rint == rint(row_score): 546 best_index.extend(row_cache_index[row - 1]) 547 if best_score_rint == rint(col_score): 548 best_index.extend(col_cache_index[col - 1]) 549 550 # Set the score and traceback matrices. 551 score = best_score + match_fn(sequenceA[row], sequenceB[col]) 552 if not align_globally and score < 0: 553 score_matrix[row][col] = 0 554 else: 555 score_matrix[row][col] = score 556 trace_matrix[row][col] = best_index 557 558 # Update the cached column scores. The best score for 559 # this can come from either extending the gap in the 560 # previous cached score, or opening a new gap from the 561 # most previously seen character. Compare the two scores 562 # and keep the best one. 563 open_score = score_matrix[row - 1][col - 1] + first_B_gap 564 extend_score = col_cache_score[col - 1] + extend_B 565 open_score_rint, extend_score_rint = \ 566 rint(open_score), rint(extend_score) 567 if open_score_rint > extend_score_rint: 568 col_cache_score[col - 1] = open_score 569 col_cache_index[col - 1] = [(row - 1, col - 1)] 570 elif extend_score_rint > open_score_rint: 571 col_cache_score[col - 1] = extend_score 572 else: 573 col_cache_score[col - 1] = open_score 574 if (row - 1, col - 1) not in col_cache_index[col - 1]: 575 col_cache_index[col - 1] = col_cache_index[col - 1] + \ 576 [(row - 1, col - 1)] 577 578 # Update the cached row scores. 579 open_score = score_matrix[row - 1][col - 1] + first_A_gap 580 extend_score = row_cache_score[row - 1] + extend_A 581 open_score_rint, extend_score_rint = \ 582 rint(open_score), rint(extend_score) 583 if open_score_rint > extend_score_rint: 584 row_cache_score[row - 1] = open_score 585 row_cache_index[row - 1] = [(row - 1, col - 1)] 586 elif extend_score_rint > open_score_rint: 587 row_cache_score[row - 1] = extend_score 588 else: 589 row_cache_score[row - 1] = open_score 590 if (row - 1, col - 1) not in row_cache_index[row - 1]: 591 row_cache_index[row - 1] = row_cache_index[row - 1] + \ 592 [(row - 1, col - 1)] 593 594 return score_matrix, trace_matrix
595 596
597 -def _recover_alignments(sequenceA, sequenceB, starts, 598 score_matrix, trace_matrix, align_globally, 599 gap_char, one_alignment_only):
600 # Recover the alignments by following the traceback matrix. This 601 # is a recursive procedure, but it's implemented here iteratively 602 # with a stack. 603 lenA, lenB = len(sequenceA), len(sequenceB) 604 tracebacks = [] # list of (seq1, seq2, score, begin, end) 605 in_process = [] # list of ([same as tracebacks], prev_pos, next_pos) 606 607 # sequenceA and sequenceB may be sequences, including strings, 608 # lists, or list-like objects. In order to preserve the type of 609 # the object, we need to use slices on the sequences instead of 610 # indexes. For example, sequenceA[row] may return a type that's 611 # not compatible with sequenceA, e.g. if sequenceA is a list and 612 # sequenceA[row] is a string. Thus, avoid using indexes and use 613 # slices, e.g. sequenceA[row:row+1]. Assume that client-defined 614 # sequence classes preserve these semantics. 615 616 # Initialize the in_process stack 617 for score, (row, col) in starts: 618 if align_globally: 619 begin, end = None, None 620 else: 621 begin, end = None, -max(lenA - row, lenB - col) + 1 622 if not end: 623 end = None 624 # Initialize the in_process list with empty sequences of the 625 # same type as sequenceA. To do this, take empty slices of 626 # the sequences. 627 in_process.append( 628 (sequenceA[0:0], sequenceB[0:0], score, begin, end, 629 (lenA, lenB), (row, col))) 630 if one_alignment_only: 631 break 632 while in_process and len(tracebacks) < MAX_ALIGNMENTS: 633 seqA, seqB, score, begin, end, prev_pos, next_pos = in_process.pop() 634 prevA, prevB = prev_pos 635 if next_pos is None: 636 prevlen = len(seqA) 637 # add the rest of the sequences 638 seqA = sequenceA[:prevA] + seqA 639 seqB = sequenceB[:prevB] + seqB 640 # add the rest of the gaps 641 seqA, seqB = _lpad_until_equal(seqA, seqB, gap_char) 642 643 # Now make sure begin is set. 644 if begin is None: 645 if align_globally: 646 begin = 0 647 else: 648 begin = len(seqA) - prevlen 649 tracebacks.append((seqA, seqB, score, begin, end)) 650 else: 651 nextA, nextB = next_pos 652 nseqA, nseqB = prevA - nextA, prevB - nextB 653 maxseq = max(nseqA, nseqB) 654 ngapA, ngapB = maxseq - nseqA, maxseq - nseqB 655 seqA = sequenceA[nextA:nextA + nseqA] + gap_char * ngapA + seqA 656 seqB = sequenceB[nextB:nextB + nseqB] + gap_char * ngapB + seqB 657 prev_pos = next_pos 658 # local alignment stops early if score falls < 0 659 if not align_globally and score_matrix[nextA][nextB] <= 0: 660 begin = max(prevA, prevB) 661 in_process.append( 662 (seqA, seqB, score, begin, end, prev_pos, None)) 663 else: 664 for next_pos in trace_matrix[nextA][nextB]: 665 in_process.append( 666 (seqA, seqB, score, begin, end, prev_pos, next_pos)) 667 if one_alignment_only: 668 break 669 670 return _clean_alignments(tracebacks)
671 672
673 -def _find_start(score_matrix, sequenceA, sequenceB, gap_A_fn, gap_B_fn, 674 penalize_end_gaps, align_globally):
675 # Return a list of (score, (row, col)) indicating every possible 676 # place to start the tracebacks. 677 if align_globally: 678 starts = _find_global_start( 679 sequenceA, sequenceB, score_matrix, gap_A_fn, gap_B_fn, penalize_end_gaps) 680 else: 681 starts = _find_local_start(score_matrix) 682 return starts
683 684
685 -def _find_global_start(sequenceA, sequenceB, 686 score_matrix, gap_A_fn, gap_B_fn, penalize_end_gaps):
687 # The whole sequence should be aligned, so return the positions at 688 # the end of either one of the sequences. 689 nrows, ncols = len(score_matrix), len(score_matrix[0]) 690 positions = [] 691 # Search all rows in the last column. 692 for row in range(nrows): 693 # Find the score, penalizing end gaps if necessary. 694 score = score_matrix[row][ncols - 1] 695 if penalize_end_gaps[1]: 696 score += gap_B_fn(ncols, nrows - row - 1) 697 positions.append((score, (row, ncols - 1))) 698 # Search all columns in the last row. 699 for col in range(ncols - 1): 700 score = score_matrix[nrows - 1][col] 701 if penalize_end_gaps[0]: 702 score += gap_A_fn(nrows, ncols - col - 1) 703 positions.append((score, (nrows - 1, col))) 704 return positions
705 706
707 -def _find_local_start(score_matrix):
708 # Return every position in the matrix. 709 positions = [] 710 nrows, ncols = len(score_matrix), len(score_matrix[0]) 711 for row in range(nrows): 712 for col in range(ncols): 713 score = score_matrix[row][col] 714 positions.append((score, (row, col))) 715 return positions
716 717
718 -def _clean_alignments(alignments):
719 # Take a list of alignments and return a cleaned version. Remove 720 # duplicates, make sure begin and end are set correctly, remove 721 # empty alignments. 722 unique_alignments = [] 723 for align in alignments: 724 if align not in unique_alignments: 725 unique_alignments.append(align) 726 i = 0 727 while i < len(unique_alignments): 728 seqA, seqB, score, begin, end = unique_alignments[i] 729 # Make sure end is set reasonably. 730 if end is None: # global alignment 731 end = len(seqA) 732 elif end < 0: 733 end = end + len(seqA) 734 # If there's no alignment here, get rid of it. 735 if begin >= end: 736 del unique_alignments[i] 737 continue 738 unique_alignments[i] = seqA, seqB, score, begin, end 739 i += 1 740 return unique_alignments
741 742
743 -def _pad_until_equal(s1, s2, char):
744 # Add char to the end of s1 or s2 until they are equal length. 745 ls1, ls2 = len(s1), len(s2) 746 if ls1 < ls2: 747 s1 = _pad(s1, char, ls2 - ls1) 748 elif ls2 < ls1: 749 s2 = _pad(s2, char, ls1 - ls2) 750 return s1, s2
751 752
753 -def _lpad_until_equal(s1, s2, char):
754 # Add char to the beginning of s1 or s2 until they are equal 755 # length. 756 ls1, ls2 = len(s1), len(s2) 757 if ls1 < ls2: 758 s1 = _lpad(s1, char, ls2 - ls1) 759 elif ls2 < ls1: 760 s2 = _lpad(s2, char, ls1 - ls2) 761 return s1, s2
762 763
764 -def _pad(s, char, n):
765 # Append n chars to the end of s. 766 return s + char * n
767 768
769 -def _lpad(s, char, n):
770 # Prepend n chars to the beginning of s. 771 return char * n + s
772 773 _PRECISION = 1000 774 775
776 -def rint(x, precision=_PRECISION):
777 return int(x * precision + 0.5)
778 779
780 -class identity_match(object):
781 """identity_match([match][, mismatch]) -> match_fn 782 783 Create a match function for use in an alignment. match and 784 mismatch are the scores to give when two residues are equal or 785 unequal. By default, match is 1 and mismatch is 0. 786 787 """
788 - def __init__(self, match=1, mismatch=0):
789 self.match = match 790 self.mismatch = mismatch
791
792 - def __call__(self, charA, charB):
793 if charA == charB: 794 return self.match 795 return self.mismatch
796 797
798 -class dictionary_match(object):
799 """dictionary_match(score_dict[, symmetric]) -> match_fn 800 801 Create a match function for use in an alignment. score_dict is a 802 dictionary where the keys are tuples (residue 1, residue 2) and 803 the values are the match scores between those residues. symmetric 804 is a flag that indicates whether the scores are symmetric. If 805 true, then if (res 1, res 2) doesn't exist, I will use the score 806 at (res 2, res 1). 807 808 """
809 - def __init__(self, score_dict, symmetric=1):
810 self.score_dict = score_dict 811 self.symmetric = symmetric
812
813 - def __call__(self, charA, charB):
814 if self.symmetric and (charA, charB) not in self.score_dict: 815 # If the score dictionary is symmetric, then look up the 816 # score both ways. 817 charB, charA = charA, charB 818 return self.score_dict[(charA, charB)]
819 820
821 -class affine_penalty(object):
822 """affine_penalty(open, extend[, penalize_extend_when_opening]) -> gap_fn 823 824 Create a gap function for use in an alignment. 825 826 """
827 - def __init__(self, open, extend, penalize_extend_when_opening=0):
828 if open > 0 or extend > 0: 829 raise ValueError("Gap penalties should be non-positive.") 830 self.open, self.extend = open, extend 831 self.penalize_extend_when_opening = penalize_extend_when_opening
832
833 - def __call__(self, index, length):
834 return calc_affine_penalty( 835 length, self.open, self.extend, self.penalize_extend_when_opening)
836 837
838 -def calc_affine_penalty(length, open, extend, penalize_extend_when_opening):
839 if length <= 0: 840 return 0 841 penalty = open + extend * length 842 if not penalize_extend_when_opening: 843 penalty -= extend 844 return penalty
845 846 863 864
865 -def format_alignment(align1, align2, score, begin, end):
866 """format_alignment(align1, align2, score, begin, end) -> string 867 868 Format the alignment prettily into a string. 869 870 """ 871 s = [] 872 s.append("%s\n" % align1) 873 s.append("%s%s\n" % (" " * begin, "|" * (end - begin))) 874 s.append("%s\n" % align2) 875 s.append(" Score=%g\n" % score) 876 return ''.join(s)
877 878 879 # Try and load C implementations of functions. If I can't, 880 # then just ignore and use the pure python implementations. 881 try: 882 from .cpairwise2 import rint, _make_score_matrix_fast 883 except ImportError: 884 pass 885 886
887 -def _test():
888 """Run the module's doctests (PRIVATE).""" 889 print("Running doctests...") 890 import doctest 891 doctest.testmod(optionflags=doctest.IGNORE_EXCEPTION_DETAIL) 892 print("Done")
893 894 if __name__ == "__main__": 895 _test() 896