2019 5/13 コメントおよびパラメータ追記
2020 1/4 conda インストール追記
ますます多くのバイオインフォマティクスツールがリリースされており、特定の実験に最適なツールや最適なツールを知ることは困難になっている。ゲノミクスとメタゲノミクスのデータのシミュレーションは、実験の計画と新しい方法の開発の両方において重要な役割を担っている。実際のデータとは対照的に、シミュレートされたデータは、例えばメタゲノミクスの場合には、サンプル中に存在する種の存在量など、制御されたパラメータで生成することができる。パラメータを固定することで、制御された条件でのツールのベンチマークとテストが可能になり、新しいツールやパイプラインをテストするための模擬データも提供されます。このようなテスト環境は、メタゲノミクス（Escalona et al、2016）のような急成長中のサブフィールドにとって特に重要になる。加えて、シミュレートされたデータは授業で非常に有用であることが証明されており、教師はしばしば模擬データセットを必要としている。このようなデータセットは、迅速に分析され、学生にとって解釈しやすい有意義で明確な結果をもたらす（Halley、1991）。
驚くべきことに、メタゲノミクスのためのこのようなシミュレーションソフトウェアはほんのわずかしか存在せず、既存のソリューションはしばしば使用が困難で不便であり、保守が不十分である。ここでは、InSilicoSeqについて説明する。InSilicoSeqは、（メタ）ゲノムからの現実的なIlluminaのリードをシミュレートするソフトウェアである。 InSilicoSeqはマルチスレッド化されており、Document化されており、Pythonのパッケージマネージャーpip経由で簡単にインストールできる。 InSilicoSeqは、（メタ）ゲノムソフトウェアのベンチマークとテストを容易にすることを目指している。
InSilicoSeqはPythonで書かれており、PHRED scoresを正確にモデル化し、置換、挿入および欠失エラーをサポートするだけでなく、サイズ分布とGCバイアスを導入できる。
InSilicoSeqは、base qualityと挿入サイズをモデル化するためにKDE（Kernel Density Estimation）を実装している。簡潔に言えば、KDEはノンパラメトリックな推定値クラスであり、一般的にヒストグラムよりも分布の推定がより滑らかになる（Silverman、1986）。挿入、欠失、置換、およびGCバイアスのシミュレーションは、empirically に行われ、揃えられたリードから計算される。
現在のリリースでは、InSilicoSeqにはMiSeq、HiSeq、NovaSeq用に事前に構築されたエラーモデルが付属しているが、アラインメントされたリードを含むすべてのbamファイルからエラーモデルを生成するコマンドを提供している。提供されたエラーモデルは、MiSeq計測器用のPRJEB20178と、HiSeqおよびNovaSeq計測器用のIllumina Basespaceのパブリックデータセットの3つのデータセットから生成された、bam形式のアライメントされたリードから計算される。 3つのデータセットはMEGAHIT（Li et al。、2015）でアセンブルされ、リードはデフォルトパラメータでbowtie2（Langmead and Salzberg、2012）を使用してアセンブリにマップされた。
InSilicoSeqはメタゲノミクス用に設計されており、デフォルトの対数正規分布に基づいて複数のゲノムからのシーケンシングリードを生成する。他のディストリビューションが組み込まれているだけでなく、ソフトウェアに各入力ゲノムの正確な量を提供する可能性もある。
documentation
http://insilicoseq.readthedocs.io/en/latest/
インストール
mac os10.12のpython2.7.10でテストした。
本体 Github
conda/pipで導入できる。dockerイメージも用意されている。
#bioconda(link)
conda install -c bioconda -y insilicoseq
pip install InSilicoSeq #python2の時。python3ならpip3
> iss -h
$ iss -h
usage: iss [subcommand] [options]
InSilicoSeq: A sequencing simulator
optional arguments:
-h, --help show this help message and exit
-v, --version print software version and exit
available subcommands:
model generate an error model from a bam file
generate simulate reads from an error model
——
> iss generate -h
$ iss generate -h
usage: iss generate [-h] [--quiet | --debug] [--cpus <int>]
[--genomes <genomes.fasta> | --ncbi [<str> [<str> ...]]]
[--n_genomes [<int> [<int> ...]]]
[--abundance <str> | --abundance_file <abundance.txt> | --coverage <coverage.txt>]
[--n_reads <int>] [--mode <str>] [--model <npz>]
[--gc_bias] --output <fastq>
simulate reads from an error model
arguments:
-h, --help show this help message and exit
--quiet, -q Disable info logging. (default: False).
--debug, -d Enable debug logging. (default: False).
--cpus <int>, -p <int>
number of cpus to use. (default: 2).
--genomes <genomes.fasta>, -g <genomes.fasta>
Input genome(s) from where the reads will originate
--ncbi [<str> [<str> ...]], -k [<str> [<str> ...]]
Download input genomes from NCBI. Requires
--n_genomes/-u option. Can be bacteria, viruses,
archaea or a combination of the three (space-
separated)
--n_genomes [<int> [<int> ...]], -u [<int> [<int> ...]]
How many genomes will be downloaded from NCBI.
Required if --ncbi/-k is set. If more than one kingdom
is set with --ncbi, multiple values are necessary
(space-separated). Can also be set with --genomes/-g
in which case random genomes will be taken from the
input multifasta
--abundance <str>, -a <str>
abundance distribution (default: lognormal). Can be
uniform, halfnormal, exponential, lognormal or zero-
inflated-lognormal.
--abundance_file <abundance.txt>, -b <abundance.txt>
abundance file for coverage calculations (default:
None).
--coverage <coverage.txt>, -C <coverage.txt>
file containing coverage information (default: None).
--n_reads <int>, -n <int>
Number of reads to generate (default: 1000000). Allows
suffixes k, K, m, M, g and G (ex 0.5M for 500000).
--mode <str>, -e <str>
Error model. If not specified, using kernel density
estimation (default: kde). Can be kde or basic.
--model <npz>, -m <npz>
Error model file. (default: None). Use HiSeq, NovaSeq
or MiSeq for a pre-computed error model provided with
the software, or a file generated with iss model. If
you do not wish to use a model, use --mode basic. The
name of the built-in models is case insensitive.
--gc_bias, -c If set, may fail to sequence reads with abnormal GC
content. Does not guarantee --n_reads (default: False)
--output <fastq>, -o <fastq>
Output file prefix (Required)
——
> iss model -h
$ iss model -h
usage: iss model [-h] [--quiet] [--debug] --bam <bam> --output <npz>
generate an error model from a bam file
arguments:
-h, --help show this help message and exit
--quiet, -q Disable info logging. (default: False).
--debug, -d Enable debug logging. (default: False).
--bam <bam>, -b <bam>
aligned reads from which the model will be inferred
(Required)
--output <npz>, -o <npz>
Output file prefix (Required)
——
ラン
２つのモードがある。 iss generateはイルミナのシミュレーションデータを発生させるモードで、 iss modelは指定されたデータからエラープロファイルを構築するモードとなる。オーサーらが準備したエラーファイルで大半のケースは問題ないとされる。
iss generate: simulate reads from an error model（リンク）
Mock メタゲノムコミュニティゲノム配列をリファレンスにしてペアエンドリードを100万発生させる。
# Mock bacteria communityのゲノムデータSRS121011.fasta（リンク）のダウンロード。
curl -O -J -L https://osf.io/thser/download
#SRS121011をテンプレートにペアエンドリードを発生。
iss generate --genomes SRS121011.fasta -p 8 --model miseq -o miseq -n 1000000 --abundance lognormal
- --model <npz> Error model file. (default: None). Use HiSeq, NovaSeq or MiSeq for a pre-computed error model provided with the software
-
-p <int> number of cpus to use. (default: 2)
-
-o <fastq> Output file prefix (Required)
- -n Number of reads to generate (default: 1000000). Allows suffixes k, K, m, M, g and G (ex 0.5M for 500000)
- --abundance abundance distribution (default: lognormal). Can be uniform, halfnormal, exponential, lognormal or zero-inflated-lognormal.
デファルトではリード数の分布は対数正規分布（log normal distribution）になっている。分布は "--abundance"で変更できる（default: lognormal）。エラープロファイルに従い、塩基置換やindelも発生する。indel変異を除きたければ"--mode basic"を使う。 "--gc_bias"をつけると、GC含量が極端な配列についてGC biasが起きる。
ジョブが終わるとmiseq_reads_R1.fastq、miseq_reads_R2.fastq、miseq_reads_abundance.txtができる。今回の例ではMock communityの21ゲノムのシーケンスデータが含まれる。それぞれの存在量はレポートファイルmiseq_reads_abundance.txtにまとめられる。
存在量が0.147なら、トータルリード長X bpとして、X * 0.147がそのゲノムのトータルリード長（bp）になる。ゲノムサイズからカバレッジを計算できる。
seqkitで出力のmiseq_reads~を分析する。
seqkit stats miseq_reads_*
$ seqkit stats miseq_reads_*
file format type num_seqs sum_len min_len avg_len max_len
miseq_reads_R1.fastq FASTQ DNA 499,986 150,495,786 301 301 301
miseq_reads_R2.fastq FASTQ DNA 499,986 150,495,786 301 301 301
リード長はmax 301bp。
fastqcでクオリティを分析。
fastqc
R1
R2
手持ちのゲノムがない場合、NCBIからランダムにゲノムをダウンロードしてリードを発生させることもできる。ダウンロードにはNCBIのE-utilitiesが使われている。
ランダムに10ゲノムダウンロードしてリードを発生させる。
iss generate --ncbi bacteria -u 10 --model MiSeq --output ncbi_reads -p 8
time outした時は数を減らしてみてください。
abundance.txtファイル。分布にはモデルに従う。
defaultの分布はlognormal
iss model: generate an error model from a bam file
モデルを作る場合、手持ちのraw NGSリードをマッピングして、できたbamからモデルシーケンスプロファイルを構築する。詳細はGithubで確認してください。
InSilicoSeq/README.md at master · HadrienG/InSilicoSeq · GitHub
引用
Simulating Illumina Metagenomic Data with InSilicoSeq
Hadrien Gourlé, Oskar Karlsson-Lindsjö, Juliette Hayer, Erik Bongcam-Rudloff
Bioinformatics, bty630, 17 July 2018
