WGSデータの参照ゲノム配列へのマッピング (4): BWAによるマッピング

2015-08-06

ヒトゲノム参照配列に対して、NGSのリード（ERR251633）をbwaを使ってマッピングします。この辺りの解析は、以前のエントリと大体同じです。

ただし、今回はmappingの前にNGSのraw dataに前処理をします。MappingするFASTQ fileに対して、qualityの低いリードを除き、長さの短い（<50bp）リードも除去します。前処理を行うことで、bwaの解析にかかる時間の短縮と、精度の向上が期待できます（150805 --> 本当に？　どれくらい？）。具体的な方法は、井上さんによる以下のサイトを参考にしています。

参考URL: Trinity - 井上 潤

必要な解析ソフト（SolexaQA, DynamicTrim, LengthSort）
SolexaQA

(1) Index fileの作成

まずリファレンス配列のindexを作成します。ここでは時間の都合上、ヒトゲノム参照配列の6番染色体を使用します。データの入手はこちらを参照。

# bwa index実行 // 3分くらいかかる
bwa index -a bwtsw Homo_sapiens.GRCh38.dna_rm.chromosome.6.fa

(2) NGS reads (.fastq) ファイルの前処理

.sra fileからFASTQ fileを抽出

以下のコマンドを実行します。ここではpaired endの右側・左側のデータをそれぞれ別のfileに分けます。

fastq-dump --split-files ERR251633.sra

DynamicTrimによるlow quality readsの除去

以下のようにコマンドを実行します。

DynamicTrim.pl ERR251633_1.fastq -h 20 &
DynamicTrim.pl ERR251633_2.fastq -h 20 &

LengthSortによる短すぎるreadsの除去

以下のようにコマンドを実行します。ここでは50bpより短いreadsを除きます。

LengthSort.pl ERR251633_1.fastq.trimmed ERR251633_2.fastq.trimmed -length 50

ここで生成される"ERR251633_1.fastq.trimmed.paired1"と"ERR251633_1.fastq.trimmed.paired2"が、以下のbwaでのmappingに使用するFASTQ fileになります。

(3) BWAによるmappingの実行

ここではbwa memを使用します。以下のコマンドを実行します。

# bwa mem実行 // 約21分 (MacPro 2 x 3.06GHz 6 core Intel Xeon, 32GB RAM)
bwa mem -M -t 12 Homo_sapiens.GRCh38.dna_rm.chromosome.6.fa ../ERR251633_1.fastq.trimmed.paired1 ../ERR251633_1.fastq.trimmed.paired2 > ERR251633.chr6.sam

(4) samtoolsによるSAM->BAM変換、sorting, indexing

以下のコマンドを順次実行します。

samtools view -bS ERR251633.chr6.sam > ERR251633.chr6.bam    # .sam -> .bam変換
samtools sort ERR251633.chr6.bam ERR251633.chr6.sort    # sorting
samtools index ERR251633.chr6.sort.bam

.bam fileはsamtools tviewコマンドや、Tabletなどのviewerを使って見ることができます。