質問:
BAMファイルでのINDEL呼び出し/修正を高速化するにはどうすればよいですか?
gringer
2017-06-05 04:25:12 UTC
view on stackexchange narkive permalink

samtools mpileup コマンドには、INDELのミスアライメントに関連するマッピングエラーを修正できる非常に優れた機能があります。デフォルトでは、 mpileup コマンドは、リファレンスゲノムのカバレッジが250倍を超える読み取りでは機能しません。この制限を増やすことはできますが、カバレッジが非常に高いとmpileupプログラムが停止するため、これを高速化する簡単な方法があるかどうかを知っておくと便利です。

もう少し追加するにはコンテキストでは、イルミナの全ゲノムシーケンス(カバレッジ〜1000X)とIonTorrentで実行されたターゲットアンプリコンシーケンス(カバレッジ〜4000X)の両方から抽出されたミトコンドリアゲノムリードを使用してこれを行ってきました。

@rightskewedが、 samtools view -s <float> こちらを参照)でsamtoolsのダウンサンプリング機能について言及しているようです。これは、この解決策として機能する可能性があります。 mpileup操作の前に使用した場合。

1 回答:
#1
+4
gringer
2017-06-05 04:34:22 UTC
view on stackexchange narkive permalink

数年前にこの問題が発生したとき、samtoolsのサブサンプリングに気づかなかったため、マップされた読み取りを処理するための独自のデジタル正規化メソッドを作成することになりました。この方法では、ゲノムカバレッジが減少しますが、カバレッジが低い読み取りは保持されます。

IonTorrent読み取り(可変長)を使用していたため、にマップされた最長の読み取りを選択するというアイデアを思いつきました。ゲノム内の各場所(そのような読み取りが存在すると仮定)。これは、さまざまなサンプルの非常に可変的なカバレッジ(200Xの場合もあれば、4000Xの場合もある)が、約100〜200Xのはるかに一貫したカバレッジにフラット化されることを意味しました。私が書いたPerlコードのコアは次のとおりです。

  if(($ F [2] ne $ seqName)||($ F [3]!= $ pos)||(長さ($ bestSeq)< = length($ F [9]))){if(length($ bestSeq)== length($ F [9])){##リザーバーサイズ1のリザーバーサンプリング## httpsを参照://en.wikipedia.org/wiki/Reservoir_sampling ## *確率1 / iで、現在のアイテムの代わりに新しいアイテムを保持します$ seenCount ++; if(!rand($ seenCount)){##つまり、rand($ seenCount)== 0の場合、次に置換を続行します。 }} else {$ sawenCount = 1; } if(($ F [2] ne $ seqName)||($ F [3]!= $ pos)){if($ output eq "fastq"){printSeq($ bestID、$ bestSeq、$ bestQual); } elsif($ output eq "sam"){print($ bestLine); }} $ seqName = $ F [2]; $ pos = $ F [3]; $ bestLine = $ line; $ bestID = $ F [0]; $ bestFlags = $ F [1]; $ bestSeq = $ F [9]; $ bestQual = $ F [10];  

完全なコード(非圧縮SAMファイルのフィルターとして機能)はここにあります。

明確にするために:これはPerlコードですか?
はい、Perlコードです。どこでもドル記号でわかりますが、Perlに慣れていない人にはわかりにくいと思います。
私はPerlを話しますが、識別子の前にドルがある言語はPerlだけではありません。 ;-)
確かに。私の心は空白でしたが、妻は私にbash / shを思い出させました。


このQ&Aは英語から自動的に翻訳されました。オリジナルのコンテンツはstackexchangeで入手できます。これは、配布されているcc by-sa 3.0ライセンスに感謝します。
Loading...