Predicting Geographic Population using Genome Variants and K-Means - The Databricks Blogの翻訳です。
本書は抄訳であり内容の正確性を保証するものではありません。正確な内容に関しては原文を参照ください。
2016年の記事です。
これはNorthwest Genome CenterとUniversity of WashingtonのDeborah SiegelとDatabricksのDenny LeeによるSparkを用いたADAMによるゲノム変異体解析のコラボレーションから生まれたゲスト記事です。
これはK-Means、ADAM、Apache Sparkを用いたゲノム変異体解析に関する3パートのシリーズのパート1です。
イントロダクション
過去数年間を通じて、我々はゲノムシーケンシングのコストと時間の急激な削減を目撃してきました。ゲノムシーケンスのバリエーションを理解することのポテンシャルは、一般的な病気の傾向がある人の特定、希少疾病の解決、医師による個人にパーソナライズした処方と投薬の実現と多岐にわたります。
この3パートのブログ記事では、ゲノムシーケンシングの入門書とそのポテンシャルをご説明します。ここでは、ゲノムシーケンスの差異を見るゲノム変異体解析にフォーカスし、Databricksコミュニティエディションを用い、Apache SparkとADAM(ゲノム処理のためのスケーラブルなAPIとCLI)を活用することでどのように高速に行えるのかを説明します。最後に、ゲノム変異データに対してk-meansアルゴリズムを実行し、これらの変異体に基づいて個人の起源となる地理グループを予測するモデルを構築します。
この記事では、ゲノム変異体とk-meansを用いた地域グループの予測にフォーカスします。ゲノムシーケンシングの振り返りとして、つまりゲノムシーケンシングとはを読んでいただくこともできますし、ゲノム変異体解析の並列化にある詳細を参照いただけます。
ゲノム変異体とK-Meansを用いた地理的グループの予測
Databricks Community EditionでApache Sparkを用いてADAMデータに対してK-Meansを実行することでゲノム変異体解析をデモンストレーションします。ノートブックでは、ビッグデータゲノミクスADAMプロジェクト(0.19.0リリース)を用いて、1000 genomes projectからのデータをどのように解析するのかを説明します。それぞれの人物がどの地域のグループなのかを予測するためにk-meansを試み、結果を可視化します。
準備
多くのデータサイエンスプロジェクトと同じように、最初に完了すべき数多くの準備タスクが存在します。この例では、サンプルノートブックで以下を説明します。
- サンプルVCFファイルをADAM parquetファイルに変換
- サンブルVCF / ADAM parquet内のデータを説明するパネルファイルのロード
- ADAMデータをRDDに読み込み、遺伝子型の並列処理をスタート
ADAM Parquetファイルの作成
VCFからADAM parquetファイルを作成するために、ADAMのSparkContext loadGenotypesメソッドを使用して、最初にVCFファイルをロードします。adamParquetSaveメソッドを用いることで、VCFをADAM parquetフォーマットで保存します。
val gts:RDD[Genotype] = sc.loadGenotypes(vcf_path)
gts.adamParquetSave(tmp_path)
パネルファイルのロード
VCFデータにはサンプルIDが含まれていますが、予測したいと考えている人口コードは含まれていません。この分析では、教示なしアルゴリズムを使用しますが、サンプルをフィルタリングし、予測誤差を推定するためには、依然として説明変数が必要となります。1000 genomes projectのintegrated_call_samples_v3.20130502.ALL.panelパネルファイルから、サンプルごとの人口コードを取得することができます。
ソース: 1000-genomes-map_11-6-12-2_750.jpg
以下のコードスニペットは、Sparkデータフレームpanelを作成するために、SparkのCSVリーダーを用いてパネルファイルをロードします。
val panel = sqlContext.read
.format("com.databricks.spark.csv")
.option("header", "true")
.option("inferSchema", "true")
.option("delimiter", "\\t")
.load(panel_path)
ここで使用するk-meansクラスタリングアルゴリズムにおいては、3クラスター向けにモデリングを行いますので、3つの人口グループのフィルターを作成します。イングランド、スコットランドからのイングリッシュ(GBR)、アメリカ南西のアフリカ人の先祖(ASW)、中国の北京の漢民族(CHB)です。これら3つの人口グループのみを含むfilterPanelデータフレームを作成することで、これを達成します。これは小さなパネルなので、すべてのエグゼキューターにこれをブロードキャストし、さらなるオペレーションを行う際のデータシャッフルを軽減することができ、より効率的になります。
// Create filtered panel of the three populations
val filterPanel = panel.select("sample", "pop").where("pop IN (‘GBR’, ‘ASW’, ‘CHB’)")
遺伝子型の並列処理
以下のコマンドを用いることで、これら3つの人口グループの遺伝子型をロードします。遺伝子型を含むParquetファイルがファイルレベルの述語プッシュダウンを実現しつつも、フィルタリングされたパネルはメモリーにロードされ、すべてのノードにブロードキャストされるので、より効率的に並列に処理することができます。すなわち、興味のあるレコードのみがファイルから読み込まれます。
// Create filtered panel of the three populations
val popFilteredGts : RDD[Genotype] = sc.loadGenotypes(tmp_path).filter(genotype => {bPanel.value.contains(genotype.getSampleId)})
ノートブックには、以下のような追加のステップが数多く含まれています。
- データの探索 - ここで使用するデータは約50万の塩基対をカバーする第6染色体の変異体の小規模なサブセットです。
- データのクレンジングとフィルタリング - 変異体が3対立の場合、あるいは欠けているデータ
- k-meansクラスタリング向けのデータの準備 - (全く同じ順番で変異体を含む)サンプルごとのMLベクトルの作成、モデルが処理を行う特徴量ベクトルの抽出
最終的には、データに含めた805の変異体の遺伝子型が、地理的グループを予測するために使用する特徴量となります。次のステップは、k-meansクラスタリングを実行するための特徴量ベクトルとデータフレームの作成です。
KMeansクラスタリングの実行
上の準備ステップを経ることで、ゲノムシーケンスに対するk-meansクラスタリングの実行はSpark Programming Guideのk-meansのサンプルと同じようなものとなります。
import org.apache.spark.mllib.clustering.{KMeans,KMeansModel}
これで、モデル、すなわちクラスターを手に入れることができました。人口グループを予測し、コンフュージョンマトリクスを計算しましょう。最初に、オリジナルの値(オリジナルのCHB、ASW、GBRの地理的位置をポイントする値)を持つpredictionRDDの作成タスクを実行し、特徴量(ゲノムの変異体)に基づく地域に対するモデル予測を出力するために、clusters.predictを利用します。次に、これをpredictDFデータフレームに変換し、クエリー(display()コマンドの実行、以降のセルでのRコマンドの実行など)しやすくします。最後に、オリジナルのラベル(実際の地域人口の位置)を取得するためにfilterPanelとjoinします。
// Create predictionRDD that utilizes clusters.predict method to output the model’s predicted geo location
val predictionRDD: RDD[(String, Int)] = dataPerSampleId.map(sd => {
(sd._1, clusters.predict(sd._2))
})
こちらが、予測値と実際の値の間の出力のグラフ表現となります。
コンフュージョンマトリクスをクイックに計算する例として、Rを使用します。
%r
resultsRDF <- sql(sqlContext, "SELECT pop, prediction FROM results_table")
confusion_matrix <- crosstab(resultsRDF, "prediction", "pop")
head(confusion_matrix)
以下のような出力となります。
prediction_pop CHB GBR ASW
1 2 89 30 3
2 1 14 58 17
3 0 0 3 41
ノートブック内では、オリジナルのサンプル、地域人口、人口コード、予測されたコードをjoinするための追加のSQLコードが含まれています。これによって、予測結果をそれぞれのサンプルにまでマッピングすることができます。
Lightning-Vizにおける力学モデルを用いたクラスターの可視化
これらのk-meansクラスターを可視化する興味深い方法として、Lightning-Vizによる力学モデルの利用があります。ノートブックには、lightningの可視化を生成するPythonコードが含まれています。以下のアニメーションgifでは、3つの人口グループを表現する3つのクラスター(左上:2、右上:1、下:0)が表示されています。予測されたクラスターのメンバーは、クラスターの頂点であり、色の違いが人口グループを表現します。人口をクリックすることで、サンプルのID、色(実際の人口)、予測結果(頂点に伸びる線)が表示されます。
議論
この記事では、ゲノムシーケンシングの入門(つまりゲノムシーケンシングとは)と変異体解析に関わる複雑性(ゲノム変異体解析の並列化)を説明しました。ADAMを導入することで、変異体を並列分散処理で処理することができ、解析の性能と精度を劇的に改善することができます。これは、ご自身でDatabricks Community Editionで試すことができる、Apache Sparkノートブックを用いて、ADAMデータに対するK-Meansを実行することによるゲノム変異体解析でデモンストレーションされています。ゲノム変異体解析によって、一般的な疾病の傾向がある個人を特定し、気象疾病を解決し、パーソナライズされた治療を提供することが可能となります。大規模な並列シーケンシング、大規模な並列バイオインフォマティック分析に要するコストと時間の急激な削減によって、急激に増加するシーケンスデータに対する再現可能な分析をハンドリングできるようになり、現時点では利用できるない解析手法にも貢献するかもしれません。最終的には、医療の進歩に貢献することでしょう。
リソース
ノートブックを作成する助けになった以下のリソースを特筆しておきたいと思います。
- ADAM: Genomics Formats and Processing Patterns for Cloud Scale Computing (Berkeley AMPLab)
- Andy PetrellaのLightning Fast Genomics with Spark and ADAMと関連のGitHub repo
- Neil FergusonのPopulation Stratification Analysis on Genomics Data Using Deep Learning
- Matthew ConlenのLightning-Viz project
- (Sparkによるゲノミクスにおける)Timothy Danford’s SlideShare presentations
- Centers for Mendelian Genomics uncovering the genomic basis of hundreds of rare conditions
- NIH genome sequencing program targets the genomic bases of common, rare disease
- The 1000 Genomes Project
また、Anthony Joseph、Xiangrui Meng、Hossein Falaki、Tim Hunterによる追加の貢献とレビューに感謝します。