Apache Spark SQL & Machine Learning on Genetic Variant Classifications

Introduction

Covid-19 has driven up the interest of the bioinformatics, the data science of collecting and analyzing complex biological data such as genetic codes and sequences, a combination of mathematics and biology.

Context

ClinVar is a public resource containing annotations about human genetic variants. Detail can be found on the National Institute of Health website.
These variants are usually manually classified by clinical laboratories on a categorical spectrum ranging from benign, likely benign, uncertain significance, likely pathogenic, and pathogenic. Variants that have conflicting classifications (from laboratory to laboratory) can cause confusion when clinicians or researchers try to interpret whether the variant has an impact on the disease of a given patient.

Objective

The objective is to predict whether a ClinVar variant will have conflicting classifications. This is presented as a binary classification problem, each record in the dataset is a genetic variant.
Conflicting classifications are when two of any of the following three categories are present for one variant, two submissions of one category are not considered conflicting.
Β· Likely Benign or Benign
Β· VUS (Variance of Uncertain Significance. A variation in a genetic sequence for which the association with disease risk is unclear)
Β· Likely Pathogenic or Pathogenic

Class Label

Conflicting classification has been assigned to the CLASS column. It is a binary representation of whether or not a variant has conflicting classifications, where 0 represents consistent classifications and 1 represents conflicting classifications.
The genetic variant dataset is public domain available from Kaggle

ML Project Steps

Data Acquisition

Download the csv file from Kaggle, place into a folder on Apache Spark driver node:

Load the csv file into Spark SQL DataFrame

1
import org.apache.spark._
2
import org.apache.spark.SparkContext._
3
import org.apache.spark.rdd._
4
import org.apache.spark.util.LongAccumulator
5
import org.apache.log4j._
6
import scala.collection.mutable.ArrayBuffer
7
import org.apache.spark.sql._
8
Logger.getLogger("org").setLevel(Level.ERROR)
9
10
val spark = SparkSession
11
.builder
12
.appName("genetic classifier")
13
.master("local[*]")
14
.config("spark.sql.warehouse.dir", "file:///tmp")
15
.getOrCreate()
16
/*
17
Load the csv file into SparkSQL DataFrame, with options:
18
inferSchema=true, set the datatype of the feature columns based upon value data types
19
header=true, set the DataFrame column names from the headers of the csv file.
20
*/
21
val ds = spark.read.format("csv").option("inferSchema", "true").option("header", "true")
22
.option("quote", "\"")
23
.load("file:///home/bigdata2/dataset/clinvar_conflicting.csv")
24
val df: DataFrame = ds.toDF()
Copied!

Data Exploration:

1
df.show(3)
2
/*
3
|CHROM| POS|REF|ALT|AF_ESP|AF_EXAC|AF_TGP| CLNDISDB|CLNDISDBINCL| CLNDN|CLNDNINCL| CLNHGVS|CLNSIGINCL| CLNVC| CLNVI| MC|ORIGIN|SSR|CLASS|Allele| Consequence| IMPACT| SYMBOL|Feature_type| Feature| BIOTYPE|EXON|INTRON|cDNA_position|CDS_position|Protein_position|Amino_acids| Codons|DISTANCE|STRAND|BAM_EDIT| SIFT| PolyPhen|MOTIF_NAME|MOTIF_POS|HIGH_INF_POS|MOTIF_SCORE_CHANGE|LoFtool|CADD_PHRED| CADD_RAW|BLOSUM62|
4
| 1|1168180| G| C|0.0771| 0.1002|0.1066| MedGen:CN169374| nan| not_specified| nan|NC_000001.10:g.11...| nan|single_nucleotide...|UniProtKB_(protei...|SO:0001583|missen...| 1|nan| 0| C|missense_variant|MODERATE|B3GALT6| Transcript| NM_080605.3|protein_coding| 1/1| null| 552| 522| 174| E/D|gaG/gaC| null| 1| null| tolerated| benign| null| null| null| null| null| 1.053|-0.208682| 2|
5
| 1|1470752| G| A| 0.0| 0.0| 0.0|MedGen:C1843891,O...| nan|Spinocerebellar_a...| nan|NC_000001.10:g.14...| nan|single_nucleotide...|OMIM_Allelic_Vari...|SO:0001583|missen...| 1|nan| 0| A|missense_variant|MODERATE|TMEM240| Transcript|NM_001114748.1|protein_coding| 4/4| null| 523| 509| 170| P/L|cCg/cTg| null| -1| OK|deleterious_low_c...| benign| null| null| null| null| null| 31.0| 6.517838| -3|
6
| 1|1737942| A| G| 0.0| 1.0E-5| 0.0|Human_Phenotype_O...| nan|Strabismus|Nystag...| nan|NC_000001.10:g.17...| nan|single_nucleotide...|OMIM_Allelic_Vari...|SO:0001583|missen...| 35|nan| 1| G|missense_variant|MODERATE| GNB1| Transcript| NM_002074.4|protein_coding|6/12| null| 632| 239| 80| I/T|aTc/aCc| null| -1| OK| deleterious|probably_damaging| null| null| null| null| null| 28.1| 6.061752| -1|
7
8
only showing top 3 rows
9
*/
Copied!
Notice there are lots missing values, i.e., null values in the dataset. Dataset having significant missing values usually affect training and prediction accuracy.
Some column values are large, like 6 digits, some column values are smaller than 1. Without addressing, this can confuse ML algorithm, which may regard some small values as statistically insignificant and disregard as result. This can be addressed in the data preprocessing using scaler. I plan to use MinMaxScaler, which will make all feature column values between 0 and 1

Evaluate schema:

1
df.printSchema
2
​
3
/*
4
root
5
|-- CHROM: string (nullable = true)
6
|-- POS: integer (nullable = true)
7
|-- REF: string (nullable = true)
8
|-- ALT: string (nullable = true)
9
|-- AF_ESP: double (nullable = true)
10
|-- AF_EXAC: double (nullable = true)
11
|-- AF_TGP: double (nullable = true)
12
|-- CLNDISDB: string (nullable = true)
13
|-- CLNDISDBINCL: string (nullable = true)
14
|-- CLNDN: string (nullable = true)
15
|-- CLNDNINCL: string (nullable = true)
16
|-- CLNHGVS: string (nullable = true)
17
|-- CLNSIGINCL: string (nullable = true)
18
|-- CLNVC: string (nullable = true)
19
|-- CLNVI: string (nullable = true)
20
|-- MC: string (nullable = true)
21
|-- ORIGIN: integer (nullable = true)
22
|-- SSR: string (nullable = true)
23
|-- CLASS: integer (nullable = true)
24
|-- Allele: string (nullable = true)
25
|-- Consequence: string (nullable = true)
26
|-- IMPACT: string (nullable = true)
27
|-- SYMBOL: string (nullable = true)
28
|-- Feature_type: string (nullable = true)
29
|-- Feature: string (nullable = true)
30
|-- BIOTYPE: string (nullable = true)
31
|-- EXON: string (nullable = true)
32
|-- INTRON: string (nullable = true)
33
|-- cDNA_position: string (nullable = true)
34
|-- CDS_position: string (nullable = true)
35
|-- Protein_position: string (nullable = true)
36
|-- Amino_acids: string (nullable = true)
37
|-- Codons: string (nullable = true)
38
|-- DISTANCE: integer (nullable = true)
39
|-- STRAND: integer (nullable = true)
40
|-- BAM_EDIT: string (nullable = true)
41
|-- SIFT: string (nullable = true)
42
|-- PolyPhen: string (nullable = true)
43
|-- MOTIF_NAME: string (nullable = true)
44
|-- MOTIF_POS: integer (nullable = true)
45
|-- HIGH_INF_POS: string (nullable = true)
46
|-- MOTIF_SCORE_CHANGE: double (nullable = true)
47
|-- LoFtool: double (nullable = true)
48
|-- CADD_PHRED: double (nullable = true)
49
|-- CADD_RAW: double (nullable = true)
50
|-- BLOSUM62: integer (nullable = true)
51
52
*/
Copied!
The DataFrame has 46 columns, with one column CLASS as target/label and 45 columns as features.
Some of them are String data types, which need to be encoded into numeric datatype for ML algorithm to crunching them.
Distribution of row count of CLASS 1 (genetically conflicting) and 0 (genetically consistent) as following:
1
df.groupBy("CLASS").count().show()
2
​
3
/*
4
+-----+-----+
5
|CLASS|count|
6
+-----+-----+
7
| 1|16434|
8
| 0|48754|
9
+-----+-----+
10
*/
Copied!

Data Preprocessing

NULL replacement

Replace null values in the String columns as β€œABCD” and 0 in numeric columns

Encode the String value to integer using StringIndexer

Rearrange the DataFrame, after encoding, the first column of the resultant DataFrame is CLASS, the label.
1
import org.apache.spark.ml.feature.StringIndexer
2
var df1=df.select("CLASS")
3
val featureDF=df.drop("CLASS")
4
val typeArray=featureDF.dtypes
5
for (i<-0 until featureDF.columns.size)
6
{
7
if (typeArray(i)._2=="StringType")
8
{
9
var indexer = new StringIndexer().setInputCol(typeArray(i)._1).setOutputCol(typeArray(i)._1+"Index")
10
var indexed = indexer.fit(featureDF.na.fill("ABCD").select(typeArray(i)._1)).transform(featureDF.na.fill("ABCD").select(typeArray(i)._1))
11
var temp=indexed.select(typeArray(i)._1+"index")
12
df1=df1.withColumn("id", monotonically_increasing_id())
13
.join(temp.withColumn("id", monotonically_increasing_id()), Seq("id"))
14
}
15
else
16
{
17
var temp=featureDF.na.fill(0).select(typeArray(i)._1)
18
df1=df1.withColumn("id", monotonically_increasing_id()).join(temp.withColumn("id", monotonically_increasing_id()), Seq("id"))
19
}
20
}
21
df1=df1.drop("id")
22
df1.printSchema
23
/*
24
​
25
root
26
|-- CLASS: integer (nullable = true)
27
|-- CHROMindex: double (nullable = false)
28
|-- POS: integer (nullable = false)
29
|-- REFindex: double (nullable = false)
30
|-- ALTindex: double (nullable = false)
31
|-- AF_ESP: double (nullable = false)
32
|-- AF_EXAC: double (nullable = false)
33
|-- AF_TGP: double (nullable = false)
34
|-- CLNDISDBindex: double (nullable = false)
35
|-- CLNDISDBINCLindex: double (nullable = false)
36
|-- CLNDNindex: double (nullable = false)
37
|-- CLNDNINCLindex: double (nullable = false)
38
|-- CLNHGVSindex: double (nullable = false)
39
|-- CLNSIGINCLindex: double (nullable = false)
40
|-- CLNVCindex: double (nullable = false)
41
|-- CLNVIindex: double (nullable = false)
42
|-- MCindex: double (nullable = false)
43
|-- ORIGIN: integer (nullable = false)
44
|-- SSRindex: double (nullable = false)
45
|-- Alleleindex: double (nullable = false)
46
|-- Consequenceindex: double (nullable = false)
47
|-- IMPACTindex: double (nullable = false)
48
|-- SYMBOLindex: double (nullable = false)
49
|-- Feature_typeindex: double (nullable = false)
50
|-- Featureindex: double (nullable = false)
51
|-- BIOTYPEindex: double (nullable = false)
52
|-- EXONindex: double (nullable = false)
53
|-- INTRONindex: double (nullable = false)
54
|-- cDNA_positionindex: double (nullable = false)
55
|-- CDS_positionindex: double (nullable = false)
56
|-- Protein_positionindex: double (nullable = false)
57
|-- Amino_acidsindex: double (nullable = false)
58
|-- Codonsindex: double (nullable = false)
59
|-- DISTANCE: integer (nullable = false)
60
|-- STRAND: integer (nullable = false)
61
|-- BAM_EDITindex: double (nullable = false)
62
|-- SIFTindex: double (nullable = false)
63
|-- PolyPhenindex: double (nullable = false)
64
|-- MOTIF_NAMEindex: double (nullable = false)
65
|-- MOTIF_POS: integer (nullable = false)
66
|-- HIGH_INF_POSindex: double (nullable = false)
67
|-- MOTIF_SCORE_CHANGE: double (nullable = false)
68
|-- LoFtool: double (nullable = false)
69
|-- CADD_PHRED: double (nullable = false)
70
|-- CADD_RAW: double (nullable = false)
71
|-- BLOSUM62: integer (nullable = false)
72
73
*/
Copied!

Make all integer column to double

1
import org.apache.spark.sql.types._
2
val newDf = df1.select(df1.columns.map(c => col(c).cast(DoubleType)) : _*)
3
newDf.printSchema
4
/*
5
root
6
|-- CLASS: double (nullable = true)
7
|-- CHROMindex: double (nullable = false)
8
|-- POS: double (nullable = false)
9
|-- REFindex: double (nullable = false)
10
|-- ALTindex: double (nullable = false)
11
|-- AF_ESP: double (nullable = false)
12
|-- AF_EXAC: double (nullable = false)
13
|-- AF_TGP: double (nullable = false)
14
|-- CLNDISDBindex: double (nullable = false)
15
|-- CLNDISDBINCLindex: double (nullable = false)
16
|-- CLNDNindex: double (nullable = false)
17
|-- CLNDNINCLindex: double (nullable = false)
18
|-- CLNHGVSindex: double (nullable = false)
19
|-- CLNSIGINCLindex: double (nullable = false)
20
|-- CLNVCindex: double (nullable = false)
21
|-- CLNVIindex: double (nullable = false)
22
|-- MCindex: double (nullable = false)
23
|-- ORIGIN: double (nullable = false)
24
|-- SSRindex: double (nullable = false)
25
|-- Alleleindex: double (nullable = false)
26
|-- Consequenceindex: double (nullable = false)
27
|-- IMPACTindex: double (nullable = false)
28
|-- SYMBOLindex: double (nullable = false)
29
|-- Feature_typeindex: double (nullable = false)
30
|-- Featureindex: double (nullable = false)
31
|-- BIOTYPEindex: double (nullable = false)
32
|-- EXONindex: double (nullable = false)
33
|-- INTRONindex: double (nullable = false)
34
|-- cDNA_positionindex: double (nullable = false)
35
|-- CDS_positionindex: double (nullable = false)
36
|-- Protein_positionindex: double (nullable = false)
37
|-- Amino_acidsindex: double (nullable = false)
38
|-- Codonsindex: double (nullable = false)
39
|-- DISTANCE: double (nullable = false)
40
|-- STRAND: double (nullable = false)
41
|-- BAM_EDITindex: double (nullable = false)
42
|-- SIFTindex: double (nullable = false)
43
|-- PolyPhenindex: double (nullable = false)
44
|-- MOTIF_NAMEindex: double (nullable = false)
45
|-- MOTIF_POS: double (nullable = false)
46
|-- HIGH_INF_POSindex: double (nullable = false)
47
|-- MOTIF_SCORE_CHANGE: double (nullable = false)
48
|-- LoFtool: double (nullable = false)
49
|-- CADD_PHRED: double (nullable = false)
50
|-- CADD_RAW: double (nullable = false)
51
|-- BLOSUM62: double (nullable = false)
52
53
*/
Copied!

Feature Vectorization

Apache Spark ML requires to place feature data to be in form of vector, i.e., place all feature columns into one vector, following code is to generate a column that is a vector that contains all feature columns
1
import org.apache.spark.ml.feature.VectorAssembler
2
import org.apache.spark.ml.linalg.Vectors
3
val assembler = new VectorAssembler()
4
.setInputCols(Array("CHROMindex", "POS", "REFindex", "ALTindex", "AF_ESP", "AF_EXAC", "AF_TGP", "CLNDISDBindex"
5
, "CLNDISDBINCLindex", "CLNDNindex", "CLNDNINCLindex", "CLNHGVSindex", "CLNSIGINCLindex", "CLNVCindex"
6
, "CLNVIindex", "MCindex", "ORIGIN", "SSRindex", "Alleleindex", "Consequenceindex", "IMPACTindex"
7
, "SYMBOLindex", "Feature_typeindex", "Featureindex", "BIOTYPEindex", "EXONindex", "INTRONindex"
8
, "cDNA_positionindex", "CDS_positionindex", "Protein_positionindex", "Amino_acidsindex", "Codonsindex"
9
, "DISTANCE", "STRAND", "BAM_EDITindex", "SIFTindex", "PolyPhenindex", "MOTIF_NAMEindex", "MOTIF_POS"
10
, "HIGH_INF_POSindex", "MOTIF_SCORE_CHANGE", "LoFtool", "CADD_PHRED", "CADD_RAW", "BLOSUM62"))
11
.setOutputCol("features")
12
val output = assembler.transform(newDf)
13
Resultant output DataFrame now has a new column features that is the vector containing all feature columns.
14
output.select("features").show(2,false)
15
/*
16
17
|features
18
|(45,[0,1,2,3,4,5,6,11,14,16,18,21,23,25,27,28,29,30,31,33,35,36,42,43,44],[3.0,1168180.0,1.0,3.0,0.0771,0.1002,0.1066,326.0,27451.0,1.0,3.0,1614.0,1728.0,15.0,618.0,471.0,304.0,56.0,129.0,1.0,2.0,1.0,1.053,-0.208682,2.0])|
19
|(45,[0,1,2,3,7,9,11,14,16,18,21,23,25,27,28,29,30,31,33,34,35,36,42,43,44],[3.0,1470752.0,1.0,1.0,7922.0,8822.0,467.0,25430.0,1.0,1.0,2293.0,2054.0,9.0,1534.0,1465.0,116.0,13.0,17.0,-1.0,1.0,4.0,1.0,31.0,6.517838,-3.0]) |
20
​
21
*/
Copied!
Extract from output DataFrame into a new DataFrame, that only contains vector column features and label column CLASS
1
val whole=output.select("features","CLASS")
2
whole.printSchema
3
/*
4
root
5
|-- features: vector (nullable = true)
6
|-- CLASS: double (nullable = true)
7
8
*/
Copied!

Normalize feature values in vector column features, using MinMax Scaler

1
import org.apache.spark.ml.feature.MinMaxScaler
2
import org.apache.spark.ml.linalg.Vectors
3
val scaler = new MinMaxScaler()
4
.setInputCol("features")
5
.setOutputCol("scaledFeatures")
6
// Compute summary statistics and generate MinMaxScalerModel
7
val scalerModel = scaler.fit(whole)
8
// rescale each feature to range [min, max].
9
val scaledData = scalerModel.transform(whole)
10
println(s"Features scaled to range: [${scaler.getMin}, ${scaler.getMax}]")
11
scaledData.select("features", "scaledFeatures").show(2,false)
12
/*
13
Features scaled to range: [0.0, 1.0]
14
|features |scaledFeatures
15
|(45,[0,1,2,3,4,5,6,11,14,16,18,21,23,25,27,28,29,30,31,33,35,36,42,43,44],[3.0,1168180.0,1.0,3.0,0.0771,0.1002,0.1066,326.0,27451.0,1.0,3.0,1614.0,1728.0,15.0,618.0,471.0,304.0,56.0,129.0,1.0,2.0,1.0,1.053,-0.208682,2.0])|(45,[0,1,2,3,4,5,6,11,14,16,18,21,23,25,27,28,29,30,31,33,35,36,40,42,43,44],[0.13043478260869565,0.004713998164155383,0.0011560693641618498,0.006564551422319475,0.15450901803607214,0.2004440977014943,0.21328531412565024,0.005000997131329867,0.9926592897953279,0.001949317738791423,0.00804289544235925,0.6932989690721649,0.729421696918531,0.004595588235294118,0.044237652111667865,0.03447266339749689,0.041422537130399235,0.044374009508716325,0.05810810810810811,1.0,0.5,0.25,0.9999999999999999,0.010636363636363637,0.10125579884341004,0.8333333333333333])|
16
|(45,[0,1,2,3,7,9,11,14,16,18,21,23,25,27,28,29,30,31,33,34,35,36,42,43,44],[3.0,1470752.0,1.0,1.0,7922.0,8822.0,467.0,25430.0,1.0,1.0,2293.0,2054.0,9.0,1534.0,1465.0,116.0,13.0,17.0,-1.0,1.0,4.0,1.0,31.0,6.517838,-3.0]) |(45,[0,1,2,3,7,9,11,14,16,18,21,23,25,27,28,29,30,31,34,35,36,40,42,43],[0.13043478260869565,0.0059359829438109775,0.0011560693641618498,0.002188183807439825,0.8580093144156828,0.9528026784749973,0.007164005093040024,0.9195776379547261,0.001949317738791423,0.002680965147453083,0.9849656357388316,0.8670325031658928,0.0027573529411764703,0.1098067287043665,0.10722388933616336,0.015805968115547075,0.010301109350237718,0.0076576576576576575,0.5,1.0,0.25,0.9999999999999999,0.31313131313131315,0.2305282934974466])
17
​
18
*/
Copied!

Train/Test Split, randomly split 70% of rows for training, 30% for test

1
// Prepare training and test data.
2
val wholeScaledData=scaledData.select("scaledFeatures","CLASS")
3
val Array(training, test) = wholeScaledData.randomSplit(Array(0.7, 0.3), seed = 12345)
Copied!

Training and Prediction

Algorithm Selection

Since this is a binary classification, I plan to use Multilayer Perceptron Classifier and Logistic Regression.

Multilayer Perceptron Classifier, there are 45 feature columns, and binary classification, I construct a neural network of 4 layer as below:

1
val layers = Array[Int](45, 45, 45, 2)
2
import org.apache.spark.ml.classification.MultilayerPerceptronClassifier
3
import org.apache.spark.ml.evaluation.MulticlassClassificationEvaluator
4
val trainer = new MultilayerPerceptronClassifier()
5
.setLayers(layers)
6
.setBlockSize(128)
7
.setSeed(1234L)
8
.setMaxIter(100)
9
.setFeaturesCol("scaledFeatures")
10
.setLabelCol("CLASS")
Copied!

Train and test the neural network

1
val model=trainer.fit(training.select("scaledFeatures","CLASS"))
2
//and Test
3
val result = model.transform(test)
4
//Evaluate accuracy metrics
5
val predictionAndLabels = result.select("prediction", "CLASS")
6
val evaluator = new MulticlassClassificationEvaluator().setMetricName("accuracy").setLabelCol("CLASS")
7
println(s"Test set accuracy = ${evaluator.evaluate(predictionAndLabels)}")
8
//Test set accuracy = 0.7482470955524848
Copied!

Logistic Regression Classifier

1
import org.apache.spark.ml.{Pipeline, PipelineStage}
2
import org.apache.spark.ml.classification.{LogisticRegression, LogisticRegressionModel}
3
import scala.collection.mutable
4
val maxIter=100
5
val lr = new LogisticRegression()
6
.setFeaturesCol("scaledFeatures")
7
.setLabelCol("CLASS")
8
.setMaxIter(maxIter)
9
val model_lr = lr.fit(training)
10
val prediction_lr = model_lr.transform(test)
11
import org.apache.spark.ml.evaluation.MulticlassClassificationEvaluator
12
val my_mc_accu = new MulticlassClassificationEvaluator().setPredictionCol("prediction").setLabelCol("CLASS").setMetricName("accuracy")
13
my_mc_accu.evaluate(prediction_lr)
14
my_mc_accu: org.apache.spark.ml.evaluation.MulticlassClassificationEvaluator = mcEval_74ebac4de5a0
15
//res29: Double = 0.7463534469522494
Copied!

Take away:

Both Mutilayer Perceptron Classifier and Logistic Regression produce same prediction accuracy about 75% given the dataset having significant amount missing values (null values). To improve machine learning performance, reduce missing values and possible adding additional feature columns may help.
As always, code used in this writing is in my GitHub Repo:
Thank you for your time viewing this writing.
Last modified 1yr ago