BigQueryを使ってNCBI SRAの配列データからウイルスゲノムを釣り上げよう

2024-02-28 追記:
この記事で紹介した手法を利用して、エビに感染するウイルスのゲノム配列を公共NGSデータから収集、解析した論文が出版されました。ぜひご覧ください。
Kawato et al. (2023) "Evolutionary genomics of white spot syndrome virus"

背景

次世代シーケンサーから出力された配列データには、実験過程程で発生した汚染 (コンタミ) や自然感染に由来するウイルス配列が紛れ込んでいることがある。思いがけないソースから見つかったウイルス配列は、新たなウイルス遺伝子型や、未知の宿主の発見に繋がるかもしれない。

SRA Taxonomy Analysis Tool

NCBIは、次世代シーケンシングデータのレポジトリであるSequence Read Archive (SRA) 上の各配列ファイル (Run) についてk-mer解析を実施している (SRA Taxonomy Analysis Tool)。各Runのk-mer解析結果を調べれば、何の生物由来のk-mer (配列) がどれだけ含まれているか分かる。
例: 熱帯魚ベタのゲノムシーケンスデータ (DRR320246)

0.85%がウイルス由来の配列であることに注意。

しかし、SRAには膨大な量のRunがある。上図のようにそれら一つ一つのページをチェックしていくのは不可能である。

Cloud-based Taxonomy Analysis Table

幸いにも、NCBIは任意の生物種に由来する配列を含むRunを選抜するツールとしてCloud-based Taxonomy Analysis Tableを提供している。これはGoogle Cloud PlatformのBigQueryを介したサービスのため、利用するにはBigQueryのユーザー登録が必要である。BigQueryは一定以上利用すると以後課金制となるため、無料ユーザーは一ヶ月に数回しか検索を走らせられない。データ検索にはSQLという形式の呪文 (クエリ) を覚える必要がある。Cloud-based Taxonomy Analysis Tableを用いたウイルス探索は、新型コロナウイルスの探索に使用実績があるようだ。

やること

本稿では、Cloud-based Taxonomy Analysis Tableを使ってウイルス配列を含んだRunを特定し、配列データをダウンロードして解析する流れを紹介する。今回は魚に感染する大型DNAウイルスのメガロサイティウイルス (Megalocytivirus) に由来する配列を探してみる。

• ウイルス配列を含むRunを検索する
  • BigQueryを開く
  • クエリの実行
• データ解析
  • ダウンロード
  • トリミング
  • アセンブリ
  • 参照ゲノム配列にマップしてバリアントコールする
• 考察と注意点
• References

ウイルス配列を含むRunを検索する

BigQueryを開く

BigQueryのウェブサイトを開く。
「コンソールへ移動」をクリック。

「＋」ボタンをクリックしてクエリを新規作成する。

クエリの実行

エディタに以下のSQLクエリをコピペし実行する。

SELECT
  m.acc,
  m.assay_type,
  m.consent,
  m.sample_name,
  m.platform,
  m.instrument,
  m.librarylayout,
  m.librarysource,
  m.organism,
  m.biosample,
  m.sra_study,
  m.releasedate,
  m.bioproject,
  m.mbytes,
  m.collection_date_sam,
  m.geo_loc_name_country_calc,
  tax_id,
  tax.name,
  tax.total_count
FROM
  `nih-sra-datastore.sra.metadata` AS m,
  `nih-sra-datastore.sra_tax_analysis_tool.tax_analysis` AS tax
WHERE
  m.acc=tax.acc
  AND tax.tax_id=249585
  AND total_count>1
ORDER BY
  releasedate DESC

FROM以下でどのデータベースを検索するか指定。
WHEREでヒットを絞り込む条件を指定。
得られたヒットについてSELECT以下に指定した情報を1行で出力する。
筆者は以下の情報を出力している。その他に利用可能な情報やそれぞれの詳細についてはこちらを参照されたい。
m.acc Runのアクセッションナンバー
m.assay_type どういう解析か (全ゲノムショットガン WGS、RNA−Seqなど)
m.consent そのデータセットは自由に使えるか。Public
m.sample_name サンプル名
m.platform シーケンサーのタイプ (Illumina, PacBio, ONTなど)
m.instrument シーケンサーの機種 (Illumina MiSeq, PacBio Sequel II, ONT MinION など),
m.librarylayout ペアエンド (PAIRED) かシングルエンド (SINGLE) か
m.librarysource ゲノムDNAかトランスクリプトーム (RNA) か
m.organism 由来生物種名
m.biosample BioSample ID
m.sra_study SRA Study accession in the form of SRP######## (ERP or DRP for INSDC partners)
m.releasedate データ公開日
m.bioproject BioProject ID
m.mbytes データ量 (メガバイト)
m.collection_date_sam サンプリングした日時
m.geo_loc_name_country_calc サンプリングした地名 (国名や州、都道府県レベルが多い)
tax_id クエリ生物種Taxonomy ID
tax.name クエリ生物種名
tax.total_count クエリ生物種由来k−merの数

tax.tax_idには対象ウイルスのNCBI Taxonomy IDを入れる。経験上、属レベルのTaxonomy IDを入れると特異性とヒット数のバランスが良い感じの結果が得られる。種レベルまで絞るとヒットがかなり少なくなってしまう。
ORDER BYで並び替え順を指定。

  total_count DESC

にするとクエリ生物種由来k−merの数が多い順に並び替えられる。

以下のような結果が出た。

BigQueryではクエリを投げるごとにクレジットを消費してしまう。5回位回すと一ヶ月分の無料クレジットを使い切ってしまうようだ。それ以上回したい場合は課金するか別アカウントを使う。

CSVファイルをローカルにダウンロードする。

ダウンロードしたCSVファイルをExcelで開く。

経験上、ゲノムサイズ1万塩基あたり10000mer 程度のカバレッジがあれば全ゲノムをアセンブルできる。メガロサイティウイルスのゲノムサイズは約11万塩基なので、ウイルス由来のk−merを110000 mers以上含んでいるサンプルを解析してみる。

ランキング上位にはメガロサイティウイルスの一種RBIV (rock bream iridovirus) に感染させたサンプル由来のトランスクリプトームデータが並ぶ。トランスクリプトーム (=RNA配列のデータ) からウイルスのDNAゲノム配列をアセンブルすることが可能な場合もあるが、今回は除外する。グラスフィッシュのゲノムシーケンスデータがかなり有望だが、ファイルサイズが大きい (11000 MB = 11GB) のでこれも除外する。

メガロサイティウイルス配列を大量に含んでおり、かつファイルサイズが手頃なサンプルとして、熱帯魚ベタ (Betta splendens) のIlluminaゲノムシーケンシングファイル (DRR320246) があった。

acc	assay_type	consent	sample_name	platform	instrument	librarylayout	librarysource	organism	biosample	sra_study	releasedate	bioproject	mbytes	collection_date_sam	geo_loc_name_country_calc	tax_id	name	total_count
DRR320246	WGS	public	SAMD00407021	ILLUMINA	NextSeq 500	PAIRED	GENOMIC	Betta splendens	SAMD00407021	DRP006411	2021-12-13 00:00:00.000000 UTC	PRJDB7253	2837	2021/1/1	Singapore	308906	Megalocytivirus	166635

今回はこのファイルを解析してみる。

データ解析

ダウンロード

NCBIからNGSデータをダウンロードするとちょっとした下処理が必要になる。そこでEuropean Nucleotide Archive (ENA) から.fastq.gzファイルを直接ダウンロードする。

wgetでダウンロードする。

nohup wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR320/DRR320246/DRR320246_1.fastq.gz &
nohup wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR320/DRR320246/DRR320246_2.fastq.gz & 

トリミング

生リードは不正確な (低品質な) 配列を含んでいる場合があるため、そうした配列を取り除く作業 (クオリティトリム) を行うことが勧められる。クオリティトリムに用のプログラムにはTrimmomatic, cutadapt、Trim Galore!、RabbitQCPlusなどいろいろある。今回はFastpでトリム。

fastp -i DRR320246_1.fastq.gz -I DRR320246_2.fastq.gz -o DRR320246_1_trim.fq.gz -O DRR320246_2_trim.fq.gz --html DRR320246_fastp.html --thread 8 --trim_poly_x

アセンブリ

生リードは、元のウイルス配列を細切れにした断片でしかない。そこで、細切れの配列を元のゲノム配列に復元する作業(de novoアセンブリ)を行う。
今回用いるIlluminaリードはショートリードと呼ばれるタイプのデータ。ショートリード用のアセンブラとしてポピュラーなものにPlatanus, MEGAHIT、ABySSなどいろいろある。今回はSPAdesでアセンブルする。

spades.py --meta --only-assembler \
-o DRR320246_spades_k127_meta \
-1 DRR320246_1_trim.fq.gz \
-2 DRR320246_1_trim.fq.gz \
-t 20 -m 50 -k 23,35,57,79,91,103,115,127 

アセンブリ結果 (assembly_graph_with_scaffolds.gfa) をBandageで可視化する。

断片的な宿主配列に交じって、明らかにカバレッジの厚い約110 kbの環状配列が得られた。

この環状配列に対応するコンティグがscaffolds.fastaにあった。

scaffolds.fasta (part)

>NODE_14_length_111283_cov_30.639057
CTCTTGGGTGGCGGCGTGGTTGTAGGCGGCGCGCTGGTGATATCCTCATACTCTGGACAG
TCGTACTCGTCCTCCTCTTCCTGCGGGTCCTCAGGCTCGGCCTCGGTGGGTAGCGGTATT
GCGACTGTGCACATTAGTCCTGTGAAGCCTTCTTTACAGGCACACTGTCCGGCTTCATTA
CACTGCAGTGACATTGATTGAGTGGGATGACACATGCATGCTTCACAGCCTGTCTGTATA
TTGAAATATCCAGGCATACATTGATCACATTGCAAGCCCATCACATTAGGCTTGCACATG
CATTGTCCAGTGTACTCGTTGCACAGTGGATGCAATGAGCCGTGTCCACTGCAAACACAT
GGTTTGCACTTGTTGTTAGGGTCAGGACTCAGGGCGTTGCCATAATAGCCCCTGGTGCAT
TTGTCACAGAAGAATCCATCGGTATTATGCAGGCAGCGCAGGCATTCGCCAGTCTCAGGA
TTACACTGTCCTACTGGGTTGTCGTCAATGTTGCCTGAACATTGGCAGCGCCTGCACGGG
CGAACAGGGCCTCGCAAGCCACTAGGGTCGCCAAATGACCCATCGAGGCACTCTTCGCAT
GCGTTGTCCTTTGTGGTATTGATAGGGCAGCCCGTACATATGATGTTCTCTCCAGTCTGT
CTACACGGCCCATTGTGAGGGCACGGACATGGTATGCATCCATCTGAAACATGGAAAAAG
CCAGGTGCACACTGTTCGCAGTATCTCCCGACGAAGCCAGTTGGGCAGTCGACACAATGT
ACACCACTATTGTCGTCGGTATAGTTCATAAAGCAGGGTCCTCCTCCTCGACATTGGCTA
CATTGCATACAGACCACTTGGCCATTGTGTTGCAGCGCACATTGTCCATCGCATGGGCAG
CGGCTACAATACTCTCCTGTCCAGAAGTAGCCCATTTGACATACATCACACTGTGGGCCA
GCATATCCAACTGGACATACACAGACTATTGTACCATTCTCCGACATACACGGGCCATTG
AACGGACATGCACATTCCCGGCACACACCATCACGGGCTTTATAATATCCATTCAGACAG
CTGTCGCATGTGTCTCCAATGTGACCCGGTGGACAGCTTGAACATACCAGGCGACCCATA
TGGACTTCACACGGACCATCTTCTGGACATGGGCATTGTTGGCATGAAATCGGGTCGCCG
TAGTAGCCGTCTGGACAGCACACGTGGGCCCCATCATAGCTTACTCGACATAGTCCATTG
TCACAATGGCATGTGACATTCGATTGGTAGCACTGGAGCATATTGTAGGGCATCTCACCG
GGGACACACTCATCCAGACACTCCAAAGAGCCATTTGTGTATCTCACCACATTTTCACAG
TCACTGCAGTTGCCGCTCAAACACTCTGGCTCATCTATGTCATCGTAGTCGTCCATTCCG
CTGCCCCCATCGTCAAGCAGTGTAGGTGGTGGAGTAACGTCATCGGTGTCTGTCGGCAGC
TCACATGAGACGCCTACACAAGGCTGACTGTCAGATGAGATGCGGCTGGCGTGGCATGTG
ACGCTCTGCACAGGGTGAGGTTTCAGCTTGATGACAGACACAATGGTACCATCATACAGC
ACCACTCCATGTTTCAGGACTTCGCTGCTGTTGCGGCCTACATGTACCACCTCGCCATGT
ACAACATGCTCTGCCAACAGGCTGTTGCTGTCGCTTGCCCAAACAATCTTCACATCAGTC
...

NODE_14_length_111283_cov_30.639057をSnapGene Viewerで可視化する。

配列上にタンパク質をコードする遺伝子 (open reading frames; ORFs) が密に配置している (オレンジと緑の矢印)。ウイルスゲノムが得られたとみて間違いない。
配列の一部をNCBI BLASTNで相同性検索すると、メガロサイティウイルスのなかでも東南アジアから中国にかけて分布する伝染性脾腎臓壊死症ウイルス (infectious spleen and kidney necrosis virus: ISKNV) に99％一致した。紛れもなくISKNVの配列だ。

ISKNVはスズキ目魚類に広く感染する。ベタもスズキ目魚類の一種であり、ISKNVに感染することは十分考えられる。このことから、今回得られたISKNVのゲノム配列は実際にベタに感染していたウイルスに由来する可能性が極めて高い。

ISKNVをはじめとするメガロサイティウイルスは、養殖魚の大量死を引き起こすマダイイリドウイルス病の病原体。ISKNVは日本では未確認のようだが、輸入観賞魚を媒介として日本に侵入する可能性は十分ありそうだ。定着可能性に関わらず、観賞魚を野外に放流してはいけない。

参照ゲノム配列にマップしてバリアントコールする

すでに参照ゲノムが利用可能なウイルスの場合、

• 参照ゲノムにリードを貼りつける (マッピング)
• 貼り付いたリードの情報をもとにして、該当サンプルと参照ゲノム間で変異が生じている箇所を特定する (バリアントコール)
• バリアントコールをもとに参照ゲノム配列を書き換える (コンセンサス配列の取得)

ことで、アセンブルすることなくサンプル由来ウイルスのゲノム配列 (コンセンサス配列) を得ることができる。
今回はリードのマッピングにMinimap2、バリアントコールとコンセンサス配列の取得にはbcftoolsを使う。

# Map reads
(minimap2 --MD -ax sr -t 8 ISKNV.fasta DRR320246_1_trim.fq.gz DRR320246_2_trim.fq.gz | samtools view -@ $NUMTH -G 12 -bS - |samtools sort -@ 8 -O BAM -o DRR320246.bam )2>DRR320246.minimap2.log
samtools index DRR320246.bam

# Call variants
bcftools mpileup  -f ISKNV.fasta DRR320246.bam | bcftools call --ploidy 1 -c -Oz -o DRR320246.vcf.gz

#indexing
bcftools index DRR320246.vcf.gz 

#export consensus sequence
cat ISKNV.fasta | bcftools consensus DRR320246.vcf.gz > DRR320246.fa

この方法では塩基置換や小さな挿入/欠失といった小規模な変異を反映したコンセンサス配列を得られる。しかし、遺伝子の挿入や転座といった大規模な変異はキャプチャできない可能性があることに注意したい。

考察と注意点

k−mer解析は核酸配列レベルで既知の配列に似たウイルスしか検出できない。そのため、既知のウイルスとアミノ酸配列レベルでしか似ていない、例えば属や科レベルで新しいウイルスは発見できない。

結果が生物学的に妥当か、自然感染なのかコンタミなのか見極めることも重要だ。例えば植物のゲノム・トランスクリプトームデータから動物ウイルスの配列が大量に検出された場合、コンタミあるいはメタデータの間違いを疑うべきだ。

NCBI SRAはINSDCという国際的なデータベースの一部であるため、日本のDDBJおよび欧州のENAにアップロードされた配列もカバーしている。しかしながら、INSDCを避けて自国ローカルのデータベースにしかデータをアップロードしない事例が近年散見される。そのようなデータについてはBigQueryで掘ることはできない。

今回紹介した方法は特定のウイルスに焦点を絞った方法であり、新奇のウイルスを含め大規模に探索するには別の手法が必要となる (References参照)。とはいえ、BigQueryを用いたウイルス探索はなにかと「小回りが利く」手法であるといえるだろう。

References