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

 

SNP callingはゲノムの変異を調べる解析です。

この記事はコースになっており、前回はこちら、初回はこちらです。

マッピングとはどういうステップですか？

アライメントとは、リードの配列をみて、ゲノム上の同じ配列の場所に並べる行為です。

本コースでは大腸がんのサンプル由来のリードという設定でした。従って、このリードをヒトゲノム配列状に並べることによって、大腸がん細胞のゲノムが構築されます。

これにより、大腸がん細胞のゲノム配列と普通のヒトゲノムの配列を比べることができます💡

実際のスクリプト

マッピングは時間がかかるので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

bwa index ${REF}

リードがゲノム上のどこにアライメントされるか検索するのは、非常に時間がかかります。これはゲノムがとても長いからです。

そこで、予めゲノム配列について特徴を調べて目次を作っておきます。これにより検索速度がかなり改善されます。この目次を作成する操作をインデックスすると言います。

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
マッピングＱＶが低い	512

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ファイル

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ファイルはバイナリ形式になっており、中身を見ても理解することはできません。

samtools view -h test.bam | grep -v ^@ | less

レファレンスの取得方法

accession numberは論文や、データベースにおける種名での検索により調べます。

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

MEMオプション

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 は真核生物の RNA-Seq リードのマッピングや、大きなINDELが予想される構造変異(Structual Variant)の検出には使われません。

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

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

 

sortとindex

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)

最初からここまでの結果をまとめよう

 

samtoolsのインストール

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

動画 ↓

BWAのインストール

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