質問:
BEDファイルからクロスアラインメントを除外するにはどうすればよいですか?
SmallChess
2017-05-19 10:49:47 UTC
view on stackexchange narkive permalink

BAMファイルがあります:

  @SQ SN:chr1 LN:248956422 @ SQ SN:chrx LN:248956423ST-E00110:348:HGVKKALXX:1:1201:5822:48670 323 chr1 9999 0 67H66M16H chrx 1000 0 GATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC JJJJJJJJJJJJJJJJAJJJJJJJJJJJJFJJJJJJFJFJJJJJJFJJJJJJJJJJJA77FJFJJJJJJFJFJJJJJJFJJJJJJJJJJJA77FJFJJJ NM:i:0 MD:Z RG:Z:g1  

chr1 にアラインされた読み取りがあり、 chrx にアラインされたメイトです。

BEDファイルがあります:

  chr1 0 100000 TestOnly  

BEDの範囲外のすべてを除外したいクロスアラインメントを含む領域。私の例では、読み取りは chr1 に揃えられていますが、相手はそうではありません。これを読みたくない。

読むとき:

samtools view -L test.bed test.bam

このコマンドはクロスアライメントをチェックしないため、読み取りが行われます。

私の解決策:

samtools view -L test.bed test.bam | grep -v chrx

しかし、これは非常に遅く、不器用です。私の本番パイプラインでは、次のようなことを行う必要があります。

samtools view -L test.bed test.bam | grep -v chrx | grep -v ... | grep -v ... | grep -v ... | grep -v ...

Q:より良い解決策はありますか?

1 回答:
#1
+6
terdon
2017-05-19 22:44:29 UTC
view on stackexchange narkive permalink

SAM仕様によると、SAM行の3番目のフィールド( RNAME )は次のとおりです。

RNAME:参照シーケンスNAME配置の。 @SQヘッダー行が存在する場合、RNAME(「*」でない場合)はSQ-SNタグの1つに存在する必要があります。座標のないマップされていないセグメントには、このフィールドに「*」があります。ただし、マップされていないセグメントは、ソート後に目的の位置に配置できるように、通常の座標を持っている場合もあります。 RNAMEが「*」の場合、POSとCIGARについて仮定することはできません。

そして7番目のフィールドは(私の強調、彼らの「to」が欠落しています):

RNEXT:テンプレートで読み取られたNEXTの一次アラインメントの参照配列名。最後の読み取りの場合、次の読み取りはテンプレートの最初の読み取りです。 @SQヘッダー行が存在する場合、RNEXT(「*」または「=」でない場合)はSQ-SNタグの1つに存在する必要があります。このフィールドは、情報が利用できない場合は「*」に設定され、 RNEXTが同一のRNAMEである場合は「=」に設定されます。 「=」ではなく、テンプレートの次の読み取りに1つのプライマリマッピングがある場合(FLAGのビット0x100も参照)、このフィールドは次の読み取りのプライマリ行のRNAMEと同じです。 RNEXTが '*'の場合、PNEXTおよびビット0x20については仮定できません

したがって、7番目のフィールドが = ではない行を削除する必要があります。また、念のため、7番目のフィールドが = およびではない行は3番目のフィールドと同じではありません。したがって、次のようなものを使用できます。

  samtools view -L test.bed test.bam | awk '$ 7 == "=" || $ 3 == $ 7  

そして、bamファイルとして再度保存するには:

  samtools view -L test.bed test.bam | awk '$ 7 == "=" && $ 3 == $ 7 | samtolls view -b > fixed.bam  

別の注意点として、そのような複数のgrepコマンドをチェーンする必要はほとんどありません。 \ | (または | -E または -P オプション)を使用して、それらを区切ることができます。次のようなもの:

  samtools view -L test.bed test.bam | grep -v'chrx \ | chr2 \ | chr10 \ | chrN ' 

または

  samtools view -L test.bed test.bam | grep -Ev'chrx | chr2 | chr10 | chrN ' 
このようにすると、 `fixed.bam`ファイルにヘッダーがないため、私の経験では多くの問題が発生します。常にヘッダーを追加し直すことをお勧めします。元のBAMを読み取るときに `-h`を指定するか、個別に追加します:`(samtools view -H infile.bam; samtools view…)> samtools view -b> outfile.bam`。


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