#標準出力修正
2021 7/11 condaインストール追記
2022/09/17 -hオプション修正
wgsimはfastqをシミュレートできるツールである。Wgsimは、SNPと挿入/欠失多型と二倍体ゲノムをシミュレートできる（シーケンスエラーを再現したり、diploidゲノムの多型を想定して、一定の確率で変異を入れることができる）。もともとSAMtoolsのパッケージでリリースされたが、その後スタンドアローンになった(ref.1)。
wgsimはARTなどのツールでは不可能な300bp以上の配列を発生させることでも可能である。ただしシーケンスプロファイル等はないので、シーケンスのbias等を再現するには向いていない点には注意が必要。論文にはまだなっていない（リンク）。
*wgsimは元々はSAMtoolsの一部
インストール
git clone https://github.com/lh3/wgsim.git
cd wgsim/
gcc -g -O2 -Wall -o wgsim wgsim.c -lz -lm
#conda
mamba install -c bioconda wgsim
> wgsim
ラン
ref.faの配列を使い、デフォルトの条件でfastqを発生させる。
wgsim ref.fasta R1.fastq R2.fastq > mutation
- -e base error rate [0.020]
- -d outer distance between the two ends [500]
- -s standard deviation [50]
- -N number of read pairs [1000000]
- -1 length of the first read [70]
- -2 length of the second read [70]
- -r rate of mutations [0.0010]
- -R fraction of indels [0.15]
- -X probability an indel is extended [0.30]
- -S seed for random generator [-1]
- -A disgard if the fraction of ambiguous bases higher than FLOAT [0.05]
- -h haplotype mode
カレントディレクトリにR1.fastqとR2.fastqができる。デフォルトだと70-bpでインサートサイズが500、 SDが50、リード数が1000000などとなっている。また、上記のデフォルトに設定された確率でindelやSNV、シーケンスエラーが導入されている。カバレッジはリード数xリード長をゲノムサイズで割って計算できる（ペアエンドはx２）。
カバレッジにはある程度ばらつきができる。
リード長150-bpのfastqを30000リード数発生させる。
wgsim -1 150 -2 150 -N30000 ref.fasta R1.fastq R2.fastq > mutation
1-kbp以上のfastqを発生させることも可能であるが、シミュレーションには相当な時間がかかる。ロングリードを作るなら、pacbioかnanoporeのシミュレーターを使うことをお勧めします。
2021 7/12
indelもSNVもシーケンスエラーもゼロのfastqをシミュレートする。
wgsim -e0 -r0 -R0 -X0 ref.fasta R1.fastq R2.fastq
k-mer profile をチェック
たしかにユニークk-merがほぼない（エラーや変異がなくても、カバレッジの薄い部分でわずかにユニークk-merができている）。
-hをつけると2倍体ゲノムのシミュレートになる。1倍体ゲノムのシミュレートになる。
変異の頻度ですが、ヒトゲノムだとヘテロ接合の部位の平均存在率は0.001という数値がアフリカ系人種から出ているそうです（-rのデフォルトは0.001）。
https://www.biostars.org/p/6177/
引用
https://www.biostars.org/p/137165/
https://popmodels.cancercontrol.cancer.gov/gsr/packages/wgsim/
ref.1
https://bmcresnotes.biomedcentral.com/articles/10.1186/1756-0500-7-40
