Genomic information is increasingly being used for medical research, giving rise to the need for efficient analysis methodology able to cope with thousands of individuals and millions of variants. Catering for this need, we developed VariantSpark, a framework for applying machine learning algorithms in MLlib to genomic variant data using the efficient in-memory Spark compute engine. We demonstrate a speedup compared to our earlier Hadoop-based implementation as well as a published Spark-based approach using the ADAM framework. Adapting emerging methodologies for fast efficient analysis of large diverse data volumes to process genomic information, will be the cornerstone for precision genome medicine, which promises tailored healthcare for diagnosis and treatment.
Introduction
We apply Apache Spark to publicly available genomic data. We demonstrate the
potential Spark has to process huge amounts of data, as well as its scalability and
flexibility. We also compare it to our previous implementation of similar techniques in
Apache Hadoop, in regards to speed and performance.
1. Background
With the decreasing costs of DNA sequencing [
‘Apache Spark’ [
Recognising this, the Big Data Genomics (BDG) group has recently demonstrated
the strength of Spark in a genomic clustering application using ADAM, a set of formats
and APIs as well as processing stage implementations for genomic data [
We implemented VariantSpark in Scala, leveraging Apache Spark and MLLib. The aim
of the first step is to read in VCF files and effectively transpose them to vectors of each
individual’s variants. Using tuples to store the data in key-value pairs allows us to
achieve this task in a distributed fashion that should scale to allow for an input of
hundreds of gigabytes. Effectively, we split each line from the VCF file(s) into an array.
We then zip each item in each array with its heading from the VCF file, resulting in an
array of tuples (See Figure 1 Spark “Zip with header”). Each of these arrays now stores
a single variant for every individual in the dataset and can be evenly distributed across
the cluster in Spark’s Resilient Distributed Dataset (RDD) format [
We therefore zip each array with a unique index representing a variant (See Figure 1 Spark “Zip with index”). We now have a tuple containing an integer paired with an array of yet more tuples. We can now use the FLATMAP function to effectively flatten the data structure and get closer to the structure we require. The FLATMAP, with a custom inner function, results in a tuple for every variant across every individual. This tuple contains the individual’s ID as the ‘key’, and as the ‘value’, a secondary tuple of the variant index and the variant (See Figure 1 Spark “Flatmap”).
The value of each variant is the distance from the wild-type allele. For example, a homozygous variant is 2, a heterozygous variant is 1 and no variant is 0. We calculate this value from the variant string during the aforementioned FLATMAP stage. In addition, to keep data-sizes to a minimum, we remove any variants with a value of zero. Therefore, we can store the data as a sparse vector, rather than a dense vector. The sparse vectors store non-zero variants with their index, rather than the entire array of variants, as a dense-vector would do.
We need to group the collection of tuples by the individual using Spark’s GROUPBYKEY. Following this function, we have an array for each individual, containing the variant-tuples for said individual (See Figure 1 Spark “Group by individual”). This is essentially the structure of a sparse vector, and we can map each array to the sparse vector structure required by MLlib’s machine-learning algorithms.
We proceed to create an MLlib k-means model and train the model on this data. We report the accuracy based on the resulting cluster-set compared to the expected result, according to each individual’s population and super population. We run a maximum of ten k-means iterations and use the Manhattan distance measure.
3.1. Comparison to Hadoop
We previously implemented a similar program in Hadoop using the Apache Mahout
machine-learning library. Although falling shy of Apache’s claim, we observe a
speedup of more than 50% with our Spark implementation. In Hadoop, clustering
under 1 million variants took close to 16 minutes. Processing and clustering over 1
million variants using Spark, on the other hand, took less than 7 minutes. One
contributing factor to the speedup is in-memory caching. Unlike Hadoop, which writes
the output of its data transformations to disk, Spark can cache intermediate
datastructures in memory. This minimises disk IO and removes a huge bottleneck. Another
reason for the better performance is the different programming model. Hadoop relies
on the MapReduce paradigm, which, although easy to use with distributed systems [
We processed and clustered the human chromosome 22, which includes close to 500,000 variants, using the same options as BDG. The process took approximately 1.3 minutes to complete. This includes training and predicting, and is faster than the 2 minutes runtime reported by BDG. Our method also eliminates the time required for converting VCF files to the ADAM format. However, a trade-off for our faster approach is disk-space. Because the ADAM format uses compressed files, as opposed to uncompressed VCF files, their file-sizes would be smaller. 3.3. Scalability We demonstrate the potential of Spark to scale from a subset of chromosome 1 containing less than one hundred thousand variants, to multiple chromosomes containing millions of variants. The time required scales in a linear-fashion, based on the number of variants present in the data. However, the overhead of Spark becomes apparent with relatively small datasets, as can be seen in Table 2, where the times required for completing the two smallest jobs only differ by one second. For the smaller jobs, Spark requires only 1GB memory per executor (or core). Increasing the number of variants beyond 2 million requires a greater memory allocation for each executor. Although this is not a problem with high-memory clusters, on smaller clusters, such as our four-node Azure cluster, increasing the memory allocation beyond 1GB will limit the number of executors that can run simultaneously, due to memory-constraints. Although of no advantage to clustering quality, as we explain in the next section, we further demonstrate the scalability by increasing the dataset to include multiple chromosomes. These jobs complete successfully, albeit after more time and with higher memory requirements (Table 1).
Spark not only scales well with data size, but also with cluster size. To
demonstrate this we run identical jobs on our in-house cluster. Clustering every
chromosome 1 variant from VCF files takes only six minutes when using 80 executers.
As is expected, this is much faster than the small cluster which took over 17 minutes
for the same job. Although these two clusters are different in OS, size and hardware
configuration, a Spark jobs, being implemented in Scala (or Java), can run on any
cluster. The only change we make is to the launch script where we specify how many
executors to use.
3.4. Clustering quality
When we cluster the entire group of 1092 individuals, the majority are clustered with
individuals from the same super-population. This is portrayed in the Adjusted Rand
Index (ARI) [
We can substantially improve the accuracy to a perfect classification (ARI=1.0) by removing the fourth super-population, AMR (American). Including AMR individuals, we observe that the majority of these individuals are placed in the same group as Europeans, likely accurately reflecting their migrational backgrounds. Only a minority of AMR individuals form an independent group, likely comprising of genetic information otherwise not captured by the 26 sub-populations of the 1000 genomes project.
Interestingly, we achieved these results with as few as one million variants, i.e. less than half the variants available in chromosome 1, and increasing the size of the dataset further saw no improvement on clustering the the AMR individuals.
We demonstrate the potential applications of processing VCF files using Apache Spark. We successfully apply k-means clustering to over 1 thousand individuals with millions of variants. We explore performance on different architectures starting with a relatively small 24-node cluster and subsequently scale up to a cluster with hundreds of nodes. We demonstrate our method of transforming and clustering individuals from VCF files to be faster than clustering individuals stored in the ADAM format, albeit at the cost of disk space.
In future work we aim to improve the disk space footprint by using compressed VCF files. The 1000 genomes consortium distributes their VCF files with gzip compression, which is unsuitable for use in Spark as it is not possible to split gzip files. We therefore plan to investigate other compression methods, such as bz2, which Spark can read line by line. We expect this to not only minimise storage requirements of the actual data, but also further reduce runtime, making the approach attractive for applying the more complicated machine learning approaches to genomic data.