SNPcalling例1バイオインフォ

#3【SNPcalling例1】BWAでマッピング

 

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
相方がマップされてない
逆相補鎖にマップされている 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は論文や、データベースにおける種名での検索により調べます。

レファレンスの配列は更新されていくので、常に最新のものを使用するよう注意しましょう。

ねこ研究員
ねこ研究員

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のインストール

ねこ研究員
ねこ研究員

以下の記事または動画で紹介してますので是非ご覧ください。

記事:#3【NGS】samtoolsのインストール

動画 ↓

【Linux】samtoolsのインストール実況してみた【2分】

BWAのインストール

ねこ研究員
ねこ研究員

git clone して make するだけです。

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

 

ねこ研究員
ねこ研究員

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