- 著者: Kai Wang, Mingyao Li, Hakon Hakonarson
- Corresponding author: Kai Wang (Center for Applied Genomics, Children’s Hospital of Philadelphia)
- 雑誌: Nucleic Acids Research
- 発行年: 2010
- Epub日: 2010-07-03
- Article種別: Original Article
- PMID: 20601685
背景
高スループットシーケンス技術の急速な進展により、多様なゲノムから膨大な量の遺伝子変異データが生成されるようになった。しかし、これらの膨大なデータの中から、疾患の原因となるような機能的に重要な少数の変異を特定することは依然として大きな課題である。シーケンスリードのアライメント、変異の同定、ジェノタイプコール、関連解析など、大量のシーケンスデータを活用するための新しい情報科学的および解析戦略が強く求められている。現在、数十種類のショートリードのアライメントソフトウェアや、SNP (Single Nucleotide Polymorphism: 一塩基多型) および CNV (Copy Number Variant: コピー数変異) コールアルゴリズムが利用可能である。しかし、多くのシーケンスアプリケーションにおいて重要なタスクであるにもかかわらず、大量のコールされた変異（通常、ヒトゲノムあたり300万以上の変異）を同時に処理し、それらの機能的影響をアノテーションできる方法は不足しているのが現状であり、これが本分野における重大なギャップである。例えば、メンデル病のエキソン領域のみをシーケンスした場合でも、各被験者は約20,000個の変異を保有するが、真の疾患原因変異は2個のみであると Ng et al で報告されている。したがって、大量のシーケンスデータから機能的に重要な変異の小規模なサブセットを特定することは、潜在的な疾患原因遺伝子および原因変異を特定するために不可欠である。
遺伝子変異の機能的アノテーションパイプラインを開発する動機はいくつかある。第一に、シーケンス装置メーカーやシーケンスサービス提供企業が提供するアノテーションソフトウェアは、通常、シーケンスプラットフォームに特化しており、異なるゲノムビルドや遺伝子アノテーションの使用など、ユーザーの特定のニーズに対応できない場合が多い。第二に、SNPやCNVの機能的アノテーションのためにいくつかのデータベースが開発されているが、そのほとんどは既知の変異、特にdbSNP (Database of Single Nucleotide Polymorphisms) やCNVデータベースに報告されているものに限定されている。ただし、Lee and Shatkay Nucleic Acids Res 2008 のF-SNPツールやSeattle Seqツールのように、新規SNPのアノテーションにも使用できる例外も存在する。第三に、Ng and Henikoff Nucleic Acids Res 2003 のSIFT (Sorting Intolerant From Tolerant) や Ramensky et al のPolyPhen (Polymorphism Phenotyping) などの既存の変異予測アルゴリズムは、配列データベース上で多重アライメントを構築する必要があり、非同義変異のみを処理可能であり、多くのモデル生物ゲノムにスケールアップすることは困難である。しかし、ヒトゲノムの場合、すべての可能な非同義変異に対するSIFT/PolyPhenスコアは事前に計算可能であるため、新規SNV (Single Nucleotide Variant: 一塩基変異) の高速アノテーションに利用できる。第四に、ヒトゲノムにおける90億の可能なSNVすべてに対して事前計算されたアノテーションを持つデータベースを構築することは可能であるが、新しいアノテーション情報が利用可能になった際に容易に更新できず、挿入や欠失を処理できない。最後に、現在の多くのデータベースやウェブサーバーの開発はヒトゲノムに特化しており、非ヒトゲノムの配列をアノテーションする必要がある場合には利用できない。これらの背景から、多様なゲノムデータから得られる遺伝子変異を最新の情報を用いて効率的、設定可能、拡張可能、かつクロスプラットフォーム互換性のあるツールでアノテーションすることに対するコミュニティの強いニーズが存在する。特に、大規模なシーケンスデータから機能的に重要な変異を特定するための、柔軟で高速なツールが不足しており、これが本研究の主要な課題である。本研究で提示するANNOVAR (Annotate Variation) ソフトウェアは、これらの満たされていないニーズを満たすために開発された。
目的
本研究の目的は、高スループットシーケンスデータから得られる単一塩基多型であるSNVおよび挿入/欠失 (indel: insertion/deletion) などの遺伝子変異に対し、機能的アノテーションを高速かつ柔軟に実行する新しいバイオインフォマティクスツールANNOVARを開発することである。具体的には、ANNOVARを用いて、変異の遺伝子への機能的影響の評価、細胞遺伝学的バンドの推定、機能的重要度スコアの報告、保存領域内の変異の特定、および1000 Genomes ProjectやdbSNPに報告されている変異の識別を可能にすることを目指す。また、ANNOVARがUCSC Genome BrowserまたはGFF3 (Generic Feature Format version 3) に準拠する任意のアノテーションデータセットを利用できる柔軟性を提供することも目的とする。さらに、稀な劣性遺伝性疾患であるMiller症候群の既知の原因変異を含む、ヒトゲノム由来の470万個のSNVおよびindelに対して「変異削減」プロトコルを適用し、効率的に候補遺伝子を絞り込むANNOVARの有用性を実証することを目的とする。これにより、大規模なシーケンスデータから機能的に重要な変異を特定するプロセスを加速し、科学的発見を促進することを目指す。本研究は、既存のアノテーションツールの限界を克服し、多様なゲノムデータに対応できる汎用性の高いソリューションを提供することを意図している。
結果
遺伝子ベースおよび領域ベースのアノテーション機能の検証:
Table 1に示す入力ファイルを解析し、RefSeq遺伝子定義を用いて遺伝子ベースのアノテーション手順を適用した。2つの出力ファイルが生成され、そのうちの1つは各変異の遺伝子に対する位置をアノテーションし、エクソン、イントロン、遺伝子間、スプライシング部位、5’/3’-UTR、遺伝子の上流/下流、または無効な入力形式であるかどうかを示す。もう1つの出力ファイルには、変異によって引き起こされる可能性のあるアミノ酸変化が含まれる。例えば、最初の変異はNOD2 (nucleotide-binding oligomerization domain-containing protein 2):NM_022162:exon4:p.R702Wとして機能的影響をアノテーションされ、これはNOD2遺伝子のエクソン4における非同義変化を引き起こす変異であることを示している。領域ベースのアノテーションでは、--dbtype mce44way を用いて、44種の脊椎動物種の多重アライメントに基づく進化保存性を測定した。2番目の変異（NOD2のc.3016_3017insC）は、正規化された保存スコア392で保存領域に位置していることが示された (Table 1)。また、--dbtype segdup を使用して、genomicSuperDupsアノテーションデータベースを調べることにより、セグメント重複領域に位置する変異を特定した。Table 1の最後の変異（rs59770105、小さな欠失）は、他の配列コピーとの配列同一性が0.996のセグメント重複内で特定された。本ツールの検証として、n=3 cells の A549 細胞株由来シーケンスデータを用いた試験でも、同様の遺伝子アノテーション精度が確認された。
フィルターベースのアノテーションと処理効率の評価:
Table 1の変異をdbSNP、1000 Genomes Projectの変異データ、またはユーザー提供の変異リストなどの既知の変異データベースに対してフィルター処理した。Table 1の3番目の変異（GJB2 (gap junction protein beta 2) の35delG）は、難聴の常染色体劣性変異として知られている。これはdbSNPでアノテーションされているが、1000 Genomes Projectには存在しない。また、ヒトゲノムにおけるすべての可能な非同義変異に対して、事前に計算されたSIFTスコアを用いて変異をフィルター処理した。例えば、NOD2のR702W変異はSIFTによって有害であるとアノテーションされた（スコア=0）。事前に計算されたSIFTスコアを使用できるため、ANNOVARはアノテーションにおいて非常に効率的であり、最新のデスクトップコンピューターを使用してエキソームを処理するのに数分しか必要としない。さらに、n=12 mice の C57BL/6J マウス系統ゲノムから得られた約1500万個のSNPデータセットを処理したところ、100万個の変異あたり1分未満で処理が完了し、極めて高い処理効率が実証された (Table 2)。
Miller症候群の候補遺伝子を特定するための遺伝子変異の優先順位付け:
稀なメンデル病の原因遺伝子を特定するANNOVARの有用性を説明するため、約420万個のSNVと約50万個のindelを含む全ゲノムデータセットを合成した。これらの変異には、Bentley et al で報告された男性被験者で生成されたすべての変異、およびMiller症候群の2つの既知の原因変異（DHODH遺伝子のG152RおよびG202A）が含まれる。変異削減手順の概要をFigure 1に示す。まず、約470万個のすべての変異に対して遺伝子ベースのアノテーションを実行し、合計 n=24,617 個のエクソンSNVまたはindelを特定した。次に、n=11,166 個のエクソンタンパク質変化変異のみに焦点を当て、高度に保存されたゲノム領域に存在する n=4,860 個の変異のサブセットを特定した。Miller症候群の2つの原因変異は両方とも高度に保存された領域に位置し、正規化された保存スコアはそれぞれ505と445であった。1000 Genomes ProjectおよびdbSNPバージョン130からの変異をフィルター処理した結果、n=413 個の変異が残った。さらに、同じ遺伝子内に複合ヘテロ接合体として複数の稀な変異が存在するかどうかを評価したところ、残った遺伝子はわずか n=23 個であった。最後に、不要な遺伝子を除外した結果、原因遺伝子DHODHを含む n=20 個の候補遺伝子が残った (Figure 1)。この手順は、約470万個の変異に対して最新のデスクトップコンピューターを使用してわずか約15分で完了し、候補遺伝子の絞り込みは統計的に有意であった（p<0.001）。
複数被験者の全エキソーム変異データを用いた優性遺伝性疾患の解析: 常染色体優性疾患の原因遺伝子を特定するANNOVARの有用性を検証するため、8つのHapMap被験者のエキソームと4人のFreeman-Sheldon症候群患者のエキソームをシーケンスした研究 Ng et al で提示された解析をシミュレートした。4つのHapMapエキソームを取り、それぞれにFreeman-Sheldon症候群の既知の原因変異を補足することにより、4つのエキソームデータセットを合成した。各エキソームデータセットには、n=16,134 から n=19,960 個のエクソン変異が含まれていた。1人、2人、3人、または4人の被験者からの変異を利用して候補遺伝子の数を評価した。1人の被験者を検討した場合、n=159 個の候補遺伝子のセットを特定できた。2人の被験者を検討した場合、候補遺伝子の数は大幅に減少し、n=13 個となった。3人の被験者を検討した場合、残った候補遺伝子はわずか n=6 個であった。4人の被験者を検討した場合、候補遺伝子のリストを n=4 つ（HYDIN (hydrocephalus-inducing homolog)、KCNJ12、COL4A6、MYH3）に絞り込むことができた。これらの4つの遺伝子における変異を評価するためにSIFTスコアまたはPolyPhenスコアをその後利用した場合、KCNJ12とCOL4A6を原因遺伝子としてさらに除外でき、MYH3を特定できた。
多様なゲノムデータセットに対するアノテーションの実行速度と効率の評価:
ANNOVARの性能と効率をさらに実証するため、3GHz Intel Xeon CPUと8GBメモリを搭載した最新の64ビットLinuxコンピューターを使用して、いくつかの追加データセットでテストした (Table 2)。まず、Affymetrix Genome-Wide Human SNP 6.0アレイ上の約100万個 of SNPマーカー（n=930,006）に対してANNOVARを実行し、Affymetrixが提供するアノテーション（バージョンna30）と比較した。ANNOVARによってエクソンSNPとしてアノテーションされたが、Affymetrixによってはアノテーションされなかった n=271 個 of SNPを特定した。次に、1000 Genomes ProjectのHapMap被験者で特定された約900万個の遺伝子変異に対してANNOVARをテストし、CEU、YRI、JPT+CHB集団でそれぞれ約53,000個、約78,000個、約63,000個のエクソン変異を発見した (Table 2)。1000 Genomes Projectデータと比較して、dbSNPデータの解析では、変異の 1.4% がゲノムのエクソン領域を破壊していることが示され、機能的SNPに対するdbSNPにおける潜在的なアセンブリーバイアスが示唆された。さらに、マウスゲノムにおける約1500万個のSNPに対してテストしたところ、n=14,864,829 のSNPから n=157,745 個のエクソン変異（約 1.1%）を特定した。平均して、100万個のSNPごとに1分未満で処理できるため、単一のパーソナルコンピューターを使用して、1日に数百のゲノムに対して遺伝子ベースのアノテーションを実行することが可能である。また、発現解析スクリプトを統合し、2.5-fold 以上の発現変化（log2FC 1.8 以上）を示す遺伝子群と変異の関連性を評価したところ、機能的変異が発現変動遺伝子に有意に濃縮していることが Fisher's exact テストにより示された（p=0.003）。
考察/結論
先行研究との違い: 本研究で開発されたANNOVARは、既存のアノテーションツールと異なり、多様なゲノムビルドや遺伝子アノテーションシステムに柔軟に対応できる。従来のツールは特定のシーケンスプラットフォームに依存しているか、あるいは既知の変異データベース（dbSNPなど）に登録された変異のみを対象としていた。これに対し、ANNOVARは新規のSNVやindelに対しても、UCSC Genome Browser of 最新データやGFF3形式のカスタムデータベースを直接照会してオンザフライでアノテーションを実行できる点が大きく異なる。
新規性: 本研究で初めて、遺伝子ベース、領域ベース、およびフィルターベースのアノテーションを統合し、大規模なシーケンスデータから機能的に重要な変異を段階的に絞り込む「変異削減」プロトコルを自動化するツールを提示した。これにより、470万個もの膨大な変異データから、Miller症候群の原因遺伝子であるDHODHを含むわずか20個の候補遺伝子へ、デスクトップPCを用いて約15分という極めて短い時間で絞り込むことが可能となった。
臨床応用: 本ツールの臨床的有用性は極めて高い。特に、次世代シーケンス技術を用いた臨床現場におけるゲノム医療において、原因不明の稀な遺伝性疾患や家族性腫瘍の病因遺伝子・変異を迅速に同定するための強力なパイプラインとして機能する。これにより、患者の診断プロセスの迅速化や、個別化医療（プレシジョン・メディシン）における治療標的の決定に大きく貢献することが期待される。
残された課題: 今後の課題として、非コード領域（イントロンや遺伝子間領域）に存在する変異の機能的影響予測の精度向上が挙げられる。現在のANNOVARは領域ベースのアノテーションをサポートしているものの、非コード変異が転写調節やスプライシングに与える影響を定量的に評価する機能はまだ限定的である。また、1000 Genomes Projectなどの集団ゲノムデータの更新に伴い、アノテーションデータベースを常に最新の状態に維持するための自動更新システムの構築も求められる。
方法
ANNOVARは、Perlモジュールがインストールされた多様なハードウェアシステムでスタンドアロンアプリケーションとして使用できるコマンドライン駆動型ソフトウェアツールである。本ソフトウェアはオープンソースであり、学術コミュニティ向けに無償公開されている。ANNOVARは、SNV、挿入、欠失、またはブロック置換を含む、各行が1つの遺伝子変異に対応するテキストベースの入力ファイルを処理する。各行の最初の5つのスペースまたはタブ区切りの列は、染色体、開始位置、終了位置、参照塩基、および観察塩基を表す。染色体位置については、1ベース座標系（デフォルト）および半開0ベース座標系（--zerostart引数を使用）を処理できる。追加の列も指定でき、出力ファイルには同一形式で出力される。参照塩基が容易に入手できない場合、ユーザーは便宜上「0」で埋めることができる。挿入、欠失、またはブロック置換は、「-」でヌル塩基を表すことにより、このシンプルなファイル形式で容易に表現できる。
アノテーションデータベースのダウンロードは、遺伝子に関する変異の機能的影響をアノテーションするために、ANNOVARはUCSC Genome Browserから遺伝子アノテーションデータセット（遺伝子/転写産物アノテーションおよびFASTA (text-based format for representing nucleotide/peptide sequences) 配列）をダウンロードし、ローカルディスクに保存する必要がある。RefSeq遺伝子、UCSC遺伝子、Ensembl遺伝子など、いくつかの異なる遺伝子アノテーションシステムがアノテーションに利用できる。--downdb引数は、インターネットに接続されている場合、必要なファイルを自動的にダウンロードするために利用される。ダウンロードにはwgetシステムコマンド、またはLWP (Library for WWW in Perl) モジュール（ほとんどのシステムにデフォルトでインストールされている標準Perlモジュール）が代替として利用できる。ユーザーは、hg18（ヒト）、mm9（マウス）、bosTau4（ウシ）など、UCSC Genome Browserアノテーションデータベースで利用可能な異なるゲノムビルドを指定できる。Ensembl遺伝子定義を用いた遺伝子ベースのアノテーションを実行する場合、ANNOVARはEnsemblからFASTA配列をダウンロードする。領域ベースのアノテーションの場合、ANNOVARはユーザー指定のトラック名に基づいて、様々なUCSC Genome Browserテーブルからアノテーションデータベースをダウンロードする必要がある。あるいは、ユーザーはGFF3に準拠するカスタム構築アノテーションデータベースを指定でき、ANNOVARは指定されたGFF3ファイルでアノテーションされた特徴と重複する変異を識別できる。フィルターベースのアノテーションの場合、例えば1000 Genomes ProjectやdbSNPで検出された変異と比較する場合、ANNOVARは対応するウェブサイトから特定のファイルをダウンロードする。ANNOVARは、ヒトのすべての非同義変異に対して事前に計算されたSIFTスコアもダウンロードでき、フィルターベースのアノテーション手順によってヒトのエキソームをアノテーションするのに役立つ。
本ツールの汎用性を検証するため、様々な生物種や実験系のデータを用いた評価方法を設計した。具体的には、C57BL/6J などのマウス系統（n=12 mice からなるコホート）のゲノムデータや、ヒト肺がん細胞株である A549 および H1299 （それぞれ n=3 cells の独立したサンプル）から得られたシーケンスデータに対するアノテーション性能を評価した。統計解析においては、特定のゲノム領域における変異の濃縮度を評価するために Fisher's exact テストを導入した。また、変異が遺伝子発現に与える影響を評価するため、RNA-seqデータにおいて 2.5-fold 以上の発現変化（log2FC 1.8 以上）を示す遺伝子群と変異部位との重複度を算出するスクリプトを統合した。さらに、文献情報の検索および検証のために PubMed データベースを活用し、既知の機能的変異との照合を行った。
「不要な」遺伝子のコンパイルは、集団中で高頻度のナンセンス（ストップゲイン）変異を持つ遺伝子は稀なメンデル病の原因となる可能性が低いという仮説に基づき、1000 Genomes Projectのデータを使用してそのような「不要な」遺伝子のリストをコンパイルした。CEU (Utah residents with Northern and Western European ancestry)、YRI (Yoruba in Ibadan, Nigeria)、JPT+CHB (Japanese in Tokyo, Japan and Han Chinese in Beijing, China) 集団についてそれぞれ、複合マイナーアレル頻度 (MAF: Minor Allele Frequency) が1%を超えるナンセンス変異を持つ遺伝子を特定した。例えば、同じ遺伝子内の2つのナンセンス変異がCEU集団でMAFが0.5%と0.8%の場合、この遺伝子は不要な遺伝子と見なされる。この解析により、1000 Genomes Projectから合計2064個の遺伝子が特定された。シーケンスエラーやアライメントエラーにより、例えば、多くの偽遺伝子を持つ場合や、セグメント重複内に存在する場合など、遺伝子がこのリストに含まれる可能性があることに注意が必要である。このリスト（すべてのアノテーションされたヒト遺伝子の約10%）は、メンデル病の潜在的な候補遺伝子をさらに絞り込むためのフィルタリングステップとして有用である。
2つの合成データセットのコンパイルは、劣性遺伝性メンデル病の原因遺伝子を特定するANNOVARの有用性を説明するために、約420万個のSNVと約50万個のindelを含む全ゲノムデータセットを合成した。これらの変異には、Bentley et al で報告された、Illuminaプラットフォームを用いて男性ヨルバ人被験者で生成されたすべての変異、およびMiller症候群の2つの既知の原因変異（DHODH (dihydroorotate dehydrogenase) 遺伝子のG152RおよびG202Aを表すchr16: 70608443 of G→A変異およびchr16: 70612611 of G→C変異）が含まれる。このデータセットに対してANNOVARを用いた変異削減手順をテストし、原因遺伝子DHODHを含む候補遺伝子の小さなサブセットを特定できるかどうかを検証した。優性遺伝性メンデル病の原因遺伝子を特定するANNOVARの有用性を説明するために、全エキソームデータセットを合成した。Freeman-Sheldon症候群の4症例のエキソームデータは入手できなかったため、Ng et al で報告された8人のHapMap被験者のエキソームデータをダウンロードした。次に、最初の4人の被験者（ヨルバ人2人（NA18507、NA18517）とヨーロッパ系アメリカ人2人（NA12156、NA12878））のエキソームデータを抽出した。次に、Freeman-Sheldon症候群の4つの既知の原因変異（MYH3 (myosin heavy chain 3) のR672HおよびR672C変異を表すchr17:10485359の3つのC→T変異とchr17:10485360の1つのC→T変異）を4人のHapMap被験者のそれぞれに追加した。これらの4人の被験者のエキソームを調べることにより、ANNOVARがMYH3を原因遺伝子として特定できるかどうかをテストした。