5

たとえば... (Multiple Sequence Alignment) MSAに 50 を超える列があり、ギャップが 50% 未満であるかどうかを調べるためのスクリプトが 2 つあります。

最初にBioPythonを使用すると、609 列の 16281 シーケンスの MSA で4.2 秒かかります (fasta 形式の Pfam の PF00085)。[ Biopython の Multiple Sequence Alignment オブジェクトの getitem メソッドに時間がかかる ]

単純な IO を使用して MSA で2D Numpy Arrayを生成する 2 つ目は、同じアライメントでわずか1.2 秒しかかかりません。

MSA オブジェクトへの Numpy アプローチは、より便利で高速になると思います。たとえば、ブール型の numpy 配列を使用して、特定の行と列を選択できます。実際には、列の削除と選択 (たとえば、ギャップの 50% 以上の列を削除するため) は非常に時間がかかり、Biopython ではうまく実装されていません。これは、PDB 座標の nx3 numpy 配列にも役立つと思います。

私には 5 つのアイデアがありますが、そのうちの 1 つまたは 2 つだけが役立つ可能性があります。

1 - str の代わりに numpy に基づいてSeq および Multiple Sequence Alignment オブジェクト ( Bio.Align.MultipleSeqAlignment ) を作成します。これは互換性の問題になる可能性があります...おそらくそれは良い考えではありません。知らない。

2 - Biopython オブジェクトから numpy 配列バージョンを取得するためのより高速なメソッドを Biopython で作成します。Multiple Sequence Alignment オブジェクトの numpy 配列を生成しようとしましたが、これはgetitemメソッドを複数回呼び出すため、Biopython を単独で使用するよりも時間がかかります。しかし、もっとプログラミングのスキルがあれば、もっと良いことができるかもしれません。

3 - アライメントと PDB の IO サポートを使用して、numpy または scipy のモジュールを作成します。たぶん、よりシンプルで便利なアイデアです。

4 - 別の完全な Bio モジュールを作成しますが、numpy に基づいています。scipy または numpy の内部にある可能性があります。

5 - アイデア 2 と 3 のように、モジュールとメソッドを作成して、Biopython と numpy オブジェクト間の高速かつ効率的な互換性を実現します。

どう思いますか？どのアイデアが良いですか？もっと良いアイデアはありますか？何かできることはありますか？Biopython プロジェクトと協力したい... numpy との統合は良いスタートになると思います。

たくさんの感謝 ;）

PD: 私の 2 つのスクリプト...遅い、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]
4

1 に答える 1

3

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のような豊富な配置形式ではより便利です。

于 2012-11-25T18:43:16.150 に答える