---
title: D. Machine Learning
author:
Martin Morgan (mtmorgan@fredhutch.org)<br />
Sonali Arora (sarora@fredhutch.org)
output:
BiocStyle::html_document:
toc: true
vignette: >
%\VignetteIndexEntry{D. Machine Learning}
%\VignetteEngine{knitr::rmarkdown}
\usepackage[utf8]{inputenc}
---
```{r style, echo=FALSE, results='asis'}
BiocStyle::markdown()
```
## Introduction to Machine Learning
Lets say that we are interested in predicting the run time of an athlete
depending on his shoe size, height and weight in a study of 100 people.
We can do so using a simple linear regression model where
```{r eval=FALSE}
y = beta0 + beta1 * height + beta2 * weight + beta3 * shoe_size
```
Here y is the response variable (run time), n is the number of observations
(100 people), p is the number of variables/ features/ predictors (3 IE height,
weight, shoe size), X is a nxp matrix
This data set is a low dimensional data where n >> p but most of the biological
data sets coming out of modern biological techniques are high dimensional IE
n << p This poses statistical challenge and simple linear regression can no
longer help us.
For example,
* Identify the risk factors(genes) for prostrate cancer based on gene
expression data
* Predict the chances of breast cancer survival in a patient.
* Identify patterns of gene expression among different sub types of
breast cancer
In all of the 3 examples, listed above n, number of observations, is 30-40 patients
whereas p, number of features, is approximately 30,000 genes. Try writing a linear
regression formula for the outcome variable, y, in any of the above three
scenarios..
Listed below are things that can go wrong with high dimensional data
- some of these predictors are useful, some are not
- if we include too many predictors, we can over fit the data
This is why we need Machine Learning. Lets first introduce some basic concepts
and then dive into examples and a lab session.
**Supervised Learning** - Use a data set X to predict the association with a
response variable Y. The response variable can be continuous or categorical.
For example: Predicting the chances of breast cancer survival in a patient.
**Unsupervised Learning** - Discover the associations or patterns in X. No
response variable is present. For example: Cluster similar genes into groups.
**Training & Test Datasets** - Usually we split observation into test and
training data sets. We fit the model on the training data set and evaluate on the
test data set. The test set error rate is an estimate of the models performance
on future data sets.
**Model Selection** - We usually consider numerous models for a given problem.
For example, we are trying to identify the genes responsible for a given
disease using gene expression data set- we could have the following models
a) model 1 - Use all 30000 genes from the array to build a model
b) model 2 - we include only genes related to the pathway that we know is
upregulated in that disease to build a model
c) model 3 - include genes found in literature which are known to influence
this disease
It is highly recommended to use the test set only on our final model to see
how our model will do with new, unseen data. So how do we pick the best
model which can be tested on the test data set?
**Cross-validation**
We can use different approaches to find the best model. Lets look at the
commonly used approaches, namely, validation set, leave one out
cross-validation, k-fold cross validation.
Briefly, the __validation set approach__ deals with diving the full data sets into
3 groups - training set, validation set and the test set. We train the models on
the training set, evaluate their performance on the validation set and then the
best model is chosen to fit on the test set.
The __leave one out cross validation__ starts with fitting n models (where n is
number of observations in the training data set), each on n-1 observations,
evaluating each model on the left-out observation. The best model is the one
for which the total test error is the smallest and that is then used to predict
the test set.
Lastly the __5 fold cross validation__ (here k=5), is splitting the training
data set into 5 sets and repeatedly training the model on the other 4 sets and
evaluating the performance on the fifth.
**Bias, Variance, Overfitting** - Bias refers to the average difference between
the actual betas and the predicted betas, Variance refers to the amount by
which the betas differ across experiments. As the model complexity(no of
variables) increases, the bias decreases and the variance increases. This is
know as the Bias-Variance Tradeoff and a model that has too much of variance,
is said to be over fit.
## Datasets
For **Unsupervised learning**, we will use RNA-Seq count data from the
Biocoductor package, `r Biocpkg("airway")`. From the abstract, a brief
description of the RNA-Seq experiment on airway smooth muscle (ASM) cell
lines: “Using RNA-Seq, a high-throughput sequencing method, we characterized
transcriptomic changes in four primary human ASM cell lines that were treated
with dexamethasone - a potent synthetic glucocorticoid (1 micromolar for
18 hours).â€
```{r message=FALSE}
library(airway)
data("airway")
se <- airway
colData(se)
library("DESeq2")
dds <- DESeqDataSet(se, design = ~ cell + dex)
```
For **Supervised learning**, we will use cervical count data from the
Biocoductor package, `r Biocpkg("MLSeq")`. This data set contains
expressions of 714 miRNA's of human samples. There are 29 tumor and 29
non-tumor cervical samples. For learning purposes, we can treat these
as two separate groups and run various classification algorithms.
```{r message=FALSE}
library(MLSeq)
filepath = system.file("extdata/cervical.txt", package = "MLSeq")
cervical = read.table(filepath, header = TRUE)
```
## Unsupervised Learning
Unsupervised Learning is a set of statistical tools intended for the setting
in which we have only a set of 'p' features measured on 'n' observations.
We are primarily interested in discovering interesting
things about the 'p' features.
Unsupervised Learning is often performed as a part of Exploratory Data Analysis.
These tools help us to get a good idea about the data set. Unlike a supervised
learning problem, where we can use prediction to gain some confidence about our
learning algorithm, there is no way to check our model. The learning algorithm
is thus, aptly named "unsupervised".
**RLOG TRANSFORMATION**
Many common statistical methods for exploratory analysis of multidimensional
data, especially methods for clustering and ordination (e.g.,
principal-component analysis and the like), work best for (at least
approximately) homoskedastic data; this means that the variance of an observed
quantity (here, the expression strength of a gene) does not depend on the mean.
In RNA-Seq data, the variance grows with the mean.If one performs PCA
(principal components analysis) directly on a matrix of normalized read counts,
the result typically depends only on the few most strongly expressed genes
because they show the largest absolute differences between samples.
As a solution, DESeq2 offers the regularized-logarithm transformation, or rlog
for short. See the help for ?rlog for more information and options.
The function rlog returns a SummarizedExperiment object which contains the
rlog-transformed values in its assay slot:
```{r}
rld <- rlog(dds)
head(assay(rld))
```
To assess overall similarity between samples: Which samples are similar to each
other, which are different? Does this fit to the expectation from the
experiment's design? We use the R function dist to calculate the Euclidean
distance between samples. To avoid that the distance measure is dominated by
a few highly variable genes, and have a roughly equal contribution from all
genes, we use it on the rlog-transformed data
```{r}
sampleDists <- dist( t( assay(rld) ) )
sampleDists
```
Note the use of the function t to transpose the data matrix. We need this
because dist calculates distances between data rows and our samples constitute
the columns.
**HEATMAP**
We visualize the sample-to-sample distances in a heatmap, using the
function heatmap.2 from the gplots package. Note that we have changed the row
names of the distance matrix to contain treatment type and patient number
instead of sample ID, so that we have all this information in view when
looking at the heatmap.
```{r message=FALSE}
library("gplots")
library("RColorBrewer")
sampleDistMatrix <- as.matrix( sampleDists )
rownames(sampleDistMatrix) <- paste( rld$dex, rld$cell, sep="-" )
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
hc <- hclust(sampleDists)
heatmap.2( sampleDistMatrix, Rowv=as.dendrogram(hc),
symm=TRUE, trace="none", col=colors,
margins=c(2,10), labCol=FALSE )
```
**PCA**
Another way to visualize sample-to-sample distances is a principal-components
analysis (PCA). In this ordination method, the data points (i.e., here, the
samples) are projected onto the 2D plane such that they spread out in the two
directions which explain most of the differences in the data. The x-axis is
the direction (or principal component) which separates the data points the most.
The amount of the total variance which is contained in the direction is
printed in the axis label. Here, we have used the function plotPCA which comes
with DESeq2. The two terms specified by intgroup are the interesting groups
for labelling the samples; they tell the function to use them to choose colors.
```{r}
plotPCA(rld, intgroup = c("dex", "cell"))
```
From both visualizations, we see that the differences between cells are
considerable, though not stronger than the differences due to treatment
with dexamethasone. This shows why it will be important to account for this
in differential testing by using a paired design (“pairedâ€, because each dex
treated sample is paired with one untreated sample from the same cell line).
We are already set up for this by using the design formula ~ cell + dex when
setting up the data object in the beginning.
**MDS**
Another plot, very similar to the PCA plot, can be made using the
multidimensional scaling (MDS) function in base R. This is useful when we don't
have the original data, but only a matrix of distances. Here we have the MDS
plot for the distances calculated from the rlog transformed counts:
```{r}
library(ggplot2)
mds <- data.frame(cmdscale(sampleDistMatrix))
mds <- cbind(mds, colData(rld))
qplot(X1,X2,color=dex,shape=cell,data=as.data.frame(mds))
```
### Exercise:
Use the plotMDS function from the limma package to make a simila plot.
What is the advtange of using this function over base R's cmdscale?
**Solutions:**
A similar plot can be made using the plotMDS() function in limma where the input
is a matrix of log-fold expression values. Here the advantage is that the
distances on plot are proportional to log2-fold change and not only is the plot
created, but the object (with distance matrix) is also returned.
```{r plotMDS}
suppressPackageStartupMessages({
library(limma)
library(DESeq2)
library(airway)
})
plotMDS(assay(rld), col=as.integer(dds$dex), pch=as.integer(dds$cell))
```
## Supervised Learning
In supervised learning, along with the 'p' features, we
also have the a response Y measured on the same n observations. The goal is then
to predict Y using X (n x p matrix) for new observations.
For the cervical data, we know that the first 29 are non-Tumor samples
whereas the last 29 are Tumor samples. We will code these as 0 and 1
respectively. We will randomly sample 30% of our data and use that as a
test set. The remaining 70% of the data will be used as training data
```{r }
set.seed(9)
class = data.frame(condition = factor(rep(c(0, 1), c(29, 29))))
nTest = ceiling(ncol(cervical) * 0.2)
ind = sample(ncol(cervical), nTest, FALSE)
cervical.train = cervical[, -ind]
cervical.train = as.matrix(cervical.train + 1)
classtr = data.frame(condition = class[-ind, ])
cervical.test = cervical[, ind]
cervical.test = as.matrix(cervical.test + 1)
classts = data.frame(condition = class[ind, ])
```
MLSeq aims to make computation less complicated for a user and
allows one to learn a model using various classifier's with one single function.
The main function of this package is classify which requires data in the form of
a DESeqDataSet instance. The DESeqDataSet is a subclass of SummarizedExperiment,
used to store the input values, intermediate calculations and results of an
analysis of differential expression.
So lets create DESeqDataSet object for both the training and test set, and run
DESeq on it.
```{r}
cervical.trainS4 = DESeqDataSetFromMatrix(countData = cervical.train,
colData = classtr, formula(~condition))
cervical.trainS4 = DESeq(cervical.trainS4, fitType = "local")
cervical.testS4 = DESeqDataSetFromMatrix(countData = cervical.test, colData = classts,
formula(~condition))
cervical.testS4 = DESeq(cervical.testS4, fitType = "local")
```
Classify using Support Vector Machines.
```{r}
svm = classify(data = cervical.trainS4, method = "svm", normalize = "deseq",
deseqTransform = "vst", cv = 5, rpt = 3, ref = "1")
svm
```
It returns an object of class 'MLseq' and we observe that it successfully
fitted a model with 97.8% accuracy. We can access the slots of this S4 object by
```{r}
getSlots("MLSeq")
```
And also, ask about the model trained.
```{r}
trained(svm)
```
We can predict the class labels of our test data using "predict"
```{r}
pred.svm = predictClassify(svm, cervical.testS4)
table(pred.svm, relevel(cervical.testS4$condition, 2))
```
The other classification methods available are 'randomforest', 'cart' and
'bagsvm'.
### Exercise:
Train the same training data and test data using randomForest.
**Solutions:**
```{r}
rf = classify(data = cervical.trainS4, method = "randomforest",
normalize = "deseq", deseqTransform = "vst", cv = 5, rpt = 3, ref = "1")
trained(rf)
pred.rf = predictClassify(rf, cervical.testS4)
table(pred.rf, relevel(cervical.testS4$condition, 2))
```
## SessionInfo
```{r}
sessionInfo()
```
## References
1. Zararsiz G, Goksuluk D, Korkmaz S, Eldem V, Duru IP, Unver T and Ozturk A (2014). MLSeq: Machine learning interface for RNA-Seq data. R package version 1.3.0.
2. Himes, E. B, Jiang, X., Wagner, P., Hu, R., Wang, Q., Klanderman, B., Whitaker, M. R, Duan, Q., Lasky-Su, J., Nikolos, C., Jester, W., Johnson, M., Panettieri, A. R, Tantisira, G. K, Weiss, T. S, Lu and Q. (2014). “RNA-Seq Transcriptome Profiling Identifies CRISPLD2 as a Glucocorticoid Responsive Gene that Modulates Cytokine Function in Airway Smooth Muscle Cells.†PLoS ONE, 9(6), pp. e99625. http://www.ncbi.nlm.nih.gov/pubmed/24926665.
3. An Introduction to Statistical Learning with Applications in R, Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani
4. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Trevor Hastie, Robert Tibshirani, Jerome Friedman