Numpy 和 Biopython 必须集成？

Question

例如...如果 （多序列比对）MSA 有超过 50 列且间隔少于 50% 的列，我有两个脚本。

第一次使用 BioPython 在 609 列的 16281 个序列的 MSA 中花费 4.2 秒（fasta 格式的 Pfam 的 PF00085）。 [Biopython的Multiple Sequence Alignment对象的getitem方法耗费大量时间]

第二个使用简单的 IO 生成一个带有 MSA 的 2D Numpy Array，在同一个 Alignment 中只需要 1.2 秒。

我认为对 MSA 对象使用 Numpy 方法会更有用、更快。例如，您可以使用布尔 numpy 数组来选择特定的行和列。实际上，删除和选择列（例如，删除具有超过 50% 间隙的列）非常耗时，并且在 Biopython 中没有很好地实现。我认为这对于 PDB 坐标的 nx3 numpy 数组也很有用。

我有五个想法，也许只有一两个有用：

1 - 基于 numpy 而不是 str 创建一个 Seq 和多序列比对对象 (Bio.Align.MultipleSeqAlignment)。这可能是兼容性问题......也许这不是一个好主意。我不知道。

2 - 在 Biopython 中创建一个更快的方法来从 Biopython 对象获取 numpy 数组版本。我尝试为 Multiple Sequence Alignment 对象生成 numpy 数组，但这会多次调用 getitem 方法，而且它比单独使用 Biopython 更耗时。但是，也许有更多编程技能的人可以做得更好。

3 - 为 numpy 或 scipy 创建一个模块，并为对齐和 PDB 提供 IO 支持。也许是更简单和有用的想法。

4 - 创建另一个完整的 Bio 模块，但基于 numpy.也许在 scipy 或 numpy 内部。

5 - 与想法 2 和 3 一样，创建模块和方法以在 Biopython 和 numpy 对象之间更快、更有效地兼容。

你怎么看？哪些想法更好？你有更好的主意吗？能不能做点什么？我想与 Biopython 项目合作...我认为与 numpy 的集成可能是一个好的开始。

非常感谢 ;)

P.D.：我的两个脚本... 慢，基于 Biopython：

#!/usr/bin/python2.7

from sys import argv
from Bio import AlignIO
aln = AlignIO.read(open(argv[1],"r"), "fasta")
longitud = aln.get_alignment_length()
if longitud > 150:
    corte = 0.5 * len(aln)
    j = 0
    i = 0
    while j<50 and i<longitud:    
        if aln[:,i].count("-") < corte:
            j += 1
        i += 1
    if j>=50:
        print argv[1]

并且基于numpy数组最快：

#!/usr/bin/python2.7

from sys import argv
import numpy as np

with open(argv[1],'r') as archivo:
    secuencias=[]
    identificadores=[]
    temp=[]
    for linea in archivo:
        if linea[0]=='>':
            identificadores.append(linea[1:].replace('\n',''))
            secuencias.append(list(temp))
            temp=""
        else:
            temp += linea.replace('\n','')
    secuencias.append(list(temp))

sec = np.array(secuencias[1:])
ide = np.array(identificadores)

if len(ide)>150:
    corte = len(ide) * 0.5
    if np.sum(np.sum(sec=='-',1) < corte) >= 50:
        print argv[1]

Answer 1

如果您要对 MSA 对象进行大量操作，将它们视为字符数组很有用，那么我将使用 Biopython 的 AlignIO 加载对齐，然后将其转换为 NumPy 字符数组。例如：

import numpy as nump
from Bio import AlignIO
filename = "opuntia.aln"
format = "clustal"
alignment = AlignIO.read(filename, format)
align_array = numpy.array([list(rec) for rec in alignment], numpy.character)

这个快速示例可以很容易地作为 to_array 方法添加到对齐对象中，或者包含在教程中。有用吗？

当然，您仍在为所有对象创建（Seq 对象、SeqRecord 对象、空注释字典、对齐对象等）支付开销，但这是 AlignIO 接口的缺点 - 它适用于相对较重的对象模型。这对于 FASTA 和 Clustal 等简单格式并不真正需要，但对于 Stockholm 等丰富的对齐格式更有用。