macでインフォマティクス

macでインフォマティクス

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

メタゲノムのシミュレーター InSilicoSeq

2019 5/13 コメントおよびパラメータ追記

 

 ますます多くのバイオインフォマティクスツールがリリースされており、特定の実験に最適なツールや最適なツールを知ることは困難になっている。ゲノミクスとメタゲノミクスのデータのシミュレーションは、実験の計画と新しい方法の開発の両方において重要な役割を担っている。実際のデータとは対照的に、シミュレートされたデータは、例えばメタゲノミクスの場合には、サンプル中に存在する種の存在量など、制御されたパラメータで生成することができる。パラメータを固定することで、制御された条件でのツールのベンチマークとテストが可能になり、新しいツールやパイプラインをテストするための模擬データも提供されます。このようなテスト環境は、メタゲノミクス(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/

InSilicoSeqに関するツイート。

 

インストール

mac os10.12のpython2.7.10でテストした。

本体 Github

pipで導入できる。dockerイメージも用意されている。

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)

——

 

 ラン

2つのモードがある。 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にまとめられる。

f:id:kazumaxneo:20180720125422j:plain

存在量が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

f:id:kazumaxneo:20180720142312j:plain

R2

f:id:kazumaxneo:20180720142237j:plain

 

 手持ちのゲノムがない場合、NCBIからランダムにゲノムをダウンロードしてリードを発生させることもできる。ダウンロードにはNCBIE-utilitiesが使われている。

ランダムに10ゲノムダウンロードしてリードを発生させる。

iss generate --ncbi bacteria -u 10 --model MiSeq --output ncbi_reads -p 8

time outした時は数を減らしてみてください。

abundance.txtファイル。分布にはモデルに従う。

f:id:kazumaxneo:20181014122145p:plain

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