2017 8/29追記
2019 4/16 誤字修正
ゲノム解析の検証方法やベンチマーク方法の障害は、サンプルゲノムの突然変異の状況についての「根拠のある真実」がわかっていて完全に検証されている参照データセットがほとんどないことである。さらに、本物のヒトゲノムデータセットの無料で一般公開されている利用可能性は、ドナーのプライバシーの保護と両立しない。ゲノムデータをよりよく分析し、理解するためには、既知の生物学やシークエンシングアーティファクトを反映した、すべてのバリアントをモデル化するテストデータセットが必要である。リードシミュレータはこの要件を満たすことができるが、多くの場合、真のデータと全体的な柔軟性の欠如に類似しているとの批判がある。リードシミュレータだけでなく、バリアント比較やツール評価を容易にするためのスクリプトを含むツールのセットであるNEAT（NExt-generation sequenceing Analysis Toolkit）を紹介する。 NEATには、デフォルトモデルに手動で設定することも、実際のデータセットを使用してパラメータ化することもできる、さまざまな調整可能パラメータがある。
NEAT-genReadsは2016年に発表されたfastqをシミュレートできるツール。変異のVCFファイルなどの情報も与えて現実に近いfastqを発生させることができる。fastq以外にポジコンとして使えるbamやVCFファイルも生成されるため、indel検出ツールの妥当性をポジコンと比較して考えることが可能になっている。
NEAT-genReadsにはすでに学習済みのプロファイルが用意されているが、ユーザーがbamやVCFファイルを指定することでユーザー定義のプロファイルでシーケンスbiasやエラー発生率を再現することができる。また領域をbedファイルで明示することで、特定の領域だけシミュレートすることができる（例えばexome解析のシミュレーション）。
インストール
依存
- Python 2.7
- Numpy 1.9.1+
本体 Github
git clone https://github.com/zstephens/neat-genreads.git
cd neat-genreads/
> ./genReads.py -h
$ genReads.py -h
usage: genReads.py [-h] -r <str> -R <int> -o <str> [-c <float>] [-e <str>]
[-E <float>] [-p <int>] [-t <str>] [-to <float>] [-m <str>]
[-M <float>] [-Mb <str>] [-N <int>] [-v <str>]
[--pe <int> <int>] [--pe-model <str>] [--gc-model <str>]
[--job <int> <int>] [--nnr] [--bam] [--vcf] [--rng <int>]
[--gz] [--no-fastq]
NEAT-genReads V2.0
optional arguments:
-h, --help show this help message and exit
-r <str> * ref.fa
-R <int> * read length
-o <str> * output prefix
-c <float> average coverage
-e <str> sequencing error model
-E <float> rescale avg sequencing error rate to this
-p <int> ploidy
-t <str> bed file containing targeted regions
-to <float> off-target coverage scalar
-m <str> mutation model pickle file
-M <float> rescale avg mutation rate to this
-Mb <str> bed file containing positional mut rates
-N <int> below this qual, replace base-calls with 'N's
-v <str> input VCF file
--pe <int> <int> paired-end fragment length mean and std
--pe-model <str> empirical fragment length distribution
--gc-model <str> empirical GC coverage bias distribution
--job <int> <int> jobs IDs for generating reads in parallel
--nnr save non-N ref regions (for parallel jobs)
--bam output golden BAM file
--vcf output golden VCF file
--rng <int> rng seed value; identical RNG value should produce
identical runs of the program, so things like read
locations, variant positions, error positions, etc,
should all be the same.
--gz gzip output FQ and VCF
--no-fastq bypass fastq generation
——
実行方法
1、デフォルト条件でランすると平均カバレッジ10のfastqが作られる。
./genReads.py -r ref.fa -R 101 -o simulated_data
- -R read length. Required.
- -r ref.fa. Required.
ランが終わると"simulated_data_read1.fq" ができる。デフォルトのプロファイルに従い変異は導入されている。
２、インサートサイズが平均600で、カバレッジ x50、301bp x 2の ペアードエンドfastqを作る。アライメントのbamと変異が起きた部位を示したvcfファイルも出力させる。
./genReads.py --pe 600 50 -c 50 -R 301 -r ref.fa --bam --vcf -o simulated_data
- --vcf Output golden VCF file
- --bam Output golden BAM file
- --pe Paired-end fragment length mean and standard deviation. To produce paired end data, one of --pe or --pe-model must be specified.
- -c Average coverage across the entire dataset. Default: 10
"simulated_data_read1.fq" と "simulated_data4_read2.fq" ができる。
bamとvcfをIGVに読み込ませて確認する。(igvがない人はこちらを参照）
samtools index simulated_data_golden.bam #index作成
igv -g ref.fa simulated_data_golden.bam,simulated_data_golden.vcf #IGV起動と読み込み
一番上が変異部位。ホモとヘテロの部位がある。defaultの条件では かなりの頻度でindel変異が入る。
３、ユーザーが指定したVCFと同じ変異を導入する。また、それ以外の変異の発生率をゼロにしてfastqを作る。
./genReads.py --pe 600 50 -c 50 -R 301 -v mutation_call.vcf -M 0 -r ref.fa -o simulated_data
- -v Variants from this VCF will be inserted into the simulated sequence with 100% certainty.
- -M The mutation rate model is rescaled to make this the average value. Must be between 0 and 0.3. These random mutations are inserted in addition to the once specified in the -v option.
入力のvcfと出力のvcfをIGVに読み込んで結果を確認。
samtools index simulated_data_golden.bam #index作成
igv -g ref.fa simulated_data_golden.bam,simulated_data_golden.vcf,mutation_call.vcf #IGV起動と読み込み
導入される変異がホモかヘテロかは入力のVCFファイルの GT:AD:DP:GQ:PLのGTで決定される。GTが1/1なら変異がホモ、0/1ならヘテロ、0/0ならリファレンスのホモになる。
-参考-
ーtでbedを指定すれば、特定の領域だけシミュレートすることができます。exome解析やchip-seqのデータのシミュレートのほか、計算リソースを削減するために特定のクロモソームだけシミュレートしたい時に使えるオプションです。
vcfで情報を与えれば large indelも再現される。
large del（golden.bamをigvで開いて確認した）
４、既存データのGCバイアスに基づいてGC biasモデルを作り、それを使いfastqを作る。
bedtools genomecov -d -ibam normal.bam -g ref.fa > coverage
python2.7 computeGC.py -r ref.fa -i coverage -w 100 -o model
./genReads.py --gc-model GCcoverage -r ref.fa -R 101 -o simulated_data
はじめにbedtoolsでcoverageを計算させている（リンク）。
注意；リード長もモデルのbamと同じと仮定してシミュレートされるため、ラストのコマンドの-R の指定値は元のbamのリード長と揃える必要がある。
４と似た流れで変異、リードサイズ、シーケンスエラーのbiasを既存データからモデル化して利用することもできます。詳細はGitのHPで確認してください。
GitHub - zstephens/neat-genreads: NEAT read simulation tools
引用
Simulating Next-Generation Sequencing Datasets from Empirical Mutation and Sequencing Models
Stephens ZD, Hudson ME, Mainzer LS, Taschuk M, Weber MR, Iyer RK
PLoS One. 2016 Nov 28;11(11):e0167047
追記
リファレンスにvcfのバリアントを再現するならVCFtoolsのvcf-consensusが利用できます。
