グラフモデルと集団情報を使用した構造バリアントの正確なジェノタイピング

はじめに

構造バリアント（SV）は、個人間のヌクレオチド配列の差異のかなりの割合に寄与し、多くのヒト疾患に関与しています^{1, 2}。しかし、特にショートリード全ゲノムシーケンス（WGS）データでは、リファレンスゲノムからの大幅な偏差によりアライメントが困難になる可能性があるため、SVの正確な検出は依然として課題となっています³。最近のロングリードシーケンス技術の進歩により、SVの検出が容易になりましたが、そのアプリケーションはコスト、スループットによって制限され、相当量のDNAを必要とします（イルミナのPCRフリープロトコールでは約25ngに対して、ゲノムあたり約5ug）⁴。あるいは、サイエンティストは、テクノロジーの組み合わせを使用して構築できるリファレンスSVデータベースから構築されたグラフ法を使用して、既知のバリアントをジェノタイピングすることで、集団レベルでSVを研究することができます⁵。これらのリファレンスデータベースが進化するにつれて、SV変異のより完全な全体像を得るために、更新された情報で集団シーケンスデータを再ジェノタイプすることができます。

この記事では、ショートリードでシーケンスされた多数のサンプル集団において既知のSVをジェノタイピングできるグラフベースのジェノタイピング性能を評価します。パラグラフは、リファレンスアリルおよび可能性のある代替アリルを表すパスを持つ指向性非環状グラフを作成し、リードをこのシーケンスグラフに再アライメントし、再アライメントされたリードに基づいてSVの遺伝子型を決定します（図1）。ロングリード6でシーケンスされた3つのサンプルから構築されたブレークポイント精度SVトゥルースセットにパラグラフを適用しました。真実データの構築に用いたサンプル中のSVのジェノタイピングによる精度の評価に加えて、2,501人の無関係な個人の1,000のゲノム集団における真実セットSVのジェノタイピングを行い、集団遺伝学メトリクスを使用してSVコールの精度を定量化しました^7,15。Hardy-Weinberg equilibrium（HWE）などのメトリクスにより、集団におけるバリアントの全体的な性能を研究し、少数のサンプルの検査だけでは検出できない可能性のある系統的エラーを特定することができます⁸。1000ゲノムサンプルを用いて、集団レベルのジェノタイピング統計がコールのクリーンアップに役立つことを実証しました。さらに、異なる民族におけるこれらのSVの人口統計学的分布を定量化し、機能的ゲノム要素における精製選択のシグネチャーを特定しました。

結果

SVの真実をつくる

Genome in a Bottle（GIAB）プロジェクトデータに含まれる3つのサンプルから開始します。NA12878（HG001）、NA24385（HG002）、およびNA24631（HG005）（データおよび資料の提供可能性）。これらのサンプルのロングリードシーケンスデータは、Pacific Bioscience（PacBio）Sequelシステムを使用して生成されました。これらのサンプルは、約11,000塩基対（bp）のHiFiリードで平均30倍の深度にシーケンスされました。このデータから、PBSV ^{9を使用してSV（50bp+）を呼び出しました}。これらのSVをマージした後、38,709の固有の常染色体SVを特定しました。Paragraphは単一のSVを持つ地域で機能するように設計されているため、複数のSVが観察された地域（例：異なるブレークポイントを持つ重複または近くのSV）を除外しました。これにより、20,108のSV（9,238の欠失と10,870の挿入）が、検査および集団ジェノタイピングの目的で、当社のロングリードグラウンドトゥルース（LRGT）として残されました。これらのSVの詳細な説明は、こちらでご覧いただけます。⁶ .

単一サンプルの回収と精度

これらのSVのジェノタイピング性能を評価するため、グラフジェノタイサーのパラグラフ（図1）^{6を用いて、NA12878（63x）、NA24385（35x）、およびNA24631（40x）のショートリードデータに基づいてLRGT SVのジェノタイピングを行いました}。ロングリードSVコーラーにおけるホモ接合性アリルとヘテロ接合性アリルの分類精度は体系的に評価されていないため、ジェノタイピングの一致ではなくバリアントの有無に重点を置いてテストを行います。  したがって、パラグラフが非リファレンスジェノタイプコールを行い、SVがLRGTにも存在する場合、バリアントを真陽性（TP）と定義し、パラグラフがLRGTがそのサンプルにバリアントを含まない場所で非リファレンスアリルを呼び出す場合、偽陽性（FP）と定義します。この定義では、LRGTにはTPとして機能する38,239の個別の代替遺伝子型と、FPを計算するための22,085のリファレンス遺伝子型があります。3つのサンプルの独立したテストで平均すると、段落のリコールは0.84、精度は0.88です（表1）⁶。

タイプ	#真実のTP	リコール	#真実のFP	精度	Fスコア
削除	16,936	0.84	10,778	0.92	0.88
挿入	21,303	0.88	11,307	0.89	0.88

表1。LRGTでテストしたパラグラフの性能

問題のあるコールを調べたところ、FNの59%とFPの77%が、タンデムリピート（TR）内でオーバーラップするSVで発生することが判明しました。段落は、TR外のSVではTR内のSV（0.79）よりもはるかにリコール（0.90）が優れています。また、小規模（<200bp) SVs are much more likely to be within TRs (~75%) than larger (>1,000bp）のSV（約35%）も発見しました。さらに、当社の原著論文では、SVブレークポイントのエラー（場所やサイズが間違っているなど）がパラグラフのパフォーマンスに悪影響を与えることを示しています。

2,501の多様なヒトゲノムのジェノタイピング

集団内の多くの無関係なサンプルから得られる遺伝子型の結果により、グランドトゥルースデータが存在しない場合でも、HWEなどの集団遺伝学測定でバリアントコーラーの精度を評価できます。パラグラフジェノタイピング能力の実証として、また集団統計の研究のために、公開されている1000ゲノムシーケンスリソースから、無関係の2,501人（NA24385またはNA24631を除く）のLRGT SVのジェノタイピングを行いました。このデータセットは、アフリカ人660人（AFR）、アメリカ人347人（AMR）、東アジア人501人（EAS）、南アジア人489人（SAS）、ヨーロッパ人504人（EUR）の民族的に多様な集団を表しています¹⁰。すべてのサンプルは、150 bpペアエンドリードを用いて、Illumina NovaSeqプラットフォーム上で少なくとも30倍の深度までシーケンスしました。シーケンスデータは、DRAGEN v3.5.7bワークフロー^{14で処理されました}。

各集団およびデータセット全体におけるLRGT SVのアリル頻度（AF）分布を要約し、HWEを用いて集団遺伝学の予測と一致しないSVを特定しました（図2）。ジェノタイピングの頻度は、選択によってHWEの予測から逸脱する場合がありますが、HWEからの逸脱はジェノタイピングエラーによって引き起こされる可能性があります。集団サブストラクチャーもHWEからの逸脱につながる可能性があるため、この影響を最小限に抑えるため、各集団に対してHWEテスト（フィッシャーの正確確率検定¹¹）を個別に実施し、各SVについて集団ごとにHWEの不一致をカタログしました（図3a）。

5つの集団すべてで14,568人（73%）のSVがHWEテストに合格しました。2,904人（14%）のSVがどの集団でもHWEテストに合格していません。2,636人（13%）のSVが一部の集団でHWEテストに合格しましたが、5つすべてに合格したわけではありません（図3b）。5つの集団すべてでHWEテストに不合格となったSVのうち、78%がTRと重複しています。TRでは、ゲノムグラフの異なるパスは非常に似ている可能性があり、パラグラフで採用されている現在のアライメントはほぼ確率的であり、多数のヘテロ接合性遺伝子型につながります。これは、TRオーバーラップSVのジェノタイピング性能を改善するために、異なるグラフジェノタイピングモデルが必要であることを示唆しています。HWE不整合SVを除去した後、予想外のAFピークが約0.5でなくなることに気づきました（図2、実線）。一部の集団でHWEに一貫性がない14% SVでは、ほとんどが極めて有意なHWE p値を示さず、これらのSVが数サンプルまたは一部の集団でブレークポイント偏差を持つ可能性があることを示しています。

混合集団におけるSV分布

HWEテストで不合格となったSVをフィルタリングした後、AFの全体的な分布はすべての集団で類似しています。このAFの類似性にもかかわらず、HWEを通過するSVの遺伝子型用量に関する主成分解析（PCA）（ホモ接合性リファレンス遺伝子型および欠損遺伝子型については0、ヘテロ接合体については1、ホモ接合性代替遺伝子型については2）により、LRGT SVの集団特異性が明らかになります（図4a）。AFR、EAS、EUR、およびSASサンプルは、第1および第2の主成分の二次元空間に投影され、明確に分離されています。逆に、AMRサンプルは、混合集団で予想されるように、異なる祖先集団間でマッピングされます。

合計で525のSVが観察され、これは5つの集団すべてで固定されており、リファレンスゲノムの希少なアリルを表している可能性があります（図4b）。これらの他に、少なくとも1つの集団に固定された715のSVがあります。AFRは、最も分散した集団であることと一致して、他の集団よりも固定SVが少なくなっています。機能的アノテーションと組み合わせると、303の固定SVはエクソン性です。UBE2QL1における1,638 bpのエクソン挿入も、過去の2件の研究で非常に高頻度に報告されました^{12, 13}。特に、TOPMedによる最近の研究では、混合祖先からシーケンスされた53,581人のすべての個人でこの挿入が報告されています¹³。

また、ゲノム状況に基づいてSVの異なるAFも観察しました。一般的に、コード配列のエクソン、偽遺伝子、および非翻訳領域（UTR）内のSVは、イントロン領域および遺伝子間領域よりもAFが低いです。イントロンおよび遺伝子間領域のSVは、機能性要素（UTR、エクソン）のより極端なAFと比較して、より均一なAF分布を示します（図4c）。これらはすべて、機能的な結果をもたらす可能性のあるSVに対する精製の選択を示唆しています ⁷ 。一般的なSVは機能的要素でより枯渇していますが、TP73（全体でAF=0.16、腫瘍抑制遺伝子）、FAM110D（全体でAF=0.65、明確化すべき機能、おそらく細胞周期に関連する）、OVGP1（全体でAF=0.20、受精と早期胚発生に関連する）などの遺伝子のエクソン内には、いくつかの一般的なSVが見られます。具体的には、TP73の欠失の頻度はサブ集団で非常に異なります。AFRでは0.02と低く、EURとEASでは0.20より高い値です。

結論

この記事では、GIABコンソーシアムに含まれる3つのサンプル用に作成したSVトゥルースセットの初版を簡単にレビューしました。この真理値セットは、高品質のSVのデータベースのジェノタイピングにグラフ法を使用する能力を示しています。当社の継続的な取り組みにより、このデータベースが拡張され、より正確なジェノタイピングで段落が改善されます。多様な1000ゲノム集団でこれらのSVをジェノタイピングすることで、単一サンプルグラウンドトゥルースデータがない場合でも集団データを品質管理としてどのように使用できるかを示しました。さらに、グラウンドトゥルースデータで同定されたエクソンSVの多くが1000ゲノムサンプルでも検出され、選択の兆候を示すことを実証しました。WGSのコホートでパラグラフを用いてこれらのエクソンSVを直接ジェノタイピングすると、これらのSV/遺伝子の生物学をさらに理解するのに役立つ可能性があります。

データや資料の入手

パラグラフソフトウェアは、https://github.com/Illumina/paragraph。

NA12878, NA24385、およびNA24631 PacBioデータは、PRJNA540705, PRJNA529679およびPRJNA540706の下でSRAに沈着しました。イルミナのデータは、ENAのPRJEB35491に保管されています。

Dragenが再処理した1000ゲノムデータは、https://registry.opendata.aws/ilmn-dragen-1kgp/。