2024-02-28 追記:
全ゲノム増幅を利用してエビに感染するウイルスのゲノム配列を解析した論文が出版されました。ぜひご覧ください。
Kawato et al. (2023) "Evolutionary genomics of white spot syndrome virus"
動機
Multiple displacement amplificationで全ゲノム増幅したサンプルからウイルスの環状ゲノムをアセンブルするとき、キメラリードのせいでうまくアセンブルできないことがあった。そのときのトラブルシューティングを記しておく。
今回使うのはエビに感染するウイルスのWSSV (white spot syndrome virus)。
MAG: White spot syndrome virus E2 DNA, complete genome
大型DNAウイルスゲノムの解析 (1) アセンブリ編に沿って、事前に
- 極端に短いリード(10kb未満)の除去
- WSSV参照ゲノム (NC_003225.3) にリードを貼りつける
- 貼り付いたリードを回収 (
E2.up10k.mapped.fq.gz)
を事前に行った。
とりあえずFlyeでアセンブル
最もポピュラーなロングリードアセンブラの一つであるFlyeでアセンブルする。
conda activate flye-2.9
flye --nano-raw E2.up10k.mapped.fq.gz --out-dir flye_up10k --threads 16
アセンブリグラフ(assembly_graph.gfa)をBandageで可視化する。
アセンブリが２本のコンティグに断片化している。
#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)と小さいほうのコンティグをpairwise BLASTNで比較する。
ドットプロットを拡大する。
当該コンティグは参照配列の140kbp-235kbpあたりにヒットした。つまり140kbpと235kbpの部位でアセンブリが分断されている。この領域に貼り付いたONTリードをIGVで可視化すると、
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
色付きの部分がキメラ。これがFlyeのアセンブリに悪さをしている可能性がある。
Canuを試す
ロングリードアセンブラのCanuは塩基のエラー修正とトリミングを行ってくれる。トリミングのステップがMGAに起因するキメラリードをパージしてくれているようだ。
canu -d canu_up10k -p canu_up10k genomeSize=300k -nanopore-raw E2.up10k.mapped.fq.gz
>tig00000002 len=302145 reads=862 class=contig suggestRepeat=no suggestBubble=no suggestCircular=yes trim=8655-297985
TTGTTGCTGCCCCGACAACACCCCTCGTACCTTCTAATGTGGAGGAAGAGGAAGAGGAAGAAGAGCAGATGGAAGAAGAGGAGGAAGAGGAAGTAGAAAG
GGAAGAAGGATCTGATAAGGAAGATGACGGAGACGCACCAGCACAGGAAGAAATGGAGGAGGAGAAGGAAGAAGAACAACAACAACAGCCAGAAGAAGAA
AGCAATGGTAATGAGAACCAAGAAGAAGAACAACAACAACAACAACAACCAGAAAGAGAAGAGGAGAATAAGGATGCAGATAGTGACA...
suggestCircular=yesということで、無事環状コンティグをアセンブルできた。
余計かもしれないが、Canuが生成したエラー修正＆トリム済みリードをFlyeでアセンブルしてみる。
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
無事環状コンティグが得られた。
#seq_name length cov. circ. repeat mult. alt_group graph_path
contig_1 289335 39 Y N 1 * 1
エラー修正＆トリム済みリードを参照ゲノムに貼り付けると、
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
キメラリードがほぼ消えている（下段）。
ショートリードでポリッシュして、配列は完成。