NGSシミュレーションの古典的ソフトwgsimで遊ぶ 〜謎の生物プー〜

## はじめに

謎の生物 **プー** それは、**太古から存在するかもしれない 謎の生物**です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/22ec5a45-50d6-e446-324e-7303617d3568.png)
プーの想像図


## wgsimとはなんじゃらほい？

wgsimとはバイオインフォ界の巨星 Heng Liさんが作ったNGSのシミュレーターです。
2011年のソフトですが問題なく動きます。

https://github.com/lh3/wgsim

## インストール

Biocondaを導入している場合は

```
conda install wgsim
```

## 〜 謎の生物プー 〜

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/22ec5a45-50d6-e446-324e-7303617d3568.png)

謎の生物プーは、まったくの気まぐれから生み出されました。
神は生物を創造します。（進化など言語道断です。）
あなたは神なので、ちょっとしたウイルスと同じぐらいの配列、5000塩基の生物を創造しましょう。この生物の名前を「プー」と名付けました。

プーのリファレンス配列は、ランダムで生成しています。（筆者は絶滅危惧種の**Ruby派**なのですが、ほかのプログラミング言語でも同じことができるでしょう。）

```ruby
genome = Array.new(5000){ %w[A C T G].sample }
puts '>poo chr1'
genome.each_slice(80) do |l|
  puts l.join('')
end
```

このスクリプトを実行すると、こんな感じで配列が出力されると思います。

**これが謎の生物**プー**のリファレンスゲノムだッ！！！**

```
>poo chr1
AAAGCGCTGGATCATGCTTTAAGGCTTAGTTTTATGGCCCTAAACACGGAGACTGAAGCCGCTTGAGTCCGGTTTCAGAG
AGATTTGGTTTCACCCACAACTGCCAGCATATAATGGACCGTCGGCCGTACACGTTATGAGCCTGGCAGGGGTTACTACG
CGCCACGGGCTAGAGACCTAGTGGCGAAACGTGAAGGGCAATTATGGAATCCCGTACGCGCAATAAATGTCTACGCATTC
CCGGGGGGTGTTGTAGTCCCACGGCCCCGGACCTGCATTGGTCCAAATTCTCTTCGGTGCATATTAGCTTCATCGGAACC
GGCCTTTCTAATTCAGAACTGTTCGGTCATACACTATCTTCCAACACGCCGTGCACGAGTTGCACGGAGATCTAGTTGTT
ACTCGACGCGCGTACGACCGAGGCGGGGCTTAACACGCGAGTAGAAGCCCTGGGTTATCCACCACTTAATTTTGTAAAAA
GCGCTTTCAGTACTACACGGCTTGACCGTAGCTCTGTCGGACTTGTGACCCCTTCCGGCGAACATTAATTGTACCCCTTC
TTCTGACGGATTGGTTAGTAGTCCAAGCAGTCTTCACCCCAGGAGCACTTCTTCTTGCCGCGAGGCCATTCAGTCCAGTT
CGTGATGGCGGGGGTACCGATGCAAAAGCTCGTCTTTCTTGAGGCATCCACAGCCAGAGGCGGGGTATTTACGAGATTAA
GTTCTATCATTGCACTAGACAAAGCAAGTACAGCTCGGCAAAAAGAAACTCGTGGGCGCGTCTTGACCTGCAGCTAATAC
TGACTAGCATCCGTAGATATTACCGATCGTGCCGCACAAGGGAAATCTCGAGGCTGACAACGTCTTTGCCAAGGCAGTGG
CCGTGTCTTCTCGCCAGGGGTAAAGGTGCCGCTGGACTAGCGCTAACTTAGTGTCGTAATATCTCCATATGGACGATTTC
CCCGATCTTACCGGCTAGCCGCTGGTTGTATGCCGCGCCCAGGCCACCGAGTCTGACATGTTCCAAGGGCCTCATGTGAA
AGCGATGTGCTGATCCAGCATGGCCTTCAGCCTCACAAGACTGAGAGTTGGTGTAGAAACCCGGTTTGTTGGTAGTATGA
TGCTTAAAAATCGGTTTAATTTAGCTTGGTTAGATTTGAACGTGAAGAGCCTATGTCCAAAGACCAACTAACCTACTCAG
ACTAAGATTTACCGATTTCTCGTTTTGACTGAGATACGTAGTTAGCTATACGGGGCGTCGCCATCCTCTACCCCTTGAAT
ATGATTCGTAAGTGATGCCATTCTCGTGCTCGTGCGCTATGAATAGCAATGTCACCACAATAGTGCCCGGTGGCCGACAA
TCCGAGACCATGTCCGTCCTCATTCATGCGTAGTTAGTAGTGTATCGACTCACCTGGAGATGACCTCTATGTCCTGCTGC
GGTACGGCACGGCAGGATCTTTCTAGAAACATCGGCGATGCCCCAGAGGTTCTCTAAACTGTGTCATCTTTTTCTGCGGT
GCGTGGCTGGTCAAACCATTCAAATATATGAGAGTCATGCGCTAGCTCCCGGTAAAACCCCACAATAACGGGATTTGGGG
TTGAAAATGCGAGGCAGCGCACATGGATCAGGCCATTTAGGGTGTGATTGGTATGTTCCGAAAAAGCTGATGAGAACTAT
GCCCCGTTTACTTTGCCAGTAGAAAGAGTAGGAAGTAGTGTAACGCGGGCACCGGAAAACGACCGACACGATACTTTCCC
TTGTGCTGGAATTGTTTGATATATTAGGACGAGAATAGTTCTCTCCATTCCGCCGGCGTATAGACCAGGCCATGCCGTTC
TATTTCCAAGTTCCTTTGAAGGATCGGTGTTATTTCGGATGTCGTCGAGTAGTCGTGGGGCGTTAGGGTCTAAGCGACCA
GTTGTAGCTAATCGGGTGTCCTCCACTAACTGCTATTAGCATAGTACACCATGCTGCTTCAACTCATTTGGGATTTGCAC
TACTCTAGGCCAGAGATACAGATGCAAGCCGACTTAGCTGCACTAGCCAAGATTGAGCCCTAGGGCGTAGCGACCCTGAG
AATCCTGGCTAGCCCCCCTAGCCCCGGAAACTGCCCGTGTTGCTCTCGAAGCAGTTAGTAATCATCGTAATTTTGTCAAG
ATGGCCGTAGAGTGCAACCTGCAAGGCTGCAGGGTGTGTTCGGCATATTAATGCGGACAGTGCAGTGGACAAAGAGTCTC
GGAAGCGACGTGTATACAGTAGTGAGTAAATTGTTAAGCTTACTATCGCAGGAACTTCTGATCTGCCCGACCGTCGAGAT
CGAGTCTAGTAACAAGATTGCAGTCCTCAAGACTAGAAGTACCTAACTCCCTCAAATGCGAGTCCCAACAACTCATCCTC
ACCTCTGGAATCTCACCTCTAAGCATGACAGGCTAGGACCATTTAACATAGTTACAGTGACCTGATACGTAACGCCTCGA
AAGGAAAGGCGCGCCTAGGCAGCAATTGCAACGAACTCGGACGTCAGTACTCAACCAACCGAGTTATTTCAGCGTCGCGA
ATCTAATGCATCTTCCTGCTTAAGACGTCGTACTCCGGTGTCACCCTTCGACTAACTTAAAATCGCGCATTGGCAAGCAT
ACCTCTAATACCGGATCGCTCCCCTCAGGTTTACTTGATGAGCGCTGAACCCAACCATAATTTCGGCAAGGTAATCGCTT
ATCGTTCCACCCAGGCGTGAATAATTCATGCCCAATGCAGGTGACCGTATCTCAGGGAGATATAAGTGAAAGTTATGTGG
AAGGTATTGTAGTACATTTCAGGGTTCATCTAAGCTTCAGGACCACAGATGCAGAAGACGGGTTAAAGGACATGTCGGCA
CGAAGACATATCTTGTGGATCCCTATTCAGTGCCCGTGACTCAACAGATAGACTCATTTTCTCCTTTGGCGTGTCCGCAG
GTAAGATAACTTCAGTGTCCGAGCGGCGGAGTGGGTGTTCACGTCGTCGAGTGTCCATTTCGATTGGGGATGTTGGCTCG
CATACCACACGACCCTTTAGTGGACGGTCTAGTTTGTATACTCCCCGATTGGACTTTACCCGCGAAGCTAGCACTGTAAA
TATTGCCCCGGAGTGGCGGGATACGTTGGATATGCCTATACTGGACCGCAGCCTATTACATCAGATCGGCTTACCGCCAA
TCTTAGTGCTCAGCTACCCAGAGCTGTATCCATGGACATAGGCGCTATGACAAGCAGTACTGTAAAGCAGGGACAGCTTT
ACTATGTCGACAGCCGGCGTACGTCGTTGGGCTACCGCCGGCCTAGTTCCTAGTCAGTTCCGAGGCATAATCTACAAGTC
TAAGTAGGGCGCCATGATGTCACGAGACTTATAGGAGTCCACCGGCAAGTGCGTTGACGTAAATACATATGTGTATGCAG
CACTTTTGACAGGCCCGTAGATTTTGCACAAACAGCCGGGTGTGGTTCAAAAAGGGTGTTGCGGTCAGCCGCAGCAGCGC
AATTGTAACCAATCTACCTCCCGTTTGCCGACCGCTCGCGCAGATTCATGTAGGAACCGCTAGTTCCCCGTCCGCTACAT
TTCAATGGGCTTTTGCATGTCCTATCGCGATACCATCTCTGGTGGCGATGATATTTTATCAGGTTGCATTGGCCCAGTCC
ATCTTCGCGGCCGATTGGGCACGCCATTAGTTGCTCTCTGCGCAATCGCAGCTATACTAATGGTATGGTTGGGGAAATCA
TTGGAGTGCATTATTCCGCCCTCTTAACCTGTTGTGAACGGAACCGCAGTACGGTAGATTGCAAAGTCCAAGGTTTCTGG
ATTGTCTCTGAAAATGAAAGTTCTGCTCTGTCAAGTGATAAACGGAACGTCCTAGGAGGCGCCTACCGCGTTTAGAGAAG
TAACCCACGATATCTCGAACAAAAAATCAAAAGCGGTGTAGGTACAAAACGTATCGCTTACATCCCGTGAAGCTTAGTAG
GACGCACCCTTCATCTCTCAGGGGGTCGCGTACACGATGTCCGTCTTTTATCGTGCCACGTGTACAAGTAGGGGTGCAGA
AATAACTGTACTAATGCAACGTTTCAGTTCTTCTGATGAACATCTCACACGGTTGTCCCCTCTAATATTAGGATTATCGG
CCTAGACGTGTTATGCACGCCGACCCCTTTCCCTAAGTGCTTCCGCAATGAAGTCGCGTATTTGTAACACACATTGGGTA
TGAATATAACTGTCCCTCCTGGAATAAGGAGTGCGCTTATCCTTTCGGACGTCAGAGCCCAGGGAGGAGACGCGTCAGAT
GCCATAGCCTATCGGGATCCGTCGTGCTAGTAGAACCAAAGTAACGGTGCTGTTAAATCGTTCCCTCAGCAGTCCGGGGT
TGACAGGTACAGTGTTTCCAACTACTAGAGGATCTCGCTGCTATAGAACGATACCACTGCGGGGGACACCGCCACTTTTC
GAGCTGTGGGGCCAGAAGTGCTGTCTGGATAATTCGGTAGGCATCAGCGGACAGCGGAGCTCGTAAGTGCCATGCACTTG
TAGCCCATTACGATGAGATCGGTTCCGCCAACTCCGTCAGCACAGGTATCCACGAGGAGTCACCGTCCGGCATGACGCAG
GACGGGGGCGGTGGAAATAATCGAGGCGGATACTCTTCATTTGTAGGCGCCAGCTACAGCGGTCCGGCGGGTCAATACCG
CATCCGCAGACATAATTGACCTCACTCAGCTCGCACGAACCTCCCTACTCTTTTCTACTCGAGGGTGCCTGATCAGAGAA
CGGAAGACGCTATAGGCAAGGCGGGCACGCCTGGAATACAGTCCCTTTCCGAATATTCACACCCGCTACCTATGTCGTAG
CGTTGCCGGGCCCGTTATCCAGACCGAACCCGTTGTCTGTTATAACATATTGAGTAATGGGCAGACAAACCCGTGTAGAG
TCCTCTCGGTTCCACATATCGAAGAGGGCGATCCGTACGA
```

ながい。

`poo.fa` として保存します。

## wgsimを動かす 〜プーの全ゲノムシークエンス〜

　さて、謎の生物**プー**は　ヒツジ用の歯磨き粉から発見されました。ところでヒツジってハミガキするんでしょうかね？さあ、私は詳しいことはわかりません。ただ、プーは若いヒツジがこっそり使用していた歯磨き粉から発見されたと報告されてますので、それを信じることにしましょう。学名は、顕微鏡でみた様子がぷるぷるしていることから、プル・プルプルスと言うそうです。（今適当につけた）

　さてプーを発見した人はその後人生に山あり谷ありプーに対する興味を失って、その後50年間プーは放置していました。しかし天と地の気のめぐり、赤道の気温の上昇、町内会長のおののき、ほのかな春のおとづれ、会長と社長の和解など、たまたまいろいろな偶然の成り行きがあって、プーが全ゲノムシークエンスされる時がやってきました。さて、あなたは神なのでゲノムシークエンスの結果をwgsimで求めます。

wgsimの実行に先立って、オプションをみてみます。

```
wgsim --help
```

```
Options: -e FLOAT      base error rate [0.020]
         -d INT        outer distance between the two ends [500]
         -s INT        standard deviation [50]
         -N INT        number of read pairs [1000000]
         -1 INT        length of the first read [70]
         -2 INT        length of the second read [70]
         -r FLOAT      rate of mutations [0.0010]
         -R FLOAT      fraction of indels [0.15]
         -X FLOAT      probability an indel is extended [0.30]
         -S INT        seed for random generator [0, use the current time]
         -A FLOAT      discard if the fraction of ambiguous bases higher than FLOAT [0.05]
         -h            haplotype mode
```

ところどころ私にはよくわからないオプションがありますが（-s が何の標準偏差なのか、-A とか -h がどう動くのかとか）、それは後で勉強することにしましょう。（しない）

ここではリード数を1000本にします。

```
wgsim -N 1000 poo.fa poo_1.fq poo_2.fq
```

実行すると突然変異が表示されます。`poo_1.fq` と `poo_2.fq` が生成されたと思います。簡単ですね。

```
poo	977	A	W	+
poo	1210	T	W	+
poo	1655	G	K	+
poo	1874	T	K	+
poo	2043	C	S	+
poo	2378	G	T	-
```

生成されたfastaファイルの中味を見てみましょう。

```
head -n 12 poo_1.fq
```

```
@poo_418_960_2:0:0_1:0:0_0/1
TCGTCGAGTAGTCGTGGGGCGTTAGGGTCTAAGCGACCAGTTGTAGCTAATCGGGTGTCCTCCACTAACT
+
2222222222222222222222222222222222222222222222222222222222222222222222
@poo_624_1182_3:0:0_1:0:0_1/1
CTTAGGACATAGGCTCTTCACGTTCAAATCTAACCAAGCTAAATTAAACCGATTTTTAAGCATCATACTA
+
2222222222222222222222222222222222222222222222222222222222222222222222
@poo_2608_3166_3:0:0_0:0:0_2/1
GTCCAGTATAGGCATATCCAACGTATCCCGCCACTCCGGGGCAATATTTACAGTGCTAGCTTCGCGGGTA
+
2222222222222222222222222222222222222222222222222222222222222222222222
```

まあ、普通リードの後半の方がクオリティスコアが落ちてきたりすると思いますが、そうはなっていませんね。
Fastqcでsequence qualityを見て確認してみましょう。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/5774428c-af4c-72b4-c74e-8b6d3a2a9a02.png)

クオリティが一定（しかも低い）ことがわかります。

一般的には下図のようにリードの後方にいくにつれて、クオリティが下がっていくはずです。（画像は Fastqcの公式ホームページから　https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/04abb0e4-2056-f35b-2188-5adc62223db3.png" width="50%" height="50%">


まあそこのへんは昔のソフトなので仕方がないということで。

## BWAを使ってマッピング

さて、全ゲノムシークエンスなので、bwaでマッピングします。なになに。minimap だと? ここは10年前の世界観で、wgsimで遊ぶ記事なんだぞ、ロングリードなんて存在しないんだ！！おっと取り乱してしまいました。神は取り乱さないように気をつけなければなりません。

まずはともあれindexの作成します。

```
bwa index poo.fa
```

プーのゲノムサイズは小さいので一瞬で実行が完了すると思います。
それではbwaをかけてみましょう。（bwa も Heng Li氏の作ったツールです。まさに神といって良いでしょう。）

```
bwa mem poo.fa poo_1.fq poo_2.fq > poo_wgs.sam
```

さて、SAM → BAM、BAM → sorted BAM、インデックスまでつけてやりましょう。


```
samtools view -bS poo_wgs.sam > poo_wgs.bam
samtools sort     poo_wgs.bam > poo_wgs.sorted.bam
samtools index    poo_wgs.sorted.bam 
```

## IGVでリードの様子を見てみよう

IGVで、`poo.fa` をリファレンスゲノムに設定し、`poo_was.sorted.bam` を開きましょう。

うおぉぉぉ！できたよ！できたよ兄さん！（うるさいな）

![igv_snapshot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/8c717409-dbae-3f2d-6894-e3e85278773c.png)

さてIGVには、実は便利なスプリットビューの機能がついています。これを利用して、突然変異の場所を観察してみましょう。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/41c2a358-0a3d-6bb1-86e0-72e9dcdfe79d.png)

うおぉぉぉ！兄さん（略）

ちゃんと突然変異が検出されていますね。

## バリアントの検出

人間の検体を使う場合はGATKとかで頑張るのかもしれないが、pooは適当な生物なので、bcftoolsとかでコールしておけばOK。

```
bcftools mpileup -f poo.fa poo_wgs.sorted.bam | bcftools call -mv -Ob -o poo.bcf
```

などとして、

```
bcftools view poo.bcf
```

とすればよし。

## おわりに

さて、謎の生物 **プー** に関する記事を気まぐれで書いてみましたが、思いのほか書いていて楽しいのと、勉強にもなることが判明しました。プーシリーズ、また気が向いたら書くかも。