例如...我有两个脚本用于查看(多序列比对)MSA是否有超过 50 列且间隙少于 50% 的列。
第一次使用BioPython在 16281 个序列的 MSA 中需要4.2 秒,有 609 列(Pfam 的 PF00085,fasta 格式)。[Biopython的Multiple Sequence Alignment对象的getitem方法耗费大量时间]
第二个使用简单的 IO 生成带有 MSA的2D Numpy 数组,在相同的对齐中只需要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 的集成可能是一个好的开始。
非常感谢 ;)
PD:我的两个脚本...慢的,基于 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]