使用生物蟒蛇计算对齐的相同站点百分比的更快方法

Question

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

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 ...