Glowを用いた大規模ゲノミクスデータセットに対する変異体正規化の効率化

[Streamlining Variant Normalization on Large Genomic Datasets](https://databricks.com/blog/2019/12/05/streamlining-variant-normalization-on-large-genomic-datasets-with-glow.html)の翻訳です。

:::note warn
本書は抄訳であり内容の正確性を保証するものではありません。正確な内容に関しては原文を参照ください。
:::

:::note info
**注意**
2019年の記事です。
:::

[Glow blog](https://glow.readthedocs.io/en/latest/blogs/variant-normalization-blog/variant-normalization-blog.html)とのクロスポストです。

ゲノミクスの世界における数多くの創薬研究プロジェクトには、大規模な遺伝子変異体のデータセットが関わっており、過去10年間で指数関数的にボリュームが増大しています。しかし、これらのデータセットに対して抽出、変換、ロード(ETL)、そして分析を行うためのツールは、この成長について行けていません。シングルノードのコマンドラインツールやスクリプトは、これらのプロジェクトのテラバイトの遺伝子データを取り扱うには非効率的すぎます。今年の10月に、[Databricks](https://databricks.com/jp)と[Regeneron Genetics Center](https://www.regeneron.com/science/genetics-center)は、この問題に取り組むために、Apache Spark™をベースとした大規模遺伝子分析のためのオープンソースツールキットである[Glow](https://projectglow.io/)プロジェクトを導入するためにパートナーシップを結びました。Glowは、いくつかの機能の最適化に加え、スケーラブルかつマネージドなクラウドサービスの上で数多くの2次解析、3次解析の機能を提供する[Databricks Unified Data Analytics Platform \(UDAP\) for Genomics](https://databricks.com/product/genomics)と連携しており、DatabricksがGlowの最適なプラットフォームとなっています。

大規模なチーム横断研究や創薬プロジェクトにおいては、コンピュータに関わる生物学者やバイオインフォマティシャンは、後段の分析を実行するために通常非常に大きな変異体コールセットをマージする必要があります。以前の記事では、[ETLにおけるGlowの力とシンプルさをご説明](https://databricks.com/blog/2019/06/26/scaling-genomic-workflows-with-spark-sql-bgen-and-vcf-readers.html)し、GlowのVCFとBGENデータソースを用いて、これまでに無いような規模の異なるソースの変異体コールセットのマージをお見せしました。異なるソースの変異体コールセットは別の大きな問題を突きつけます。これらのセットが異なる変異体コーリングツールや方法で生成されることはよくあることです。結果として、異なる変異体コールセットにおいて、同じ遺伝子変異体が(遺伝子の位置や対立遺伝子の観点で)異なる表現となることがあり得ます。これらの変異体表現の不一致は、更なるデータ分析の前に解決する必要があります。これは以下の理由から重要なこととなります。

1. 変異体のマージされたセットに対する後段の分析の結果に適切でないバイアスが生じるのを防ぐため。あるいは、正規化の欠如によって、新たな変異体と見えるものに対する分析の労力が無駄となり、実際これは冗長なものとなります(1000のゲノムプロジェクトの変異体コールとdbSNPにおける冗長性の例については[Tanらによる論文](https://academic.oup.com/bioinformatics/article/31/13/2202/196142)をご覧ください)。
1. マージされたデータセットと解析後の派生物に互換性があり、他の公開・非公開変異体データベースと比較できるようにするため。

これは、異なるデータセットにおいて同じ変異体が同じ表現となることを保証するプロセスである変異体正規化と呼ばれる処理で実現されます。人気のあるシングルノードのツールを用いて、大規模プロジェクトのテラバイトの変異体データに対して変異体正規化を実行しようとすると、これらのツールの入出力は変異体コールの格納に用いられるフラットなファイルフォーマット(VCFやBGEN)であるため、非常に困難なものとなります。この問題に対応するために、正規化された変異体のデータフレームを生成するために、Sparkデータフレーム上で直接動作するGlowに変異体正規化変換処理を組み込み、たった1行のPython、Scalaコードを追加するだけで、高速かつスケーラブルな方法で数十万のサンプルの変異体を正規化するためにSparkのパワーと組み合わせるようにしました。我々のノーマライザーに取り組む前に、変異体正規化が実際に何をするのかについて技術的に見ていきましょう。

# 変異体正規化は何をするのか？

変異体正規化は、変異体の表現が「倹約で」かつ「左に揃っている」ことを確実にします。いかなる対立遺伝子の長さをゼロに削減することなしに、可能な限り少ないヌクレオチドで表現できる時に変異体は倹約的であると言えます。図1に例を示します。

![](https://databricks.com/wp-content/uploads/2019/12/Figure-1.png)
図1: 変異体の倹約性

全ての対立遺伝子の長さを同じに保ちながらも、これ以上左に動かせない位置にある時に変異体が左に揃っていると言えます。図2に例を示します。

![](https://databricks.com/wp-content/uploads/2019/12/Figure-2.png)
図2: 左に揃った変異体

**Tanら**は、正規化は唯一の結果になることを証明しました。言い換えると、実際に異なる変異体であるときのみ、2つの変異体は異なる正規化表現を持つということになります。

# Glowにおける変異体正規化

我々はGlowにnormalize_variantsトランスフォーマーを導入しました(図3)。VCF、BGENあるいはDeltaリーダーを用いて変異体コールをSparkデータフレームに取り込んだ後で、ユーザーは全ての変異体を正規化するために1行のPythonあるいはScalaコードを呼び出すことができます。これによって、全ての変異体が正規化された状態で表現される別のデータフレームを生成します。正規化されたデータフレームは、[我々のビルトインの回帰関数](https://databricks.com/blog/2019/09/20/engineering-population-scale-genome-wide-association-studies-with-apache-spark-delta-lake-and-mlflow.html)や[効率的に並列化されたGWASツール](https://databricks.com/blog/2019/10/02/parallelizing-saige-across-hundreds-of-cores.html)を用いたGWASのような後段の解析に使用することができます。
![](https://databricks.com/wp-content/uploads/2019/12/Figure-3.png)
図3: Glowを用いたスケーラブルな変異体正規化

normalize_variantsトランスフォーマーはこの重要な上流プロセスに、これまでにないようなスケーラビリティとシンプルさをもたらし、これによって、Databricks UDAP for GenomicsとGlowが、数十万人のグループに対する遺伝変異と疾病の関連性の研究のようなバイオバンク規模の遺伝解析における理想的なプラットフォームとなっています。

# 内部の正規化アルゴリズムと精度

異なる正規化アルゴリズムを使用する変異体正規化のシングルノードツールがいくつか存在しています。広く使われている変異体正規化ツールには、[vt normalize](https://genome.sph.umich.edu/wiki/Vt)、[bcftools norm](http://www.htslib.org/doc/bcftools.html#norm)、GATKの[LeftAlignAndTrimVariants](https://gatk.broadinstitute.org/hc/en-us)があります。

我々の調査、そして[Bayatら](https://academic.oup.com/bioinformatics/article/33/7/964/2623048)、[Tanら](https://academic.oup.com/bioinformatics/article/31/13/2202/196142)によって示されているように、GATKのLeftAlignAndTrimVariantsは、いくつかの変異体を完全に左揃えにすることに頻繁に失敗します。例えば、GATK自身のテストファイルである[test\_left\_align\_hg38\.vcf](https://github.com/broadinstitute/gatk/tree/master/src/test/resources/org/broadinstitute/hellbender/tools/walkers/variantutils/LeftAlignAndTrimVariants)では、LeftAlignAndTrimVariantsの適用によってファイルの16変異体のうち、chr20:63669973、chr20:64012187、chr21:13255301の位置を含む3つで不適切な正規化が行われていることに気づきました。これらの変異体はvt normalizeやbcftools normでは適切に正規化されます。

結果として、我々のnormalize_variantsトランスフォーマーでは、根本的には類似しているbcftools normやvt normalizeアルゴリズムの改善版を使用しました。対象の変異体に対して、右側のヌクレオチドが同じである限り、変異体の全ての対立遺伝子の右側をトリミングするところからスタートします。どれかの対立遺伝子の長さがゼロになった場合、リファレンスゲノムから固定ブロックのヌクレオチドで左側に追加します(リファレンスゲノムのリファラルの数を制限するために1つづつ追加するのではなく、ヌクレオチドがブロックで追加されます)。右側のトリミングが終了すると、全ての対立遺伝子で共通する左側のヌクレオチド(場合によっては事前の左側の追加によって生成されます)を削除するために、潜在的な左側のトリミングが実行されます。このプロセスを通じて変異体の開始、終了、対立遺伝子が適切にアップデートされます。

複数のテストファイルに対してvt normalizeとbcftools normとの精度のベンチマークを行い、これらのツールの結果とマッチすることを検証しました。

# 複対立変異体の分割のオプション

normalize_variantsのトランスフォーマーは、オプションでmultiallelic変異体を二対立遺伝子に分割することができます。このトランスフォーマーに指定されるモードオプションでコントロールすることができます。以下のモードオプションを指定することができます：正規化のみを実行するnormalize(デフォルト)、正規化を実行する前にmultiallelic変異体を二対立遺伝子に分割するsplit_and_normalize、正規化を行わずにmultiallelicを分割するsplitがあります。

トランスフォーマーの分割ロジックは、GATKのsplitMultiallelicsオプションを用いたLeftAlignAndTrimVariants[ツール](https://gatk.broadinstitute.org/hc/en-us)による分割ロジックによるものと同じものです。より詳細には、VCFファイルからロードされたmultiallelic変異体を分割する際には、このトランスフォーマーは可能であれば結果の二対立遺伝子に対するGTブロックを再計算し、AC、AN、AF以外のすべてのINFOフィールドを削除します。存在する場合には新たに計算されたGTブロックに基づいてこれらの3つのフィールドが挿入され、そうでない場合にはこれらのフィールドも削除されます。

# normalize_variantトランスフォーマーの使用

ここでは、Glowを用いてどのように非常に大きな変異体コールセットを正規化、分割できるのかを簡単にデモします。最初に、[以前の記事](https://databricks.com/blog/2019/06/26/scaling-genomic-workflows-with-spark-sql-bgen-and-vcf-readers.html)でデモしたように、VCFやBGENファイルをSparkデータフレームに読み込むことができます。これは、`/databricks-datasets/genomics/call-sets`というフォルダに格納されているVCFセットに対する処理のPythonコードを示します。

```py:Python
original_variants_df = spark.read\
  .format("vcf")\
  .option("includeSampleIds", False)\
  .load("/databricks-datasets/genomics/call-sets")
```

図4にoriginal_variants_dfデータフレームのサンプルを示します。

![](https://databricks.com/wp-content/uploads/2019/12/Figure-4.png)
図4: 変異体データフレームoriginal_variants_df

この変異体は、以下のようにnormalize_variantsトランスフォーマーを用いて正規化されます。

```py:Python
import glow

ref_genome_path = '/mnt/dbnucleus/dbgenomics/grch38/data/GRCh38.fa'

normalized_variants_df = glow.transform(\
  "normalize_variants",\
  original_variants_df,\
  reference_genome_path=ref_genome_path\
)
```

正規化には、`reference_genome_path`オプションを用いて指定されるリファレンスゲノムの`.fasta`や`.fa`ファイルが必要となります。`.dict`や`.fai`ファイルは同じフォルダにリファレンスゲノムを格納する必要があります(これらのファイルフォーマットの詳細については[こちら](https://gatk.broadinstitute.org/hc/en-us/articles/360035531652?id=11013)を参照ください)。

正規化後のサンプルデータフレームは図5のようになります。
![](https://databricks.com/wp-content/uploads/2019/12/Figure-5.png)
図5: original_variants_dfに対してnormalize_variantsトランスフォーマーを適用した後に得られるnormalized_variants_dfデータフレーム

いくつかの変異体が正規化され、これらの開始、終了、対立遺伝子がそれぞれ変更されていることに注意してください。

デフォルトでは、トランスフォーマーは図5にあるように正規化の前にmultiallelic変異体を分割せずに、それぞれの変異体を正規化します。モードオプションを`split_and_normalize`に設定することで、二対立遺伝子に対しては変更はありませんが、multiallelic変異体は最初に適切な数の二対立遺伝子に分割し、結果の二対立遺伝子が正規化されます。これは以下のように実行されます。

```py:Python
split_and_normalized_variants_df = glow.transform(\
  "normalize_variants",\
  original_variants_df,\
  reference_genome_path=ref_genome_path,\
  mode=“split_and_normalize”
)
```

結果のデータフレームは図6のようになります。
![](https://databricks.com/wp-content/uploads/2019/12/Figure-6.png)
図6: original_variants_dfに`mode=“split_and_normalize”`でnormalize_variantsトランスフォーマーを適用することで得られるsplit_and_normalized_variants_dfデータフレーム

original_variants_dfのtriallelic変異体(`chr20, start=19883344, end=19883345, REF=T, ALT=[TT,C]`)が、二対立遺伝子に分割され、正規化されることで2つの二対立遺伝子(`chr20, start=19883336, end=19883337, REF=C, ALT=CT`)と(`chr20, start=19883344, end=19883345, REF=T, ALT=C`)が得られていることに注意してください。

上で述べたように、モードオプションをsplitに設定することで、トランスフォーマーは正規化を行わずにmultiallelic変異体の分割のみに使用することもできます。

# まとめ

Glowのnormalize_variantsトランスフォーマーを使うことで、コンピュータに関わる生物学者やバイオインフォマティシャンは、非常に高速かつスケーラブルな方法で数十万のサンプルの非常に大規模な変異体データセットを正規化することができます。1行のコードで、VCF、BGENリーダーを用いてソースの異なるコールセットを取り込み、マージし、このトランスフォーマーを用いて正規化を行うことができます。オプションとして、トランスフォーマーをmultiallelic変異体を二対立遺伝子へ分割する目的で使用することもできます。

# Glowを使ってみる - 変異体正規化の効率化

normalize_variantsは、非常に小規模なコード([Azure](https://glow.readthedocs.io/en/latest/etl/variant-normalization.html)|[AWS](https://glow.readthedocs.io/en/latest/etl/variant-normalization.html))で大規模な変異体データセットの正規化(そして分割)を容易にします。Glowの機能については、[こちら](https://glow.readthedocs.io/en/latest/index.html)を参照いただくか、[Databricks Unified Data Analytics for Genomics](https://databricks.com/product/genomics)をチェックするか、[すぐにでもプレビューを開始](https://pages.databricks.com/genomics-preview.html)してください。

# 次のステップ

大規模な複数のサンプルの変異体コールをDatabricksとGlowがどのようにシンプルにするのかに関してはウェビナーAccelerate and Scale Joint Genotyping in the Cloud(リンク切れ)をご覧下さい。

**リファレンス**
Arash Bayat, Bruno Gaëta, Aleksandar Ignjatovic, Sri Parameswaran, [Improved VCF normalization for accurate VCF comparison](https://academic.oup.com/bioinformatics/article/33/7/964/2623048), Bioinformatics, Volume 33, Issue 7, 2017, Pages 964–970

Adrian Tan, Gonçalo R. Abecasis, Hyun Min Kang, [Unified representation of genetic variants](https://academic.oup.com/bioinformatics/article/31/13/2202/196142), Bioinformatics, Volume 31, Issue 13, 2015, Pages 2202–2204

**追加リソース**
- [Docs > Glow](https://glow.readthedocs.io/en/latest/index.html)
- [Scaling Genomic Workflows with Spark SQL BGEN and VCF Readers](https://databricks.com/blog/2019/06/26/scaling-genomic-workflows-with-spark-sql-bgen-and-vcf-readers.html)
- [Engineering population scale Genome\-Wide Association Studies with Apache Spark™, Delta Lake, and MLflow](https://databricks.com/blog/2019/09/20/engineering-population-scale-genome-wide-association-studies-with-apache-spark-delta-lake-and-mlflow.html)
- [Parallelizing SAIGE Across Hundreds of Cores](https://databricks.com/blog/2019/10/02/parallelizing-saige-across-hundreds-of-cores.html)

### Databricks 無料トライアル

[Databricks 無料トライアル](https://databricks.com/jp/try-databricks)