SNP callingはゲノムの変異を調べる解析です。
この記事はコースになっており、前回はこちら、初回はこちらです。

今回の記事ではBWAによるマッピングを行います!
マッピングとはどういうステップですか?

fastqに記載されたリードをゲノムにアライメントする操作です。
アライメントとは、リードの配列をみて、ゲノム上の同じ配列の場所に並べる行為です。
本コースでは大腸がんのサンプル由来のリードという設定でした。従って、このリードをヒトゲノム配列状に並べることによって、大腸がん細胞のゲノムが構築されます。
これにより、大腸がん細胞のゲノム配列と普通のヒトゲノムの配列を比べることができます💡
実際のスクリプト

まずは答えからみてましょう。
マッピングは時間がかかるのでqsubで行きましょう💡
#! /bin/bash
#$ -cwd
#$ -V
#$ -o out.txt
#$ -e error.txt
#$ -m abe
TARGET=${ARG1}
CPU=20
DIR=analysis
REF=genome/genome.fa
# mapping
bwa index ${REF}
mkdir -p ${DIR}/bam
fastq_path `ls ${DIR}/fastq_clean/${TARGET}*R1*clean*`
out_name=${DIR}/bam/`basename ${fastq_path} | cut -d_ -f1`
bwa mem -M -t ${CPU} ${REF} ${fastq_name} `echo ${fastq_path} | awk -F '_R1' '{print $1}'`_R2.clean.fastq.gz | samtools view -@ ${CPU} -Sb - > ${out_name}.bam
samtools sort -@ ${CPU} -o ${out_name}.refsort.bam ${out_name}.bam
samtools index ${out_name}.refsort.bam

なお、上記のコードは、前回のTrimmomaticのコードから続けて動作するように記述してあります。
インデックスとは?

下記の1行は「インデックス」を行なっています。
bwa index ${REF}
リードがゲノム上のどこにアライメントされるか検索するのは、非常に時間がかかります。これはゲノムがとても長いからです。

ヒトゲノムは約3.3*10^9もの長さがあり、この中からリード(数十〜数百の長さ)と同じ領域を探す。これを全リード(数千万本)に付いて実行。大変すぎる…。
そこで、予めゲノム配列について特徴を調べて目次を作っておきます。これにより検索速度がかなり改善されます。この目次を作成する操作をインデックスすると言います。

なお、-p オプションをつけることで、インデックスファイルに名前をつけることもできます。
bwa index -p index_name ${RF}
samファイル

マッピングするとどんなファイルが出力されるの?
*.samというフォーマットのファイルが出力されます。
このファイルは「ヘッダ」と「本体」の二部構成になっています。
samファイルのヘッダ
ヘッダはsamファイルの一番最初から数十〜数千行を占める領域で、必ず@から始まります。
従って、
grep xxx.sam | ^@ > header.txt
のようにすればヘッダだけ抽出することもできます。

ヘッダには何が記載されているの?
以下はヘッダの例です。
@SQ SN:1 LN:37713152
@SQ SN:2 LN:25379070
@SQ SN:3 LN:38248663
...中略...
@SQ SN:21 LN:31148813
@SQ SN:22 LN:28976614
@SQ SN:23 LN:24400806
@SQ SN:24 LN:23682337
@SQ SN:MT LN:16714
@PG ID:bwa PN:bwa VN:0.7.17-r1198-dirty CL:bwa mem -M -t 20 genome/genome.fa result/fastq_clean/test.clean.fastq.gz
@SQの行は、リファレンスに用いたゲノム(リードではなくgenome.faの方)に関する情報を示しています。
ゲノム配列は1本のDNAではありません。例えば、ヒトは23本の染色体を2対保持している、つまり46本のDNAからなります。
@SQでは(基本的に)各1本の配列に付いて、その長さなどの情報が記載されています。
@PGの行では、そのsamファイルが生成されたときに使われたコマンドが記載されています。
samファイルの本体
以下は、ヘッダーに続く本体の例です。
E00441:179:HHNNFCCXY:3:1101:5741:1590 0 5 25408598 60 149M * 0 0 CNTGCCATAAAACATAAATAAATATAAATACATTTTCATCCAAATTAGTGCAAATTTCACATTGGAAGACCATAAATACATCTATGAAATAAAAAAATCTCTCTGGAAAAATAAAAGGTTCACAACATTAAATCTGGTTTAAATTTTCC A#AFFJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJ<FJJFJJJJJJJJJAJFAJJJJJJJJA<FJJJJJJJAJJJJJAJJ<JAJJJAFFJJFFJ<AFA-FFFFJJJ7JA-<AJAJF-AFJAJAAFJJ-7F<7 NM:i:1 MD:Z:1T147 AS:i:147 XS:i:0
E00441:179:HHNNFCCXY:3:1101:5883:1590 0 18 5275589 0 149M * 0 0 CNCCGCTGACCTCACCGCCCTCACAACCCGCAAATTAAAAGCTAGTCAAGGCTGGATTTGCTAATTATGCTAACAATGATAACTATGGTAATGTGGCTAAAAGTGCTAGCATTGCAATCATCATCATCAACTAACTCAGAAAAACACAT A#AFFJJJJFJJJJJJJJJJJJ<JJJJJFJJJJJJJJJJJJJFJJJJJJJJJJJJFJFJA7FJJJJJJJJJJJJJJJJ<JJJJJJFJJJF<FJ-<--7FFFJJFFFFAJJJJ<FAJJFJJ-JFJJ7-F7--F-7-<-<<-<AFJ--7FF NM:i:5 MD:Z:1A23G52C8C32G28 AS:i:127 XS:i:127
E00441:179:HHNNFCCXY:3:1101:6167:1590 16 23 22268452 60 98M * 0 0 CATGGGTTTTTTCAGCCTCCAGGCTCGGGCGACCAAGGAATGTAGATGGATGGAAGCTTGTGGTTCCACCACCCTGGTTCGCGGCTATCTGATGTGNG 7-AFFF7<-J-7A7A-A77-A7A7<A-JA-7-<-J<A<F-FA77F-FJFA7JF7FJ<--F--<77JA<FA<JFJFJ<A-7-A<JA-7-AAF<F-7A#A NM:i:6 MD:Z:1C13A10A6G46T15C1 AS:i:74 XS:i:0
...続く
1行がとても長いです。上記に付いて、横にスクロールして確かめてみてください。
各行には、各リードについてゲノム上のどの位置にアライメントされたか記載されています。

どうやって解読すれば良いの?
各行はタブで区切られた、およそ10列からなる構成をとっています。
各列に記載される内容は決まっており、こちらのsamtoolsのマニュアルを調べるとわかります。
基本的には、
- 1列目:リードの名前
- 2列目:FLAG (DNAのどちらの鎖にアライメントされたか等の情報)
- 3列目:アライメントされた染色体の番号
- 4列目:アライメントされた実際の位置
- 5列目:アライメントのスコア
- 6列目:ミスマッチ数やINDENのレポート
- 7列目:ペアエンド時の他方のリードの名前
- 8列目:ペアエンド時の他方のリードの場所
- 9列目:ペアエンド時のインサートの長さ
- 10列目:リードのシーケンス
- 11列目:リードのクオリティ
- 12列目以降:その他の情報
です。
2列目のFLAGについては、以下のようになっています。条件を満たすものの数値の合計値がFLAG欄に記載されます。
ペアリードがある | 1 |
両方ともマップされている | 2 |
両方ともマップされていない | 4 |
相方がマップされてない | 8 |
逆相補鎖にマップされている | 16 |
相方は逆鎖にマッピングされてる | 32 |
最初のフラグメント | 64 |
最後のフラグメント | 128 |
複数個所にマップされている | 256 |
マッピングQVが低い | 512 |

FLAGを用いることで、例えばsamファイルをレファレンスの鎖ごとに分けた2つに分割することもできます。。
6列目はCIGARと呼ばれ、下記のようになっています。
M | マッチしている塩基数 |
I | 挿入がある |
D | 欠損がある |
N | スキップ配列がある |
S | soft clipping (clipped sequences present in SEQ) |
H | hard clipping (clipped sequences NOT present in SEQ) |
P | padding (silent deletion from padded reference) |
= | sequence match |
X | sequence mismatch |
この値を利用することで、次のようにしてミスマッチ数の平均などを調べることができます。
samtools view test.bam | awk -F "nM:i:" '{print $2}' | cut -f 1 | sort | uniq -c > mismatch.txt
bamファイル

最初に紹介したコードでは、samファイルではなく、最終的に*.bamなるファイルが生成されています。
bwa mem -M -t ${CPU} ${REF} ${fastq_name} `echo ${fastq_path} | awk -F '_R1' '{print $1}'`_R2.clean.fastq.gz | samtools view -@ ${CPU} -Sb - > ${out_name}.bam
基本的にマッピングにより生成されるファイルはsamですが、出力結果をsamtoolsにパイプで渡すことで.samファイルのさらに次の段階のファイルを生成しています。
生成される*.bamファイルはsamファイルの圧縮版です。
samファイルは数十GBにも及ぶ超巨大なファイルなので、普通は圧縮した形でしか保存しません。
ただし、圧縮したbamファイルはバイナリ形式になっており、中身を見ても理解することはできません。

中身をじっくり見たいときはsamに戻しましょう!ただ、下記のようにすれば、中身を手軽に見ることもできます。巨大なファイルが生成される状況は防いでおきましょう。
samtools view -h test.bam | grep -v ^@ | less
レファレンスの取得方法

NCBIやEnsambleなどのデータベースサイトから、accession numberで検索してダウンロードするだけです。
accession numberは論文や、データベースにおける種名での検索により調べます。
レファレンスの配列は更新されていくので、常に最新のものを使用するよう注意しましょう。

最近はコマンドで簡単にNCBIからゲノム配列をDLできるとか…?
MEMオプション

BWAにはMEMを含めて、3つのモードがあります。
bwa mem -M -t ${CPU} ${REF} ${fastq_name} `echo ${fastq_path} | awk -F '_R1' '{print $1}'`_R2.clean.fastq.gz | samtools view -@ ${CPU} -Sb - > ${out_name}.bam
3つのオプションはアルゴリズムが異なり、それぞれ下記のケースに適しています。
- BWA-backtrack:ショートリードをマッピングするとき
- BWA-SWロングリードをマッピングするとき
- BWA-MEM:ロングリードのマッピングで、ギャップが多いと予想されるとき
詳細については、大元の論文やBWAのサイトを参考にしてください。
真核生物のRNA-SeqやSV検出が目的の場合は注意

イントロンのような長い塩基挿入がある場合に、BWA はマッピングに失敗してしまいます!
よって、BWA は真核生物の RNA-Seq リードのマッピングや、大きなINDELが予想される構造変異(Structual Variant)の検出には使われません。

BWAはSNPや小さめのINDELを検出するときに良く使います。
Rタグについて

以降の解析で用いるGATKというソフトウェアを使用するためには、samファイルにRタグというものが付いている必要があります。付いていないとエラーが生じて停止してしまいます。
以下は、bwaの出力によるsamファイルです(上記と同じもの)。
E00441:179:HHNNFCCXY:3:1101:5741:1590 0 5 25408598 60 149M * 0 0 CNTGCCATAAAACATAAATAAATATAAATACATTTTCATCCAAATTAGTGCAAATTTCACATTGGAAGACCATAAATACATCTATGAAATAAAAAAATCTCTCTGGAAAAATAAAAGGTTCACAACATTAAATCTGGTTTAAATTTTCC A#AFFJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJ<FJJFJJJJJJJJJAJFAJJJJJJJJA<FJJJJJJJAJJJJJAJJ<JAJJJAFFJJFFJ<AFA-FFFFJJJ7JA-<AJAJF-AFJAJAAFJJ-7F<7 NM:i:1 MD:Z:1T147 AS:i:147 XS:i:0
以下は、上記samファイルのRタグが付いたバージョンです。
E00441:179:HHNNFCCXY:3:1101:5741:1590 0 5 25408598 60 149M * 0 0 CNTGCCATAAAACATAAATAAATATAAATACATTTTCATCCAAATTAGTGCAAATTTCACATTGGAAGACCATAAATACATCTATGAAATAAAAAAATCTCTCTGGAAAAATAAAAGGTTCACAACATTAAATCTGGTTTAAATTTTCC A#AFFJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJ<FJJFJJJJJJJJJAJFAJJJJJJJJA<FJJJJJJJAJJJJJAJJ<JAJJJAFFJJFFJ<AFA-FFFFJJJ7JA-<AJAJF-AFJAJAAFJJ-7F<7 MD:Z:1T147 RG:Z:1 NM:i:1 AS:i:147 XS:i:0

最後の方のオプション部分に RG:Z:1 という部分が追加されているのが分かりますでしょうか。これがRタグです。
Rタグの付与は次回で紹介するPicardというソフトウェアによって可能ですが、下記のようにBWAで-Rオプションをつけることでも可能です。
bwa mem -M -R '@RG\tID:sample_1\tLB:sample_1\tPL:ILLUMINA\tPM:HISEQ\tSM:sample_1' ref input_1 input_2 > aligned_reads.sam

上記の方法は、SNPcallingのベストプラクティスなパイプラインの1つであるこちらのサイトで紹介されている方法です。
sortとindex

下記部分のコードはbamファイルのsortとindexという操作に相当し、SNPcallingに限らず多くの場合で必要な操作です。
samtools sort -@ ${CPU} -o ${out_name}.refsort.bam ${out_name}.bam
samtools index ${out_name}.refsort.bam
sortは、bamファイル内のリードを、染色体上の順番に並び替えるような操作に相当します。これにより、下流の解析が早くなります。また、bamファイル自体の容量も小さくなります。
indexは、bamファイルの検索が効率的になるように目次をつける操作です。これにより*.bai というファイルが生成されます。例えば、IGVなどでこの *.bai ファイルが必要になります。
結果をまとめよう

マッピングにより、数〜数十パーセントのリードがドロップアウト(アライメントされない)します。ですから、マッピングの後は下記のようなコードでリードの損失状況を記録して起きましょう。
mkdir -p rep
a1=/groups2/gca50068/tatsu/*15/bam/*refsort.bam
a2=rep/flagstat.txt
for i in `ls $a1`;do
echo =================== >> $a2
echo `basename $i` >> $a2
samtools flagstat $i >> $a2
done
less rep/$a2
flagstatは、生成されたbamファイルに関するまとめ情報を出力してくれるオプションです。falgstatにより、下記のような結果が出力されます。
116888737 + 0 in total (QC-passed reads + QC-failed reads)
1529806 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
114554248 + 0 mapped (98.00% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
最初からここまでの結果をまとめよう

次からSNPcallingのプロセスに入ります。そこで、ここまででリードがどれくらいドロップアウトしたか下記のような表を作成してまとめておきましょう。
samtoolsのインストール

以下の記事または動画で紹介してますので是非ご覧ください。
動画 ↓
BWAのインストール

git clone して make するだけです。
GitHubのサイトはこちらです。BWAのサイトの Repository のリンクからたどり着けます。後は、READMEに記載されている通りにインストールします。

今回はここまでです!お疲れ様でした!次回はSNPcallingの最初のステップであるPicardによる操作を行います。