用生物技术更快地计算对齐相同位点的百分比

Question

我开发了以下代码来计算对齐中相同站点的数量。不幸的是，代码是缓慢的，我必须对数百个文件进行迭代，处理1000多个对齐需要花费近12个小时，这意味着更快10倍的代码是合适的。如能提供任何帮助，将不胜感激：

import os
from Bio import SeqIO
from Bio.Seq import Seq
from Bio import AlignIO
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import generic_dna
from Bio.Align import MultipleSeqAlignment
import time

a = SeqRecord(Seq("CCAAGCTGAATCAGCTGGCGGAGTCACTGAAACTGGAGCACCAGTTCCTAAGAGTTCCTTTCGAGCACTACAAGAAGACGATTCGCGCGAACCACCGCAT", generic_dna), id="Alpha")
b = SeqRecord(Seq("CGAAGCTGACTCAGTGGGCGGAGTCACTGAAACTGGAGCACCAGTTCCTCAGAGTCCCCTTCGAGCACTACAAGAAGACAATTCGTGCGAACCACCGCAT", generic_dna), id="Beta")
c = SeqRecord(Seq("CGAAGCTGACTCAGTTGGCAGAATCACTGAAACTGGAGCACCAGTTCCTCAGAGTCCCCTTCGAGCACTACAAGAAGACGATTCGTGCGAACCACCGCAT", generic_dna), id="Gamma")
d = SeqRecord(Seq("CGAAGCTGACTCAGTTGGCAGAGTCACTGAAACTGGAGCACCAGTTCCTCAGAGTCCCCTTCGAGCACTACAAGAAGACGATTCGTGCGAACCACCGCAT", generic_dna), id="Delta")
e = SeqRecord(Seq("CGAAGCTGACTCAGTTGGCGGAGTCACTGAAACTGGAGCACCAGTTCCTCAGAGTCCCCTTCGAGCACTACAAGAAGACGATTCGTGCGAACCACCGCAT", generic_dna), id="Epsilon")

align = MultipleSeqAlignment([a, b, c], annotations={"tool": "demo"})

start_time = time.time()
if len(align) != 1:
    for n in range(0,len(align[0])):
        n=0
        i=0
        while n<len(align[0]):    #part that needs to be faster
            column = align[:,n]
            if (column == len(column) * column[0]) == True:
                i=i+1
            n=n+1

    match = float(i)
    length = float(n)
    global_identity = 100*(float(match/length))
    print(global_identity)

print("--- %s seconds ---" % (time.time() - start_time))

Answer 1

那么，您要检查5个字符串中的每个字符串在列中是否都有相同的字符？如果列中的字符都匹配，则增加i，否则增加n。你对代码的解释是正确的。

基于上述情况，我建议使用下面的代码作为更快的替代方案。

我认为align是这样的结构：

align = [
  'AGCTCGCGGAGGCGCTGCT....',
  'ACCTCGGAGGGCTGCTGTAC...',
  'AGCTCGGAGGGCTGCTGTAC...',
  # possibly more ...
]

我们尝试检测其中相同字符的列。上面，第一列是AAA (匹配)，下一列是GCG (不匹配)。

def all_equal(items):
    """Returns True iff all items are equal."""
    first = items[0]
    return all(x == first for x in items)


def compute_match(aligned_sequences):
    """Returns the ratio of same-character columns in ``aligned_sequences``.

    :param aligned_sequences: a list of strings or equal length.
    """
    match_count = 0
    mismatch_count = 0
    for chars in zip(*aligned_sequences):
        # Here chars is a column of chars, 
        # one taken from each element of aligned_sequences.
        if all_equal(chars):
            match_count += 1
        else:
            mismatch_count += 1
    return float(match_count) / float(mismatch_count)
    # What would make more sense:
    # return float(matches) / len(aligned_sequences[0])

更短的版本：

def compute_match(aligned_sequences):
    match_count = sum(1 for chars in zip(*aligned_sequences) if all_equal(chars))
    total = len(aligned_sequences[0])
    mismatch_count = total - match_count  # Obviously.
    return ...