# Genomic Analysis Using ADAM, Spark and Deep Learning

Can we use deep learning to predict which population group you belong to, based solely on your genome?

Yes, we can - and in this post, we will show you exactly how to do this in a scalable way, using Apache Spark. We

will explain how to apply [deep learning](https://en.wikipedia.org/wiki/Deep_learning) using [artifical neural networks]

(https://en.wikipedia.org/wiki/Artificial_neural_network) to predict which population group an individual belongs to -

based entirely on his or her genomic data.

This is a follow-up to an earlier post:

[Scalable Genomes Clustering With ADAM and Spark](http://bdgenomics.org/blog/2015/02/02/scalable-genomes-clustering-with-adam-and-spark/)

and attempts to replicate the results of that post. However, we will use a different machine learning technique.

Where the original post used [k-means clustering](https://en.wikipedia.org/wiki/K-means_clustering), we will use

deep learning.

We will use [ADAM](https://github.com/bigdatagenomics/adam) and [Apache Spark](https://spark.apache.org/) in

combination with [H2O](http://0xdata.com/product/), an open source predictive analytics platform, and

[Sparking Water](http://0xdata.com/product/sparkling-water/), which integrates H2O with Spark.

## Code

In this section, we'll dive straight into the code. If you'd rather get something working before looking at the code

you can skip to the "Building and Running" section.

The complete Scala code for this example can be found in

[the PopStrat.scala class on GitHub](https://github.com/nfergu/popstrat/blob/master/src/main/scala/com/neilferguson/PopStrat.scala)

and we'll refer to sections of the code here. Basic familiarity with Scala and

[Apache Spark](https://spark.apache.org/) is assumed.

### Setting-up

The first thing we need to do is to read the names of the Genotype and Panel files that are passed into our program.

The Genotype file contains data about a set of individuals (referred to here as "samples") and their genetic

variation. The Panel file lists the population group (or "region") for each sample in the Genotype file; this is

what we will try to predict.

```scala

val genotypeFile = args(0)

val panelFile = args(1)

```

Next, we set-up our Spark Context. Our program permits the Spark master to be specified as one of its arguments.

This is useful when running from an IDE, but is omitted when running from the ```spark-submit``` script (see below).

```scala

val master = if (args.length > 2) Some(args(2)) else None

val conf = new SparkConf().setAppName("PopStrat")

master.foreach(conf.setMaster)

val sc = new SparkContext(conf)

```

Next, we declare a set called ```populations``` which contains all of the population groups that we're interested

in predicting. We then read the Panel file into a Map, filtering it based on the population groups in the

```populations``` set. The format of the panel file is described [here](http://www.1000genomes.org/faq/what-panel-file).

Luckily it's very simple, containing the sample ID in the first column and the population group in the second.

```scala

val populations = Set("GBR", "ASW", "CHB")

def extract(file: String, filter: (String, String) => Boolean): Map[String,String] = {

  Source.fromFile(file).getLines().map(line => {

    val tokens = line.split("\t").toList

    tokens(0) -> tokens(1)

  }).toMap.filter(tuple => filter(tuple._1, tuple._2))

}

val panel: Map[String,String] = extract(panelFile, (sampleID: String, pop: String) => populations.contains(pop))

```

### Preparing the Genomics Data

Next, we use [ADAM](https://github.com/bigdatagenomics/adam) to read our genotype data into a Spark RDD. Since we've

imported ```ADAMContext._``` at the top of our class, this is simply a matter of calling `loadGenotypes` on the

Spark Context. Then, we filter the genotype data to contain only samples that are in the population groups which we're

interested in.

```scala

val allGenotypes: RDD[Genotype] = sc.loadGenotypes(genotypeFile)

val genotypes: RDD[Genotype] = allGenotypes.filter(genotype => {panel.contains(genotype.getSampleId)})

```

Next, we convert the ADAM ```Genotype``` objects into our own ```SampleVariant``` objects. These objects contain just the data

we need for further processing: the sample ID (which uniquely identifies a particular sample), a variant ID (which

uniquely identifies a particular genetic variant) and a count of alternate

[alleles](http://www.snpedia.com/index.php/Allele), where the sample differs from the

reference genome. These variations will help us to classify individuals according to their population group.

```scala

case class SampleVariant(sampleId: String, variantId: Int, alternateCount: Int)

def variantId(genotype: Genotype): String = {

  val name = genotype.getVariant.getContig.getContigName

  val start = genotype.getVariant.getStart

  val end = genotype.getVariant.getEnd

  s"$name:$start:$end"

}

def alternateCount(genotype: Genotype): Int = {

  genotype.getAlleles.asScala.count(_ != GenotypeAllele.Ref)

}

def toVariant(genotype: Genotype): SampleVariant = {

  // Intern sample IDs as they will be repeated a lot

  new SampleVariant(genotype.getSampleId.intern(), variantId(genotype).hashCode(), alternateCount(genotype))

}

val variantsRDD: RDD[SampleVariant] = genotypes.map(toVariant)

```

Next, we count the total number of samples (individuals) in the data. We then group the data by variant ID and filter

out those variants which do not appear in all of the samples. The aim of this is to simplify the processing of the data and, since

we have a very large number of variants in the data (up to 30 million, depending on the exact data set), filtering out

a small number will not make a significant difference to the results. In fact, in the next step we'll reduce the

number of variants even further.

```scala

val variantsBySampleId: RDD[(String, Iterable[SampleVariant])] = variantsRDD.groupBy(_.sampleId)

val sampleCount: Long = variantsBySampleId.count()

println("Found " + sampleCount + " samples")

val variantsByVariantId: RDD[(Int, Iterable[SampleVariant])] = variantsRDD.groupBy(_.variantId).filter {

  case (_, sampleVariants) => sampleVariants.size == sampleCount

}

```

When we train our machine learning model, each variant will be treated as a

"[feature](https://en.wikipedia.org/wiki/Feature_(machine_learning))" that is used to train the model.

Since it can be difficult to train machine learning models with very large numbers of features in the data

(particularly if the number of samples is relatively small), we first need to try and reduce the number of variants

in the data.

To do this, we first compute the frequency with which alternate alleles have occurred for each variant. We then

filter the variants down to just those that appear within a certain frequency range. In this case, we've chosen a

fairly arbitrary frequency of 11. This was chosen through experimentation as a value that leaves around 3,000 variants

in the data set we are using.

There are more structured approaches to

[dimensionality reduction](https://en.wikipedia.org/wiki/Dimensionality_reduction), which we perhaps could have

employed, but this technique seems to work well enough for this example.

```scala

val variantFrequencies: collection.Map[Int, Int] = variantsByVariantId.map {

  case (variantId, sampleVariants) => (variantId, sampleVariants.count(_.alternateCount > 0))

}.collectAsMap()

val permittedRange = inclusive(11, 11)

val filteredVariantsBySampleId: RDD[(String, Iterable[SampleVariant])] = variantsBySampleId.map {

  case (sampleId, sampleVariants) =>

    val filteredSampleVariants = sampleVariants.filter(variant => permittedRange.contains(

      variantFrequencies.getOrElse(variant.variantId, -1)))

    (sampleId, filteredSampleVariants)

}

```

### Creating the Training Data

To train our model, we need our data to be in tabular form where each row represents a single sample, and each

column represents a specific variant. The table also contains a column for the population group or "Region", which is

what we are trying to predict.

Ultimately, in order for our data to be consumed by H2O we need it to end up in an H2O `DataFrame` object. Currently,

the best way to do this in Spark seems to be to convert our data to an RDD of Spark SQL

[Row](http://spark.apache.org/docs/1.4.0/api/scala/index.html#org.apache.spark.sql.Row) objects, and then this can

automatically be converted to an H2O DataFrame.

To achieve this, we first need to group the data by sample ID, and then sort the variants for each sample in a

consistent manner (by variant ID). We can then create a header row for our table, containing the Region column,

the sample ID and all of the variants. We then create an RDD of type `Row` for each sample.

```scala

val sortedVariantsBySampleId: RDD[(String, Array[SampleVariant])] = filteredVariantsBySampleId.map {

  case (sampleId, variants) =>

    (sampleId, variants.toArray.sortBy(_.variantId))

}

val header = StructType(Array(StructField("Region", StringType)) ++

  sortedVariantsBySampleId.first()._2.map(variant => {StructField(variant.variantId.toString, IntegerType)}))

val rowRDD: RDD[Row] = sortedVariantsBySampleId.map {

  case (sampleId, sortedVariants) =>

    val region: Array[String] = Array(panel.getOrElse(sampleId, "Unknown"))

    val alternateCounts: Array[Int] = sortedVariants.map(_.alternateCount)

    Row.fromSeq(region ++ alternateCounts)

}

```

As mentioned above, once we have our RDD of `Row` objects we can then convert these automatically to an H2O

DataFrame using Sparking Water (H2O's Spark integration).

```scala

val sqlContext = new org.apache.spark.sql.SQLContext(sc)

val schemaRDD = sqlContext.applySchema(rowRDD, header)

val h2oContext = new H2OContext(sc).start()

import h2oContext._

val dataFrame = h2oContext.toDataFrame(schemaRDD)

```

Now that we have a DataFrame, we want to split it into the training data (which we'll use to train our model), and a

[test set](https://en.wikipedia.org/wiki/Test_set) (which we'll use to ensure that

[overfitting](https://en.wikipedia.org/wiki/Overfitting) has not occurred).

We will also create a "validation" set, which performs a similar purpose to the test set - in that it will be used to

validate the strength of our model as it is being built, while avoiding overfitting. However, when training a neural

network, we typically keep the validation set distinct from the test set, to enable us to learn

[hyper-parameters](http://colinraffel.com/wiki/neural_network_hyperparameters) for the model.

See [chapter 3 of Michael Nielsen's "Neural Networks and Deep Learning"](http://neuralnetworksanddeeplearning.com/chap3.html)

for more details on this.

H2O comes with a class called `FrameSplitter`, so splitting the data is simply a matter of calling creating one

of those and letting it split the data set.

```scala

val frameSplitter = new FrameSplitter(dataFrame, Array(.5, .3), Array("training", "test", "validation").map(Key.make), null)

water.H2O.submitTask(frameSplitter)

val splits = frameSplitter.getResult

val training = splits(0)

val validation = splits(2)

```

### Training the Model

Next, we need to set the parameters for our deep learning model. We specify the training and validation data sets,

as well as the column in the data which contains the item we are trying to predict (in this case, the Region).

We also set some [hyper-parameters](http://colinraffel.com/wiki/neural_network_hyperparameters) which affect the way

the model learns. We won't go into detail about these here, but you can read more in the

[H2O documentation](http://docs.h2o.ai/h2oclassic/datascience/deeplearning.html). These parameters have been

chosen through experimentation - however, H2O provides methods for

[automatically tuning hyper-parameters](http://learn.h2o.ai/content/hands-on_training/deep_learning.html) so

it may be possible to achieve better results by employing one of these methods.

```scala

val deepLearningParameters = new DeepLearningParameters()

deepLearningParameters._train = training

deepLearningParameters._valid = validation

deepLearningParameters._response_column = "Region"

deepLearningParameters._epochs = 10

deepLearningParameters._activation = Activation.RectifierWithDropout

deepLearningParameters._hidden = Array[Int](100,100)

```

Finally, we're ready to train our deep learning model! Now that we've set everything up this is easy:

we simply create a H2O `DeepLearning` object and call `trainModel` on it.

```scala

val deepLearning = new DeepLearning(deepLearningParameters)

val deepLearningModel = deepLearning.trainModel.get

```

Having trained our model in the previous step, we now need to check how well it predicts the population

groups in our data set. To do this we "score" our entire data set (including training, test, and validation data)

against our model:

```scala

deepLearningModel.score(dataFrame)('predict)

```

This final step will print a [confusion matrix](https://en.wikipedia.org/wiki/Confusion_matrix) which shows how

well our model predicts our population groups. All being well, the confusion matrix should look something like this:

```

Confusion Matrix (vertical: actual; across: predicted):

       ASW CHB GBR  Error      Rate

   ASW  60   1   0 0.0164 =  1 / 61

   CHB   0 103   0 0.0000 = 0 / 103

   GBR   0   1  90 0.0110 =  1 / 91

Totals  60 105  90 0.0078 = 2 / 255

```

This tells us that the model has correctly predicted 253 out of 255 population groups correctly (an accuracy of

more than 99%). Nice!

## Building and Running

### Prerequisites

Before building and running the example, please ensure you have version 7 or later of the

[Java JDK](http://www.oracle.com/technetwork/java/javase/downloads/index.html) installed.

### Building

To build the example, first clone the GitHub repo at [https://github.com/nfergu/popstrat](https://github.com/nfergu/popstrat).

Then [download and install Maven](http://maven.apache.org/download.cgi). Then, at the command line, type:

```

mvn clean package

```

This will build a JAR (target/uber-popstrat-0.1-SNAPSHOT.jar), containing the `PopStrat` class,

as well as all of its dependencies.

### Running

First, [download Spark version 1.2.0](http://spark.apache.org/downloads.html) and unpack it on your machine.

Next you'll need to get some genomics data. Go to your

[nearest mirror of the 1000 genomes FTP site](http://www.1000genomes.org/data#DataAccess).

From the `release/20130502/` directory download

the `ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz` file and

the `integrated_call_samples_v3.20130502.ALL.panel` file. The first file file is the genotype data for chromosome 22,

and the second file is the panel file, which describes the population group for each sample in the genotype data.

Unzip the genotype data before continuing. This will require around 10GB of disk space.

To speed up execution and save disk space, you can convert the genotype VCF file to [ADAM](https://github.com/bigdatagenomics/adam)

format (using the ADAM `transform` command) if you wish. However,

this will take some time up-front. Both ADAM and VCF formats are supported.

Next, run the following command:

```

YOUR_SPARK_HOME/bin/spark-submit --class "com.neilferguson.PopStrat" --master local[6] --driver-memory 6G target/uber-popstrat-0.1-SNAPSHOT.jar <genotypesfile> <panelfile>

```

Replacing <genotypesfile> with the path to your genotype data file (ADAM or VCF), and <panelfile> with the panel file

from 1000 genomes.

This runs the example using a local (in-process) Spark master with 6 cores and 6GB of RAM. You can run against a different

Spark cluster by modifying the options in the above command line. See the

[Spark documentation](https://spark.apache.org/docs/1.2.0/submitting-applications.html) for further details.

Using the above data, the example may take up to 2-3 hours to run, depending on hardware. When it is finished, you should

see a [confusion matrix](http://en.wikipedia.org/wiki/Confusion_matrix) which shows the predicted versus the actual

populations. If all has gone well, this should show an accuracy of more than 99%.

See the "Code" section above for more details on what exactly you should expect to see.

## Conclusion

In this post, we have shown how to combine ADAM and Apache Spark with H2O's deep learning capabilities to predict

an individual's population group based on his or her genomic data. Our results demonstrate that we can predict these

very well, with more than 99% accuracy. Our choice of technologies makes for a relatively straightforward implementation,

and we expect it to be very scalable.

Future work could involve validating the scalability of our solution on more hardware, trying to predict a wider

range of population groups (currently we only predict 3 groups), and tuning the deep learning hyper-parameters to

achieve even better accuracy.