BAM file (.bam)から、reference sequencesにmapされなかったreadsを抽出する方法について。以下、ちょっとした覚え書き。
BAM fileはbwaなどでreference sequence(s)にNGS readsをマッピングした結果のoutput file であるSAM file (.sam)を圧縮したファイル。基本的には、samtools view コマンドで作成します。
samtools view -bS example.sam > example.bam
マッピング結果から多型サイトを抽出したり、contigをexportしたり、などの作業は、BAM fileをもとに行うことが多いです。また、BAM indexを作成して、samtools tviewやTabletなどでmapping結果を見る際にも使用します。
samtools sort example.bam example.sort
samtools index example.sort.bam
samtools view example.sort.bam reference.fasta
以下は、BAM fileから"mappingされなかった" readsをFASTQ formatで抽出する方法について。
具体的には、novoalignの解説サイト
Extracting Unmapped Reads from a BAM File Produced by NovoAlign | Novocraft
を参考にしています。ていうか、以下はこのサイトの（ほぼ）日本語訳です。
-
-
-
-
- -
-
-
-
以下ここから
Novoalignで作られたBAM fileからunmapped readsを抽出する (extracting)方法について
Introduction
Referenceにmapされなかった配列をBAM fileから抽出することがが必要となるケースはしばしば生じる。Samtoolsを使うことで、この作業は容易に行うことができる。
（中略）
この作業は以下の2つのステップで行われる。
1. Unmapped readsを、readnameでsortされたBAM fileから抽出する
2. BAM fileをFASTQ read fileに変換する
FASTQ fileの生成をskipして、(unmapped) readsを直接、BAM fileから再度アラインメントすることも可能である。BAMをFASTQに変換するために、Hydra packageの"bamToFastq"を使用する。
（中略）
Novoalignによるアラインメント
# Illumina paired-end: mean=500, sd=100 novoalign -d database.nix -f reads1.fastq -reads2.fastq -oSAM -i 500 100 | samtools view -bS - alignments
Unmapped readsに対するフィルタリング
次に、以下3つのカテゴリーでreadsのフィルタリングを行う。（注：ここではIllumina NGSによるpaired end readsを想定している）
1. Paired endの一方はマッピングされているが、それ自体はマッピングされていない
2. それ自体はマッピングされているが、paired endのもう一方はマッピングされていない
3. Paired endの両方ともマッピングされていない
1-3の各カテゴリは、以下のsamtools filtering commandで抽出することができる。
samtools view -u -f 4 -F264 alignments.bam > tmps1.bam samtools view -u -f 8 -F 260 alignments.bam > tmps2.bam samtools view -u -f 12 -F 256 alignments.bam > tmps3.bam
上のコマンドは、sortされたBAM fileでも、sortされていないBAM fileでも可能である。
次に、これら各カテゴリのBAM fileを結合 (merge)し、readsの名前順にsortする。ここでのsortはpaired endの相手を正しく得るために必要である。
samtools merge -u - tmps[123].bam | samtools sort -n - unmapped
ReadsをBAM fileから直接再アラインメントする
novoalign -F BAMPE -f unmapped.bam -d …
または
ReadsをFASTQ formatで抽出する
次に、unmapped read pairsをそれぞれのFASTQ filesに抽出する（注：paired end readsを想定）。ここではHydra-SV packageのbamToFastqプログラムを利用する。
bamToFastq -bam unmapped.bam -fq1 unmapped_reads1.fastq -fq2 unmapped_reads2.fastq
Single endの場合はもっと簡単に、以下の1行コマンドで抽出できる。
samtools view -uf 4 alignments.bam | bamToFastq -bam stdin -fq1 unmapped.fastq -fq2 empty.fastq
ここまで
-
-
-
-
- -
-
-
-
わからないのは、3つのカテゴリをsamtoolsで抽出するところの
samtools view -u -f 4 -F 264 alignments.bam > tmps1.bam samtools view -u -f 8 -F 260 alignments.bam > tmps2.bam samtools view -u -f 12 -F 256 alignments.bam > tmps3.bam
の部分ですが、samtoolsにおけるoption -fや-Fを理解するためには、SAM fileの形式を理解する必要がありそうです。SAM fileがどのように記述されているか、については、日本語で以下のサイトが参考になりました。
http://crusade1096.web.fc2.com/sam.html
SAM fileの2列目は、FLAGという形式で、readの状況を表わしているようです（2進数の10進法表示！？）。いずれにしても、samtoolsを使いこなすには、まずSAM fileの意味を理解する必要がありそうです・・・
