全ゲノム増幅したサンプルがうまくアセンブルできないときはCanuを試すと良いかもしれない

:::note info
2024-02-28 追記:
全ゲノム増幅を利用してエビに感染するウイルスのゲノム配列を解析した論文が出版されました。ぜひご覧ください。
[Kawato et al. (2023) "Evolutionary genomics of white spot syndrome virus"](https://doi.org/10.1007/s12562-023-01715-4)
:::
# 動機
[Multiple displacement amplification](https://en.wikipedia.org/wiki/Multiple_displacement_amplification)で全ゲノム増幅したサンプルからウイルスの環状ゲノムをアセンブルするとき、キメラリードのせいでうまくアセンブルできないことがあった。そのときのトラブルシューティングを記しておく。

今回使うのはエビに感染するウイルスのWSSV (white spot syndrome virus)。
[MAG: White spot syndrome virus E2 DNA, complete genome](https://www.ncbi.nlm.nih.gov/nuccore/AP027287)
![E2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2689564/a579bc5d-67f7-8b03-1855-b566ca31ff73.png)

[大型DNAウイルスゲノムの解析 (1) アセンブリ編](https://qiita.com/satoshi_kawato/items/c22a1ec469cb974ddad3)に沿って、事前に
- 極端に短いリード(10kb未満)の除去 
- WSSV参照ゲノム ([NC_003225.3](https://www.ncbi.nlh.gov/nuccore/NC_003225.3)) にリードを貼りつける
- 貼り付いたリードを回収 (```E2.up10k.mapped.fq.gz```)
を事前に行った。
# とりあえずFlyeでアセンブル
最もポピュラーなロングリードアセンブラの一つである[Flye](https://github.com/fenderglass/Flye)でアセンブルする。
```sh
conda activate flye-2.9
flye --nano-raw E2.up10k.mapped.fq.gz --out-dir flye_up10k  --threads 16
```
アセンブリグラフ(```assembly_graph.gfa```)を[Bandage](https://github.com/rrwick/Bandage)で可視化する。
![E2_flye_up10k_graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2689564/f1b5af47-3135-f003-2f20-35034ccb9d11.png)
アセンブリが２本のコンティグに断片化している。

```sh:flye_up10k/assembly_info.txt
#seq_name       length  cov.    circ.   repeat  mult.   alt_group       graph_path
contig_1        194659  90      N       N       1       *       *,1,*
contig_2        98037   101     N       N       1       *       *,2,*
```
カバレッジは90‐101Xと十分なようだ。
# 何が起きたのか
なぜカバレッジが十分にもかかわらずアセンブリが断片化してしまったのか。
WSSVの参照配列([NC_003225.3](https://www.ncbi.nlh.gov/nuccore/NC_003225.3))と小さいほうのコンティグを[pairwise BLASTN](https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome)で比較する。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2689564/cc08d64b-aa41-f83a-aaac-9fa59abdb86b.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2689564/56448640-0149-6745-7070-7081f709ae74.png)

ドットプロットを拡大する。
![hit_matrix.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2689564/82cd3f0b-bcb8-c2db-988e-4141f18284db.png)
当該コンティグは参照配列の140kbp-235kbpあたりにヒットした。つまり140kbpと235kbpの部位でアセンブリが分断されている。この領域に貼り付いたONTリードを[IGV](https://software.broadinstitute.org/software/igv/)で可視化すると、
```sh
  minimap2 --MD -Y -ax map-ont -t 8 CN01.fa E2.up10k.mapped.fq.gz | samtools view -Sb -  |samtools sort -@ $NUMTH -O BAM -o E2.up10k.bam
  samtools index E2.up10k.bam
```
![igv_snapshot2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2689564/707b1753-1790-6046-f122-2d63d14c91e2.png)
色付きの部分がキメラ。これがFlyeのアセンブリに悪さをしている可能性がある。
## Canuを試す
ロングリードアセンブラの[Canu](https://canu.readthedocs.io/en/latest/index.html)は塩基のエラー修正とトリミングを行ってくれる。トリミングのステップがMGAに起因するキメラリードをパージしてくれているようだ。
```sh
canu -d canu_up10k -p canu_up10k genomeSize=300k -nanopore-raw E2.up10k.mapped.fq.gz
```

```canu_up10k/canu_up10k.contigs.fasta
>tig00000002 len=302145 reads=862 class=contig suggestRepeat=no suggestBubble=no suggestCircular=yes trim=8655-297985
TTGTTGCTGCCCCGACAACACCCCTCGTACCTTCTAATGTGGAGGAAGAGGAAGAGGAAGAAGAGCAGATGGAAGAAGAGGAGGAAGAGGAAGTAGAAAG
GGAAGAAGGATCTGATAAGGAAGATGACGGAGACGCACCAGCACAGGAAGAAATGGAGGAGGAGAAGGAAGAAGAACAACAACAACAGCCAGAAGAAGAA
AGCAATGGTAATGAGAACCAAGAAGAAGAACAACAACAACAACAACAACCAGAAAGAGAAGAGGAGAATAAGGATGCAGATAGTGACA...
```
``suggestCircular=yes``ということで、無事環状コンティグをアセンブルできた。
余計かもしれないが、Canuが生成したエラー修正＆トリム済みリードをFlyeでアセンブルしてみる。
```sh
conda activate flye-2.9
### Flye using Canu-corrected reads ###
flye --nano-corr ./canu_up10k/canu_up10k.trimmedReads.fasta.gz --out-dir flye_up10k_canu_trimmedReads  --threads 16
```
無事環状コンティグが得られた。
```flye_up10k_canu_trimmedReads/assembly_info.txt
#seq_name       length  cov.    circ.   repeat  mult.   alt_group       graph_path
contig_1        289335  39      Y       N       1       *       1
```
![E2_flye_up10k_canu_trimmedReads_graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2689564/7ceba416-0626-a68f-670a-88581f426061.png)

エラー修正＆トリム済みリードを参照ゲノムに貼り付けると、
```sh
minimap2 --MD -Y -ax map-ont -t 20 ../CN01.fa canu_up10k/canu_up10k.trimmedReads.fasta.gz | samtools view -Sb -  |samtools sort -@ $NUMTH -O BAM -o E2.canu_up10k.trimmedReads.bam
samtools index E2.canu_up10k.trimmedReads.bam
```
キメラリードがほぼ消えている（下段）。
![igv_snapshot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2689564/4402aaff-e6b8-e383-3eee-10e486851c70.png)

ショートリードでポリッシュして、配列は完成。