Amazon Web Services ブログ

Amazon EMRにおいてApache SparkによるADAMおよびMangoを使用したゲノムデータセットの探索的データ分析

ゲノムシーケンスのコストは急速に低下し、公的に利用可能なゲノムデータ量はここ数年間で上昇しています。新しいコホートおよび研究では、100,000を超える個人で構成される膨大なデータセットが生成されています。同時に、これらのデータセットは、各コホート向けの様々なデータの膨大な量を生成する人口の中で遺伝的変異を抽出するよう処理されています。

ビッグデータのこの時代において、Apache Spark のようなツールはラージデータセットのバッチ処理向けの使いやすいプラットフォームを提供してきました。ただし、現在のバイオインフォマティクスパイプラインへの十分な代替としてこのようなツールを使用するには、ゲノムデータ処理向けのアクセス可能かつ包括的なAPIが必要になります。さらに、これらの処理されたデータセットの相互探査向けのサポートも必要です。

ADAMおよびMangoは、Apache Sparkにおけるラージゲノムデータセットを処理、フィルタリング、視覚化するための統合環境を提供します。ADAMにより、Apache Spark内のデータを収集および選択するためのSpark SQL、SQLインターフェイスを使用して、未加工のゲノムおよび各種データをプログラム的にロード、処理、選択することができます。Mangoは、Jupyterノートパソコン環境において、未加工および収集されたゲノムデータの視覚化をサポートしています。複数の解像度のラージデータセットから結果を出すことができます。

ADAMおよびMangoの組み合わせにより、統合環境におけるデータセットをロード、参照、拡張することができます。単一ノードのバイオインフォマティクスツールを使用して以前には不可能であったゲノムデータを相互に調べることができます。この投稿では、Amazon EMRにおけるADAMおよびMangoを設定および実行する方法を説明しています。1000 Genomesデータセットを調べるために、対話型のノートブック環境内でこれらのツールを使用する方法を実証します。これは、パブリックデータセットとしてAmazon S3で利用可能です。

Amazon EMRでADAMおよびMangoを構成

最初に、EMRクラスターを起動し、構成します。MangoではDockerコンテナを使用してAmazon EMRで簡単に実行します。クラスターの起動により、EMRでは以下のブートストラップアクションを使用して、Dockerおよび必要な起動スクリプトをインストールします。スクリプトは、/home/hadoop/mango-scripts で入手できます。

aws emr create-cluster 
--release-label emr-5.14.0 \   
--name 'emr-5.14.0 Mango example' \   
--applications Name=Hadoop Name=Hive Name=Spark \   
--ec2-attributes KeyName=<your-ec2-key>,InstanceProfile=EMR_EC2_DefaultRole \   
--service-role EMR_DefaultRole \     
--instance-groups \     InstanceGroupType=MASTER,InstanceCount=1,InstanceType=c5.4xlarge \     InstanceGroupType=CORE,InstanceCount=4,InstanceType=c5.4xlarge \   --region <your-aws-region> \   
--log-uri s3://<your-s3-bucket>/emr-logs/ \   
--bootstrap-actions \     
Name='Install Mango', Path="s3://aws-bigdata-blog/artifacts/mango-emr/install-bdg-mango-emr5.sh"

Mangoノートブックを起動するには、以下を実行します:

/home/hadoop/mango-scripts/run-notebook.sh

このファイルは、Amazon EMRのDockerでMangoを実行する必要がある環境変形の全てをセットアップします。ターミナルで、Mangoノートブックセッション向けのポートおよびJupyterノートブックトークンを確認します。EMRクラスター向けのマスターノードの公開DNS URL上のこのポートへ進みます。

1000 Genomes Projectからデータをロード

作業環境が整うと、ADAMおよびMangoを使用して、トリオ(母、父、子からのデータ)のゲノムシーケンスデータから子の中に興味深い変形を発見します。このデータは、1000 Genomes Project AWS 公開データセットから利用可能です。この分析内で、トリオ(NA19685, NA19661、およびNA19660)を表示し、親内に存在しないが子内に存在する変形を検索します。

特に、デノボ変形として知られている、親内ではなく子内に存在する遺伝子変数を特定したいと考えています。興味深いリージョンがあります。複数の混乱が継続する可能性があるデノボ変形のサイトを示すことができます。

Mango’s GitHub レポジトリ、または /opt/cgl-docker-lib/mango/example-files/notebooks/aws-1000genomes.ipynb in the running Docker container for Mango 内にこれらの例が含まれているJupyterノートブックがあります。

最初に、ADAM およびMangoモジュール、必要なSparkモジュールをインポートします。

# Import ADAM modules
from bdgenomics.adam.adamContext import ADAMContext
from bdgenomics.adam.rdd import AlignmentRecordRDD, CoverageRDD
from bdgenomics.adam.stringency import LENIENT, _toJava

# Import Mango modules
from bdgenomics.mango.rdd import GenomicVizRDD
from bdgenomics.mango.QC import CoverageDistribution

# Import Spark modules
from pyspark.sql import functions as sf

次に、Sparkセッションを作成します。変形でSQLクエリを実行するためにこのセッションを使用します。

# Create ADAM Context
ac = ADAMContext(spark)
genomicRDD = GenomicVizRDD(spark)

Spark SQLを使用した変形分析

染色体17から変形データのサブセット内にロード:

genotypesPath = 's3://1000genomes/phase1/analysis_results/integrated_call_sets/ALL.chr17.integrated_phase1_v3.20101123.snps_indels_svs.genotypes.vcf.gz'
genotypes = ac.loadGenotypes(genotypesPath)

# repartition genotypes to balance the load across memory
genotypes_df  = genotypes.toDF()

データフレーム内のカラムをプリントすることでスキーマで見ることができます。

# cache genotypes and show the schema
genotypes_df.columns

この下のタイプデータセットには、1000 Genomes Projectからの全てのサンプルが含まれています。そのため、NA19685 トリオ内にあるサンプルのみを考慮した下のタイプをフィルターし、メモリ内に結果をキャッシュします。

# trio IDs
IDs = ['NA19685', 'NA19661','NA19660']

# Filter by individuals in the trio
trio_df = genotypes_df.filter(genotypes_df["sampleId"].isin(IDs))

trio_df.cache()
trio_df.count()

次に、各変形のゲノムロケーションを決定するデータフレームに新しいカラムを追加します。これは、染色体 (contigName) および変形の開始や終了位置により定義されます。

# Add ReferenceRegion column and group by referenceRegion
trios_with_referenceRegion = trio_df.withColumn('ReferenceRegion',
                    sf.concat(sf.col('contigName'),sf.lit(':'), sf.col('start'), sf.lit('-'), sf.col('end')))

現在、デノボ変形を発見するためにデータセットを参照できます。しかし、最初に、Spark SQL でデータフレームを登録しなければなりません。

#  Register df with Spark SQL
trios_with_referenceRegion.createOrReplaceTempView("trios")

データフレームが登録されると、SQLクエリを実行できます。最初のクエリで、少なくとも1つの代替 (ALT) 対立遺伝子があるサンプルNA19685に所属する変形の名前を選択します。

# filter by alleles.これは、子向けの代替対立遺伝子がある変形名のリストです
alternate_variant_sites = spark.sql("SELECT variant.names[0] AS snp FROM trios \
                                    WHERE array_contains(alleles, 'ALT') AND sampleId == 'NA19685'") 

collected_sites = map(lambda x: x.snp, alternate_variant_sites.collect())

次のクエリでは、親が参照対立遺伝子があるサイトをフィルターします。次に、子から事前に生成されたセットによりこれらの変形をフィルターします。

# get parent records and filter by only REF locations for variant names that were found in the child with an ALT
filtered1 = spark.sql("SELECT * FROM trios WHERE sampleId == 'NA19661' or sampleId == 'NA19660' \
            AND !array_contains(alleles, 'ALT')")
filtered2 = filtered1.filter(filtered1["variant.names"][0].isin(collected_sites))
snp_counts = filtered2.groupBy("variant.names").count().collect()

# collect snp names as a list
snp_names = map(lambda x: x.names, snp_counts)
denovo_snps = [item for sublist in snp_names for item in sublist]
denovo_snps[:10]

興味深い変形がいくつかあります。メモリからゲノタイプが存続しないようにできます。

trio_df.unpersist()

アラインメントデータでの作業

多くの潜在的なデノボ変数サイトを見つけました。次に、未加工のアラインメントがデノボヒットに一致することを確認するために、このようなサイトで視覚的に検証できます。

最初に、NA19685トリオ向けのアラインメントデータ内にロードします:

# load in NA19685 exome from s3a

childReadsPath = 's3a://1000genomes/phase1/data/NA19685/exome_alignment/NA19685.mapped.illumina.mosaik.MXL.exome.20110411.bam'
parent1ReadsPath = 's3a://1000genomes/phase1/data/NA19685/exome_alignment/NA19660.mapped.illumina.mosaik.MXL.exome.20110411.bam'
parent2ReadsPath = 's3a://1000genomes/phase1/data/NA19685/exome_alignment/NA19661.mapped.illumina.mosaik.MXL.exome.20110411.bam'

childReads = ac.loadAlignments(childReadsPath, stringency=LENIENT)
parent1Reads = ac.loadAlignments(parent1ReadsPath, stringency=LENIENT)
parent2Reads = ac.loadAlignments(parent2ReadsPath, stringency=LENIENT)

この例では、 s3a:// instead of s3:// style URLs を使用することに注意してください。これは、ADAMフォーマットで、BAMファイルにアクセスするためにJava NIOを使用するためです。これを実施するには、これらのファイルにアクセスするためにHadoop Distributed File System 向けのJSR 203実装を使用しています。これには、 s3a:// protocol が必要です。この GitHub レポジトリ内の実装を表示できます。

トリオ内に3つの個人向けのデータアラインメントデータがあります。ただし、データはメモリへロードされていません。データへ速く後続アクセスするためのデータセットをキャッシュするには、キャッシュ() ファンクションを実行します:

# cache child RDD and count records
# takes about 2 minutes, on 4 c3.4xlarge worker nodes 
childReads.cache()

# Count reads in the child
childReads.toDF().count()
# Output should be 95634679

アラインメントデータの品質コントロール

ゲノムアラインメントデータの品質を視覚的に再確認するための1つの一般的な分析は、範囲分布を表示することです。範囲分布では、サンプルにある読み込み範囲のアイディアを提供します。

次に、染色体17で子アラインメント向けのサンプル範囲分布を作成します。

# Calculate read coverage
# Takes 2-3 minutes
childCoverage = childReads.transform(lambda x: x.filter(x.contigName == "17")).toCoverage()

childCoverage.cache()
childCoverage.toDF().count()

# Output should be 51252612

範囲データが計算、キャッシュされ、染色体17の範囲分布がコンピュータ計算され、範囲分布が計画されます:

# 範囲分布の計算

# 以下に進んで、SparkUI 内の進捗を確認できます。 
# <PUBLIC_MASTER_DNS>:8088 および、Spark アプリケーションの現在の実行をクリックします。
cd = CoverageDistribution(sc, childCoverage)
x = cd.plot(normalize=True, cumulative=False, xScaleLog=True, labels="NA19685")

 

これは、表示しているデータexomeデータであるため、非常にスタンダードに見えます。そのため、低範囲の多数サイトおよび100リード以上の少数のゲノムポジションを見ることができます。範囲で実行すると、次の分析向けのメモリ内のスペースをクリアするためにこれらのデータセットを存続しないようにできます。

childCoverage.unpersist()

発端者内のミスセンス変形でサイトを表示

アラインメントデータを検証し、変形をフィルターした後に、YBX2、ZNF286B、KSR1、GNA13などを含む、発端者内に潜在的なミスセンス変形がある4つの遺伝子があります。子および親の未加工読み取りのフィルターおよび表示によりこれらのサイトを視覚的に確認できます。

最初に、子の読み取りを表示します。GNA13 変数 (63052580-63052581) のロケーションにズームすると、ヘテロ接合 T から A コールが確認できます:

# missense variant at GNA13: 63052580-63052581 (SNP rs201316886)
# 作業者からデータを収集するには2分かかります
contig = "17"
start = 63052180
end = 63052981

genomicRDD.ViewAlignments(childReads, contig, start, end)

この位置に変形があり、代替対立遺伝子Aがあるヘテロ接合SNPの可能性があります。この変形が親内に現れないことを確認するために親データを見てみます:

# view missense variant at GNA13: 63052580-63052581 in parent 1
contig = "17"
start = 63052180
end = 63052981

genomicRDD.ViewAlignments(parent1Reads, contig, start, end)

これにより、この変形が、発端者のみでなく親でもあるというフィルターが確認されます。

まとめ

まとめとして、この投稿により、Amazon EMRのADAMおよびMangoを設定および実行する方法を実証しました。1000 Genomesデータセットを探索するために、対話型のノートブック環境内でこれらのツールを使用する方法を実証しました。これは、パブリックデータセットとして Amazon S3 で利用可能です。1000 Genomesデータ品質の検証、ゲノム内の興味深い変形の検証、未加工データの視覚化を通した結果の検証にこれらのツールを使用しました。

Mangoに関する詳細は、Mangoユーザーガイドをご覧ください。ご質問またはご提案については、以下でコメントを残してください。

 


その他の参考資料

この投稿が有用だとお考えの場合、 Amazon EMR および Amazon AthenaのHailによるゲノム分析, Amazon Athenaを使用したゲノムデータセットの対話型分析AWS Compute ブログAWSにおける高スループットゲノムバッチワークフローの構築: はじめに (パート4-1)をご覧ください。

 


著者について

Alyssa Marrow は、カリフォルニア大学、バークレイのRISELab および Yosef Lab の大学院生です。彼女の研究は、主にシステムの交差およびコンピュータを利用した生命工学です。これには、オミクスデータの処理およびコンピュータ計算するための拡張可能なシステムの構築および簡単な並行処理アルゴリズムが含まれています。