BAM fileからunmapped readsを抽出する（samtoolsを使用する）

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