4

コンテキストを少し説明すると、sam ファイルを bam に変換しようとしています。

samtools view -bT reference.fasta sequences.sam > sequences.bam

次のエラーで終了します

[E::sam_parse1] CIGAR and query sequence are of different length
[W::sam_read1] parse error at line 102
[main_samview] truncated file

問題のある行は次のようになります。

SRR808297.2571281       99      gi|309056|gb|L20934.1|MSQMTCG   747     80      101M    =       790     142     TTGGTATAAAATTTAATAATCCCTTATTAATTAATAAACTTCGGCTTCCTATTCGTTCATAAGAAATATTAGCTAAACAAAATAAACCAGAAGAACAT      @@CFDDFD?HFDHIGEGGIEEJIIJJIIJIGIDGIGDCHJJCHIGIJIJIIJJGIGHIGICHIICGAHDGEGGGGACGHHGEEEFDC@=?CACC>CCC      NM:i:2  MD:Z:98A1A

私のシーケンスの長さは 98 文字ですが、CIGAR で 101 と報告された sam ファイルを作成する際のバグの可能性があります。いくつかの読み取りを失う余裕があり、現時点では sam ファイルを生成したソース コードにアクセスできないため、バグを見つけてアライメントを再実行する機会はありません。言い換えれば、先に進むには実用的な解決策が必要です (今のところ)。したがって、ヌクレオチドの文字列の長さを数え、それを CIGAR に登録されているものと比較し、「正常な」行を新しいファイルに保存する Python スクリプトを考案しました。

#!/usr/bin/python
import itertools
import cigar

with open('myfile.sam', 'r') as f:
    for line in itertools.islice(f,3,None): #Loop through the file and skip the first three lines
            cigar=line.split("\t")[5]
            cigarlength = len(Cigar(cigar)) #Use module Cigar to obtain the length reported in the CIGAR string
            seqlength = len(line.split("\t")[9])

            if (cigarlength == seqlength):
                    ...Preserve the line in a new file...

ご覧のとおり、CIGAR を長さを示す整数に変換するために、モジュールCIGARを使用しています。正直なところ、私はその行動に少し警戒しています。このモジュールは、非常に明白なケースで長さを誤って計算しているようです。CIGAR をシーケンスの長さに変換するための別のモジュールまたはより明確な戦略はありますか?

補足:興味深いことに、控えめに言っても、この問題は広く報告されていますが、インターネットで実用的な解決策が見つかりません。以下のリンクを参照してください。

https://github.com/COMBINE-lab/RapMap/issues/9
http://seqanswers.com/forums/showthread.php?t=67253
http://seqanswers.com/forums/showthread.php?t=21120
https://groups.google.com/forum/#!msg/snap-user/FoDsGeNBDE0/nRFq-GhlAQAJ
4

2 に答える 2

2

この問題を解決するツールがないのは、この問題を示さないソフトウェアを使用してアライメントを再度実行する以外に、一般的な解決策がないためだと思います。あなたの例では、クエリ シーケンスは参照に完全に一致するため、その場合、CIGAR 文字列はあまり興味深いものではありません (Mクエリ全体の長さがプレフィックスとして付けられた単一の atch 操作のみ)。その場合、修正するには単に に変更101Mする必要があり98Mます。

ただし、より複雑な CIGAR 文字列 (挿入、削除、またはその他の操作を含むものなど) の場合I、 CIGAR 文字列のどの部分が長すぎるかを知る方法はありません。DCIGAR 文字列の間違った部分を差し引くと、ミスアライメントの読み取りが残ります。これは、読み取り全体をそのままにしておくよりも、下流の分析にとっておそらく悪いことです。

そうは言っても、それを正しく行うのがたまたま簡単な場合 (おそらく、壊れたアライメント手順により、最初または最後の CIGAR 操作に余分な塩基が常に追加される可能性があります)、知っておく必要があるのは、クエリの長​​さを計算する正しい方法です。 CIGAR 文字列。これにより、何を差し引くかがわかります。

samtoolshtslib関数bam_cigar2qlenを使用してこれを計算します。

呼び出すその他の関数はsam.hbam_cigar2qlenで定義されており、操作がクエリ シーケンスと参照シーケンスを消費する真理値表を示す役立つコメントが含まれています。

要するに、samtools (実際には htslib) が行う方法で CIGAR 文字列のクエリの長​​さを計算するには、CIGAR 操作、、、、またはの指定された長さを追加し、M他の操作の CIGAR 操作の長さを無視する必要があります。 .IS=X

python cigar モジュールの現在のバージョンは、同じ一連の操作を使用しているようで、クエリの長​​さを計算するアルゴリズム (len(Cigar(cigar))返されるもの) は私には正しいように見えます。正しい結果が得られていないと考える理由は何ですか?

orメソッドmask_leftを.mask_rightmask="H"

于 2016-10-02T01:18:55.617 に答える