--- title: D. Machine Learning author: Martin Morgan (mtmorgan@fredhutch.org)
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