macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

ユーザー定義の変異を再現可能なfastqのシミュレーター NEAT-genReads

2017 8/29追記 

 

NEAT-genReadsは2016年に発表されたfastqをシミュレートできるツール。変異のVCFファイルなどの情報も与えて現実に近いfastqを発生させることができる。fastq以外にポジコンとして使えるbamやVCFファイルも生成されるため、indel検出ツールの妥当性をポジコンと比較して考えることが可能になっている。

NEAT-genReadsにはすでに学習済みのプロファイルが用意されているが、ユーザーがbamやVCFファイルを指定することでユーザー定義のプロファイルでシーケンスbiasやエラー発生率を再現することができる。また領域をbedファイルで明示することで、特定の領域だけシミュレートすることができる(例えばexome解析のシミュレーション)。

 

インストール

Github

GitHub - zstephens/neat-genreads: NEAT read simulation tools

本体はpythonで記述されている。ダウンロードしてパスを通す。

 依存

 $ 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" ができる。デフォルトのプロファイルに従い変異は導入されている。

 

2、インサートサイズが平均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起動と読み込み

f:id:kazumaxneo:20170817191004j:plain

一番上が変異部位。ホモとヘテロの部位がある。defaultの条件では かなりの頻度でindel変異が入る。

 

  

3、ユーザーが指定した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も再現される。

f:id:kazumaxneo:20170829194446j:plain

large del(golden.bamをigvで開いて確認した)

 

 

4、既存データの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のリード長と揃える必要がある。

 

4と似た流れで変異、リードサイズ、シーケンスエラーの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