Platform for Drug Discovery


SAMTOOLS,ANNOVARを用いたexome解析チュートリアル



リシークエンス概説

Exome/Target resequence は2009年終わりには現実的な方法として提案された。
全エキソンを対象とした場合を (Whole) Exome と呼び、特定のゲノム領域だけに限定した場合を Target resequence と呼ぶ場合が多いようである。
このページでは主に Exome を紹介するが、実験原理は両者とも同じであるので解析もほぼ同じになる。


Exome のメリット・デメリット

  1. データ量が Whole genome と較べてコンパクト
  2. コストが安い (2011年2月時点で GA IIx + TruSeq であればおおよそ1サンプル30-40万と思われる)
  3. 染色体の構造異常に対して弱いと考えられる (データの検証も兼ねて SNP array 実験も行うケースが多い)

Exome 実験用キット

現在 Whole exome については3社から発売されている。
また発表当時はおよそ30 Mbをターゲット領域としていたが、2011年2月時点では50 Mb程度までその領域を拡大させている。

Target resequence 実験用キット (Target capture または Custom capture と呼ばれる)

現在 Agilent と Roche でそれぞれ製品化されている。また illumina 社も類似の製品を発売する予定。
Agilent の場合は eArray(external link) により Bait oligo (いわゆるプローブ) の配置を設計することになる。


チュートリアルに使えるデータ

Roche NimbleGen の SeqCap EZ v2 (human whole exome) を用いたデータ群が SRA/DRA に収録されています。
添付してあるサンプルスクリプトでは SRA から SRR065071 (SRX026385 のデータで、約5 Gbp) をダウンロードして動作確認をしました。

関連論文


チュートリアル (擬似コードによる解析パイプラインの一例)

ここでは BWA, Picard(external link), BEDtools, SAMtools, Annovar(external link) による処理の流れを擬似コードとして示します。
  1. BWA によるアライメント (1:40:48.88 total)
  2. BWA によるマッピング
  3. SAMtools による SAM -> BAM への変換
  4. SAMtools による BAM ファイルのソート
  5. SAMtools による BAM ファイルのインデキシング
  6. Picard による duplicated reads の除去
  7. BEDtools による on-target reads の抽出
  8. SAMtools による BAM ファイルの再インデキシング
  9. BEDtools による depth of coverage の計算
  10. SAMtools による SNPs calling
  11. Annovar による common SNP の除去と nonsynonymous や遺伝子名などの注釈付け
  12. Integrated Genome Viewer で可視化

チュートリアル (上述にコード片を足したもの)

ここでは BWA, Picard, BEDtools, SAMtools, Annovar による処理の流れをコード片で示します。
実際に動作させて試してみたいという方は、ページ最下部にスクリプト全体を添付してありますので、そちらでどうぞ (メンバー登録して、ログインしてないと添付ファイルは見えないかもしれません)。


1. BWA によるアライメント

sh "#{BWA} aln -t #{CORE} hg18.fasta #{SEQ_1}.fastq > #{SEQ_1}_aln.sai"

sh "#{BWA} aln -t #{CORE} hg18.fasta #{SEQ_2}.fastq > #{SEQ_2}_aln.sai"


2. BWA によるマッピング

sh "#{BWA} sampe hg18.fasta #{SEQ_1}_aln.sai  #{SEQ_2}_aln.sai #{SEQ_1}.fastq #{SEQ_2}.fastq > #{SEQ}.aln.sam"


3. SAMtools による SAM -> BAM への変換

puts "Compress sam file..."

sh "#{MGZIP} -t #{CORE} -c #{SEQ}.aln.sam > #{SEQ}.aln.sam.gz"

puts "Convert file from sam to bam..."

sh "#{SAMTOOLS} view -q#{SAM_MAPQ} -bt hg18.fasta.fai -o #{SEQ}.aln.bam #{SEQ}.aln.sam.gz"


4. SAMtools による BAM ファイルのソート

puts "Sort bam file..."

sh "#{SAMTOOLS} sort #{SEQ}.aln.bam #{SEQ}.aln.sorted"

puts "=> #{SEQ}.aln.sorted.bam"


5. SAMtools による BAM ファイルのインデキシング

puts "Index bam file..."

sh "#{SAMTOOLS} index #{SEQ}.aln.sorted.bam"

puts "=> #{SEQ}.aln.sorted.bam.bai"


6. Picard による duplicated reads の除去

sh "java -jar -Xmx4g #{PICARD_MARK_DUPLICATES} I=#{SEQ}.aln.sorted.bam  O=#{SEQ}.aln.sorted.removedup.bam METRICS_FILE=jeter.metrics REMOVE_DUPLICATES=true ASSUME_SORTED=true VALIDATION_STRINGENCY=LENIENT"

puts "Index bam file...(2nd)"

sh "#{SAMTOOLS} index #{SEQ}.aln.sorted.removedup.bam"

puts "=> #{SEQ}.aln.sorted.removedup.bam.bai"


7. BEDtools による on-target reads の抽出

puts "Estimate overlapping rate..."

sh "#{INTERSECTBED} -wa -abam #{SEQ}.aln.sorted.removedup.bam -b #{EXOME_TARGET_REGIONS_BED} > #{SEQ}.aln.sorted.removedup.ontarget.bam"

puts "=> #{SEQ}.aln.sorted.removedup.ontarget.bam"


8. SAMtools による BAM ファイルの再インデキシング

puts "Index bam file...(3rd)"

sh "#{SAMTOOLS} index #{SEQ}.aln.sorted.removedup.ontarget.bam"

puts "=> #{SEQ}.aln.sorted.removedup.ontarget.bam.bai"


9. BEDtools による depth of coverage の計算

sh "#{COVERAGEBED} -abam #{SEQ}.aln.sorted.removedup.ontarget.bam -b #{EXOME_TARGET_REGIONS_BED} > #{SEQ}.aln.sorted.removedup.ontarget.bam.coverage"


10. SAMtools による SNPs calling

puts "Pileup..."

sh "#{SAMTOOLS} pileup -vcf hg18.fasta #{SEQ}.aln.sorted.removedup.ontarget.bam > #{SEQ}.aln.sorted.removedup.ontarget.bam.raw.pileup"

puts "=> #{SEQ}.aln.sorted.removedup.ontarget.bam.raw.pileup"



puts "varFilter..."

sh "#{SAMTOOLSPL} varFilter #{SEQ}.aln.sorted.removedup.ontarget.bam.raw.pileup | awk '($3==\"*\"&&$6>=#{PILEUP_BASEQ_INDEL})||($3!=\"*\"&&$6>=#{PILEUP_BASEQ})' > #{SEQ}.aln.sorted.removedup.ontarget.bam.filtered.pileup"

puts "=> #{SEQ}.aln.sorted.removedup.ontarget.bam.filtered.pileup"


11. Annovar による common SNP の除去と nonsynonymous や遺伝子名などの注釈付け

puts "Change variation file format from pileup to annovar..."

sh "#{CONVERT2ANNOVAR} -format pileup #{SEQ}.aln.sorted.removedup.ontarget.bam.filtered.pileup > #{SEQ}.aln.sorted.removedup.ontarget.bam.filtered.pileup.avlist"

puts "=> #{SEQ}.aln.sorted.removedup.ontarget.bam.filtered.pileup.avlist"



puts "Annotate SNPs...(filter-based)"

sh "#{ANNOTATE_VARIATION} -filter -dbtype snp130 #{SEQ}.aln.sorted.removedup.ontarget.bam.filtered.pileup.avlist #{ANNOVAR_HUMANDB}"

puts "=> #{SEQ}.aln.sorted.removedup.ontarget.bam.filtered.pileup.avlist.hg18_snp130_filtered"



puts "Annotate SNPs...(gene-based)"

sh "#{ANNOTATE_VARIATION} -geneanno #{SEQ}.aln.sorted.removedup.ontarget.bam.filtered.pileup.avlist.hg18_snp130_filtered #{ANNOVAR_HUMANDB}"


12. Integrated Genome Viewer で可視化


その他

  • 上で示したものは BWA + SAMtools でしたが Broad inst. で開発された GATK (The Genome Analysis Toolkit)(external link) というツールもよく用いられています。

  • SRR065071 を1台のマシン (Xeon X5680 @ 3.33GHz 12コア、48Gメモリ) で解析するとおおよそ3.5時間かかりました。
Analytical time
rake 02visualizeQualities_fastqc > SRR065071.log 2>& 1  861.00s user 13.03s system 103% cpu 14:07.27 total

rake 02visualizeQualities_fastx >> SRR065071.log 2>& 1  736.29s user 5.34s system 99% cpu 12:22.09 total

rake 03aligment >> SRR065071.log 2>& 1  42371.01s user 268.57s system 704% cpu 1:40:48.88 total

rake 04prepareSNPcalling >> SRR065071.log 2>& 1  2316.37s user 19.73s system 116% cpu 33:17.99 total

rake 05removeDuplicatedReads >> SRR065071.log 2>& 1  863.25s user 79.31s system 126% cpu 12:26.82 total

rake 06estimateOverlapRate >> SRR065071.log 2>& 1  499.47s user 4.81s system 99% cpu 8:24.60 total

rake 07calcurateCoverage >> SRR065071.log 2>& 1  145.66s user 1.73s system 99% cpu 2:27.46 total

rake 08callSNPs >> SRR065071.log 2>& 1  1461.95s user 5.37s system 99% cpu 24:28.12 total

rake 09annotateSNPs >> SRR065071.log 2>& 1  196.91s user 1.24s system 99% cpu 3:18.58 total


  • SRR065068 と SRR065079 でエラーがでる件
Exception in thread "main" java.lang.RuntimeException: SAM validation error: ERROR: Record 21531654, Read name SRR065068...
Picard と BWA の仕様の相違から来るもののようで、以下を参考にして VALIDATION_STRINGENCY=LENIENT を追記しました。
http://sourceforge.net/apps/mediawiki/picard/index.php?title=Main_Page#Q:_Why_am_I_getting_errors_from_Picard_like_.22MAPQ_should_be_0_for_unmapped_read.22_or_.22CIGAR_should_have_zero_elements_for_unmapped_read.3F.22(external link)



質問やアドバイスなど

  1. このページに直接かきこむ
  2. ライフサイエンスQA(external link)上で引用しつつ質問する
  3. 作成者 (埼玉医科大学の神田将和 mkoda 日付・時間帯 saitama-med.ac.jp) に直接きく
などしてもらえるとありがたいです。
Contact us
Copyright © 2009-2017 National Institute of Genetics  [Site Policy] [Privacy Policy]