Assignment 3
Problem Statement
High dimensional data is typically time and resource-intensive to analyze, interpret and visualize. Gene expression data (including microarray data) often consists of measurements for tens of thousands of known genes for every sample. Learning methods to inspect, reduce and display high dimensional data is necessary for many machine learning and bioinformatics problems.
Learning Objectives
- Understand the basic principles of Principal Component Analysis (PCA) and hierarchical clustering / heatmaps
- Perform PCA in R and evaluate the basic outputs of PCA
- Generate basic clustered heatmaps
Background on Microarrays
A microarray consists of thousands of specifically designed single stranded DNA sequences affixed or bound to a solid surface (glass, nylon, etc.). Extracted RNA or DNA from samples of interest are labeled with fluorescent dyes, and hybridized with the bound probes by flowing across the surface. Non-hybridized molecules are washed away and b a laser will excite the attached dye to produce light to be detected by a scanner and converted to a digital image. Finally, image processing will transform each spot (bound probe) to a numerical measurement that can be used to infer expression levels for genes.
Background on Principal Component Analysis
Principal component analysis is an exploratory data analysis technique commonly used to reduce the dimensionality of data while trying to simultaneously minimize information loss. To understand the inner workings of PCA will require an in-depth review of linear algebra and is beyond the scope of this specific assignment. We will provide several links and references if you wish to do so on your own. At its core, PCA creates new uncorrelated linear combinations of the original variables that successively maximize variance. Thus, the first principle component (PC) represents the direction of the data that explains a maximal amount of variance. The second PC represents the direction that captures the second most variance and so on and so forth. As you can infer, if a small number of PCs capture a majority of variance contained within a dataset, they can be used as a lower dimensional representation of the data.
Marisa et al. Gene Expression Classification of Colon Cancer into Molecular Subtypes: Characterization, Validation, and Prognostic Value. PLoS Medicine, May 2013. PMID: 23700391
The example intensity data was taken from the listed publication. In it, the
authors proposed the use of gene expression profiling (through microarray
technology) to generate a robust and reproducible classifier to identify
subtypes of colorectal cancer samples. They identified six subtypes that they
demonstrated to be significantly associated with distinct molecular pathways and
clinical pathologies. We have provided you a sample expression matrix that was
normalized and subjected to batch correction. We will use the
example_intensity_data.csv as our data matrix going forward.
Scaling data using R scale()
Oftentimes, it is necessary to analyze or plot data that contains variables that differ by multiple orders of magnitude. In bioinformatics, this is an extremely common situation as different genes may be expressed at wildly different levels. For instance, in typical HTS experiments and as you will see later on in the course, you may detect genes with counts ranging from the single digits to those with counts in the tens of thousands. Heatmaps of gene expression data are commonly transformed by standard scaling in order to improve the visualization while retaining the underlying pattern of expression between samples.
One such way to standardize data is the scale() function in R. Scale
essentially takes a vector of values and determines the mean and standard
deviation of the entire vector. Scale() will then subtract the mean from each
element in the vector and divide each by the standard deviation. Z-scores thus
correspond to the number of standard deviations by which a data point is above
or below the mean of all the observations (i.e. a Z-score of 0 represents a
value equal to the mean and a Z-score of 1 represents a value one standard
deviation greater than the mean value). This has the effect of representing all
of your data points on the same “scale” while preserving the pattern or profile
of differences inherent between values across different observations.
The transpose function t()
is a function covered in the textbook for this course. As covered earlier, it
converts a \(m \times n\) matrix to \(n \times m\). The built-in R function
scale() operates on a column-wise basis. Transpose the matrix such that
scaling occurs within genes rather than samples. You may also use the tidyverse
pivot functions to perform this transpose operation if you
prefer.
Using what you’ve learned in the course and prior assignments, read in the
example_intensity_data.csv to a proper format. Return the scaled matrix.
#1: We want to both center and scale our matrix.
Deliverables 1. The example_intensity_data.csv read in as a dataframe.
Proportion of variance explained
As mentioned earlier, principal components represent a successively maximal amount of variance in your dataset. It is often helpful to visualize the percentage of total variance explained by each principal component. Write a function to determine the proportion of variance explained by each principal component. Use the results generated by this function to create a bar chart displaying the variance explained by successive PCs. On this same plot, make a scatter and line plot showing the cumulative proportion of variance explained by each principal component.
#1: The summary() function in R provides top-level information about model fitting and statistical objects. #2: You may access the standard deviation of the PCA results object by $sdev
Deliverables 1. A vector with the variance explained for each PC. 2. A tibble with the labels for each PC, the variance explained and the cumulative variance explained. 3. A barchart of the proportion of individual variance explained by each principal component overlayed with a scatter plot of the cumulative sum of the variance explained for each successive principal component. Label all relevant axes and provide a descriptive caption for the figure.
Plotting and visualization of PCA
In typical HTS experiments, PCA is a common tool used to analyze the similarity between samples and determine if the experimental variable of interest (genotype, knockout, etc.) represents a large source of variance. One of the most common visualization is to plot the first two principal components (which may not always represent a majority of variance in your dataset) against each other. Typically, one will examine the pattern of clustering in this plot to see if samples are grouped together in any meaningful pattern according to important experimental variables. The value or score of each sample in terms of its principal components may be found in the pca results object accessed by $x.
When plotting individual PCs against one another, remember that different principal components explain different amounts of variance and so it will often be the case that the scale of “importance” across the axes will be different.
Deliverables 1. A scatterplot of PC1 vs PC2 from the pca_results$x values. Read in the metadata CSV and use it to annotate each sample with its corresponding assignment to the SixSubtypesClassification made in the original publication (samples belong to either c3 or c4). Label the points in the plot by color for their corresponding Six Subtypes Classification.
Hierarchical Clustering and Heatmaps
Now that we’ve discussed the basics of PCA, we will move on to discuss the use of hierarchical clustering and heatmaps. Heatmaps are a common plot used to quickly visualize the expression and pattern of many relevant observations across samples. For HTS data, they essentially depict the expression of genes (rows) for each sample (columns) as colors that denote the magnitude of expression to enable quick identification of patterns and changes between samples/experimental groups. Prior to constructing a heatmap, it is common to use some form of clustering algorithm to learn potential groups and patterns in the expression data.
One such method is known as hierarchical clustering (specifically focusing on agglomerative clustering), a method which attempts to cluster data into hierarchies where single observations begin as separate points and are successively merged into larger clusters until all points have been grouped. Typically, agglomerative clustering will begin with a dissimilarity matrix that defines how “close” all observations are by calculating some “distance” metric for every pair of observations. There are multiple “distance” metrics that can be chosen, and one of the most commonly used is simple euclidean distance. A linkage function will then use these distance metrics to define how clusters are formed based on these distance values. For example, complete linkage will define clusters by the maximum value of all pairwise distances between two different clusters.
We will be using the base R function heatmap() which by default
uses euclidean distance and the hclust()`` function to produce a clustered heatmap. Theheatmap()` function will also scale your data row-wise to aid in
visualization.
To begin, we have provided you with a representative output of differentially
expressed probes (differential_expression_results.csv) between the C3 and C4
subtypes. Filter this table to return probes with an adjusted p-value of <
.01. Use what you’ve learned to extract out the normalized intensity values for
this subset of DE probes. Use the heatmap() function in R to create a heatmap
of the normalized intensity values for these DE genes.
#1: It is important to remember that commonly used color combinations in biology (red-green) are quite difficult to discern for those with color blindness. To simplify this design decision, we are going to use RColorBrewer, which both generates appropriate color palettes and has a number of built-in color blind friendly palettes. You may use the command display.brewer.all(colorblindFriendly=TRUE) to view these curated palettes.
Deliverables 1. A list of the significant probes from the differential_expression_results.csv. 2. A matrix of the example_intensity_data.csv filtered to only contain probes contained within the list of significant probes. 3. A heatmap of the normalized intensity values for the differentially expressed probes (adjusted p-value < .01).
References
https://www.genome.gov/about-genomics/fact-sheets/DNA-Microarray-Technology
Govindarajan, R., Duraiyan, J., Kaliyappan, K. & Palanisamy, M. Microarray and its applications. Journal of Pharmacy & Bioallied Sciences 4, S310 (2012).
Trevino, V., Falciani, F. & Barrera-Saldaña, H. A. DNA microarrays: A powerful genomic tool for biomedical and clinical research. Molecular Medicine 13, 527–541 (2007).
Shlens, J. A Tutorial on Principal Component Analysis.(https://arxiv.org/abs/1404.1100)
Lever, J., Krzywinski, M. & Altman, N. Points of Significance: Principal component analysis. Nature Methods 14, 641–642 (2017).
Hastie, T., Hastie, T., Tibshirani, R., & Friedman, J. H. (2001). The elements of statistical learning: Data mining, inference, and prediction. New York: Springer.https://hastie.su.domains/ElemStatLearn/printings/ESLII_print12_toc.pdf
| output: |
| html_document: |
| css: style.css |