2024-02-28 追記:
この記事で紹介した手法を利用して、エビに感染するウイルスのゲノム配列を公共NGSデータから収集、解析した論文が出版されました。ぜひご覧ください。
Kawato et al. (2023) "Evolutionary genomics of white spot syndrome virus"
多くのDNAウイルスゲノムはイントロンを持たない。遺伝子が密に配置するため、転写産物は互いにオーバーラップすることがある。また、ウイルス感染したサンプルからのRNA抽出時に (細胞由来DNAに比べて大過剰に存在する) ウイルスDNAを除去しきれず、ウイルスDNAに由来するリードが一定量発生する可能性がある。転写産物のオーバーラップとゲノムDNAの持ち越しのどちらの貢献度が大きいのかは不明だが、RNA-seqデータからDNAウイルスのドラフトゲノム配列を取得できることがあることに気づいた。
今回使うのは、ブタに感染する大型DNAウイルスであるアフリカ豚熱ウイルス (ASFV)に感染した培養細胞のRNA-seqデータ。
目次
リードの下処理
シーケンサーから出力された元データ (生リード raw reads) はアダプター配列や低品質の領域を含んでいる可能性があるため、それらの除去 (クオリティトリム quality trimming) を行う。今回はFastpを使う。この作業を行わなくても解析結果に問題は生じないかもしれないが、習慣づけておいて損はない。
# Download raw reads
nohup wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR111/052/SRR11185052/SRR11185052_1.fastq.gz &
nohup wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR111/052/SRR11185052/SRR11185052_2.fastq.gz &
wait
fastp -i SRR11185052_1.fastq.gz -I SRR11185052_2.fastq.gz -o SRR11185052_1_trim.fq.gz -O SRR11185052_2_trim.fq.gz --html SRR11185052_fastp.html --thread 8 --trim_poly_x
参照ゲノムへの貼り付け
トランスクリプトームデータの大部分は宿主由来であるため、それらを取り除くことで以後の計算量を節約できる。すでにウイルスの参照ゲノムが利用可能な場合、参照ゲノムに貼りついた (マップされた mapped) リードのみを利用するとよい。
ここで注意しなければならないのは、ウイルス由来でありながら参照ゲノム上に存在しない配列が除去されてしまう可能性があること。ウイルスの遺伝子組成に多様性が存在すると予想され、かつ宿主ゲノムのゲノム配列が利用可能な場合、宿主ゲノムにリードをマップして、マップされなかったリードを抽出するとよいだろう。ただ、宿主ゲノムはウイルスゲノムより格段に大きいため、より計算に時間がかかってしまう。
ASFVの参照ゲノム配列にリードをMinimap2でマッピングし、貼りついたリードをSAMtools fastqで抽出する。
(minimap2 --MD -ax sr -t 8 ASFV.fa SRR11185052_1_trim.fq.gz SRR11185052_2_trim.fq.gz | samtools view -@ 8 -G 12 -bS - |samtools sort -@ $NUMTH -O BAM -o SRR11185052.bam )2>SRR11185052.minimap2.log
samtools index SRR11185052.bam
samtools collate -u -O SRR11185052.bam |samtools fastq -@ 8 -1 SRR11185052_mapped_1.fq.gz -2 SRR11185052_mapped_2.fq.gz -0 /dev/null -s /dev/null -
正規化
遺伝子間でその発現量は数十～数万倍異なることがある。これはRNA-seqの読み取り深度 (sequencing depth) が遺伝子間で劇的に異なることを意味する。
これは元のゲノム配列を再構築するうえで少々問題である。たとえば、ゲノム上に読み取り深度数十× (かけ) と数万×の領域が混在している場合、大深度の領域のデータ量が膨れ上がるだけでなく、アセンブリに支障をきたす可能性がある。
シーケンス深度の差異を均す (正規化 normalize する) ことで、データを圧縮するとともに、 願わくばミスアセンブリを減らす。今回はBBTools内のBBNormを用いる。
# bbnorm coverage normalization with increasing MINDEPTH cutoffs
for MINDEPTH in 1 5 10 20 ;do
bbnorm.sh \
in=SRR11185052_mapped_1.fq.gz in2=SRR11185052_mapped_2.fq.gz \
out=SRR11185052_mapped_bbnorm_mindepth${MINDEPTH}_1.fq.gz \
out2=SRR11185052_mapped_bbnorm_mindepth${MINDEPTH}_2.fq.gz \
target=200 mindepth=${MINDEPTH} &>SRR11185052_mapped_bbnorm_mindepth${MINDEPTH}.log
done
最小深度 (minimum depth) の閾値を複数振るとよい。筆者は1, 5, 10, 20と振っている。
最小深度に閾値を設けることで、シーケンシングのエラーに起因する雑多な配列を捨てられる。しかし、あまりに閾値を厳しくすると、低発現のウイルス遺伝子由来のリードまで捨ててしまうため注意が必要だ。
正規化によってデータ量を大幅に圧縮できた。
$ seqkit stats *.gz
file format type num_seqs sum_len min_len avg_len max_len
SRR11185052_mapped_1.fq.gz FASTQ DNA 11,898,694 1,783,533,508 15 149.9 150
SRR11185052_mapped_2.fq.gz FASTQ DNA 11,898,694 1,783,210,694 15 149.9 150
SRR11185052_mapped_bbnorm_mindepth10_1.fq.gz FASTQ DNA 200,897 30,110,803 34 149.9 150
SRR11185052_mapped_bbnorm_mindepth10_2.fq.gz FASTQ DNA 200,897 30,082,732 27 149.7 150
SRR11185052_mapped_bbnorm_mindepth1_1.fq.gz FASTQ DNA 318,255 47,697,754 28 149.9 150
SRR11185052_mapped_bbnorm_mindepth1_2.fq.gz FASTQ DNA 318,255 47,637,806 22 149.7 150
SRR11185052_mapped_bbnorm_mindepth20_1.fq.gz FASTQ DNA 191,322 28,676,086 34 149.9 150
SRR11185052_mapped_bbnorm_mindepth20_2.fq.gz FASTQ DNA 191,322 28,651,519 27 149.8 150
SRR11185052_mapped_bbnorm_mindepth5_1.fq.gz FASTQ DNA 218,618 32,766,810 34 149.9 150
SRR11185052_mapped_bbnorm_mindepth5_2.fq.gz FASTQ DNA 218,618 32,733,379 27 149.7 150
配列の再構築
読み取られた断片 (リード) からもとのウイルスゲノム配列を再構築する (アセンブリ assembly)。
参照ゲノムに貼りつけたリードをアセンブルせずそのまま変異検出 (variant calling) すればよい、という考えもある。しかし、大規模な欠失や転位といった構造変異を検出するには、やはりde novo アセンブリが有効だ。
de novoアセンブリ
リードをもとの配列に再構築する作業をde novoアセンブリという。今回は代表的なショートリードアセンブラのSPAdesを用いる。
for MINDEPTH in 1 5 10 20 ;do
R1=SRR11185052_mapped_bbnorm_mindepth${MINDEPTH}_1.fq.gz
R2=SRR11185052_mapped_bbnorm_mindepth${MINDEPTH}_2.fq.gz
# SPAdes assembly
conda activate spades-3.15.5
spades.py --isolate \
-o SRR11185052_spades_k127_bbnorm_mindepth${MINDEPTH} \
-1 $R1 \
-2 $R2 \
-t 20 \
-m 50 \
-k 23,35,57,79,91,103,115,127
conda deactivate
done
スキャフォールディング
アセンブルして得られた配列 (コンティグ contigs) はいまだ断片的である。そこでコンティグを正しい順番につなぎ合わせる (スキャフォールディング scaffolding) ことで、もとのウイルス配列を復元する。今回はBLASTNで参照ゲノムに対して相同性検索を行い、Chromosomerを使ってコンティグを参照ゲノム配列に照合し、正しい順番に並び替える。
for MINDEPTH in 1 5 10 20 ;do
cd SRR11185052_spades_k127_bbnorm_mindepth${MINDEPTH}
blastn -subject ASFV.fasta -query scaffolds.fasta -outfmt 6 > scaffolds.blastn.out
awk '{if($4 >= 500 && $11 <= 1e-20 ){print $0}}' scaffolds.blastn.out >scaffolds.blastn.mod
chromosomer fastalength scaffolds.fasta scaffolds.length
chromosomer fragmentmap --ratio_threshold 0.8 scaffolds.blastn.mod 100 scaffolds.length scaffolds.map &>chromosomer.fragmentmap.log
chromosomer assemble scaffolds.map scaffolds.fasta scaffolds.chromosomer.fa
cd ../
done
特に深い意味はないが以後mindepth=1の結果を使う。
cd SRR11185052_spades_k127_bbnorm_mindepth1
seqkit stats scaffolds.chromosomer.fa
file format type num_seqs sum_len min_len avg_len max_len
scaffolds.chromosomer.fa FASTA DNA 1 176,679 176,679 176,679 176,679
cd ../
参照ゲノム上に存在しない配列が独立したコンティグになってしまった場合、参照ゲノムを用いたスキャフォールディングでは、そのコンティグを取りこぼしてしまう可能性がある。スキャフォールドに組み込まれなかったコンティグにウイルス遺伝子らしき配列が含まれていないか注意したい。
配列修正
スキャフォールディングして得られた配列 (スキャフォールド scaffold) は配列不明の領域 (ギャップ gap) を含んでいる。この領域を可能な限りなくすため、スキャフォールドにリードを貼り付け、その情報をもとに当該領域のギャップを埋める (ギャップフィル gap filling)。
今回はギャップフィルだけでなく、挿入・欠失 (インデル indels) や塩基置換など他の配列エラーの修正も同時に行う。アセンブリ上のエラーを修正する作業を一般にポリッシング (polishing) という。ポリッシングを行うソフト (ポリッシャー polisher) にはMaSuRCAパッケージのPOLCAやPilon などいろいろある。今回はギャップフィルに適したPilonを選択する。
変化がなくなるまでポリッシングを繰り返す。
cp SRR11185052_spades_k127_bbnorm_mindepth1/scaffolds.chromosomer.fa pilon0.fa
for NUM in 0 1 2 3 4 5 6 7 8 9;do
(minimap2 --MD -ax sr -t 8 pilon${NUM}.fa SRR11185052_mapped_bbnorm_mindepth1_1.fq.gz SRR11185052_mapped_bbnorm_mindepth1_1.fq.gz | samtools view -@ $NUMTH -F 4 -bS - |samtools sort -@ $NUMTH -O BAM -o pilon${NUM}.bam )2>pilon${NUM}.minimap2.log
samtools index pilon${NUM}.bam
conda activate pilon-1.24
java -jar /home/kawato/apps/pilon/pilon-1.24.jar --output pilon$((NUM + 1 )) --changes --vcf \
--genome pilon${NUM}.fa \
--frags pilon${NUM}.bam \
--mindepth 5 --minmq 20 --mingap 1 \
>pilon$((NUM + 1 )).out 2>pilon$((NUM + 1 )).err
conda deactivate
cat pilon$((NUM + 1 )).fasta >pilon$((NUM + 1 )).fa
done
file format type num_seqs sum_len min_len avg_len max_len
pilon10.fa FASTA DNA 1 176,897 176,897 176,897 176,897
参照ゲノム配列 (NC_001659.2)とのペアワイズBLASTNの結果は以下の通り。
BLASTN検索したところ、MN393476.1という配列と100％一致することが分かった。MN393476.1の長さは190576 bpで、今回アセンブルした配列より１割弱長い。
参照ゲノム配列とMN393476.1のペアワイズBLASTNの結果は以下の通り。
どうやらMN393476.1に存在し、参照ゲノム配列から欠落している配列があるらしい。今回アセンブルした配列からは、その部分がごっそり抜けて落ちているのだろう。マッピングに使う参照ゲノム配列をMN393476.1に変えると、また違った結果が得られると考えられる。
ともあれ、RNA-seqデータからまずますなASFVゲノム配列を得られたといってよいだろう。
データセットの質に依存するが、とある別の大型DNAウイルスでも良好な結果を得ている。
考察と注意点
今回紹介した手法は、重度に感染したサンプルにのみ適用可能と考えられる。RNA-seqから構築したドラフトゲノムは、全く発現していない遺伝子や長大な遺伝子間領域、そしてウイルス由来でありながら参照ゲノム上に存在しない配列が欠けてしまう可能性があることに注意が必要だ。
同じウイルス株から得たRNA-seqデータが複数利用可能な場合、それらを合わせて一つのデータとして解析すると、ドラフトゲノムの完全性、連続性を改善できる。
単一サンプル中に複数株が混合感染している場合、得られるドラフトゲノム配列はそれらのコンセンサスになると予想される。どれか一つの株がドミナントな場合、その株のゲノム配列を得られると期待される。しかし、半々の場合はキメラ状になってしまうかもしれない。IGVでマップ結果のBAMファイルを可視化し、配列組成が雑多な領域がほとんどないことを確認する必要がある。
正規化後のカバレッジが十分 (100×以上) あるにも関わらず、アセンブリして得られたコンティグがギャップだらけということが時折ある。原因はおそらくサンプルごとに異なり、よくわからない。最小深度の閾値を変更することで改善する場合がある。
ヘルペスウイルスなど一部ウイルスの遺伝子はイントロンを有しているため、今回の方法で全ゲノムをアセンブルすることはできないと予想される。また、遺伝子配列と異なる塩基配列を生み出すRNA editingの影響を受けている可能性も考慮が必要と考えられる。
以前の記事で、NCBI SRA上のデータからウイルス配列を含むデータセットを検索する方法を紹介した。本記事では、いわばその応用として、RNA-seqデータからDNAウイルスのゲノム配列を取得する方法を示した。ウイルスの比較ゲノム解析を行う際、ゲノム情報が必ずしもふんだんに整備されていないことがある。利用可能なゲノム配列の数を少しでも増やすために、RNA-seqデータから再構築したドラフトゲノムが役立つかもしれない。