Outline

Statistics in the big data era

Dimensionality reduction

Applications to single-cell rna sequencing

Background on classical statistics (1930’s)

Data today

Tables are wide - e.g. scRNAseq with few observations (cells) and many variables (genes)

Univariate analyses of multivariate data

Data from (Hastie and Tibshirani 2009)

Pairwise scatter plots for exploratory data analysis

• tedious, impractical
• risk of spurious correlations
• risk of missing information

Perform multiple linear regression using a subset of variables

• throws away information
• assumes no “noise” in \(x\)
• regression on one \(y\) at a time

Terminology

Data

• input (dependent variables): e.g. height and weight of a person, denoted by \(\mathbf{x}\); also called features, covariates, or attributes
• output (response variable): typically represented by \(y\)

Data types

Output \(y\) is typically one of

• categorical/nominal: variable takes values in finite set \(\{1; @\dots; @C\}\) of \(C\) classes
• real-valued/continuous: continuous variables
• ordinal: categorical variable that is ordered, e.g. A-F grades

Input \(x\) can be any of the above, or more complex data structures, such as images

(Murphy 2012)

Data analysis objectives

• exploratory data analysis: visualize data for quality control, discovering trends and identifying outliers among other things
• classification: given an input \(\mathbf{x}\), assign it to one of \(C\) classes in the categorical output \(y\)
• regression: fit input \(\mathbf{x}\) to continuous response \(y\); model building and prediction
• clustering: group input \(\mathbf{x}\) according to “similarity”
• discover latent factors: dimensionality reduction of \(\mathbf{x}\) to discover (unobserved) latent factors that describe most of the variability in the data; typically when features \(p >> n\) inputs

Supervised learning

Example - classification of flowers

Iris data set

Goal: learn mapping from input \(\mathbf{x}\) to output \(y \in \{1; @\dots; @C\}\) where \(C\) is number of classes

Here, scatter plots are based on features extracted from iris data set. Setosas easy to distinguish from other species.

Exploratory data analysis useful before applying any modelling/machine learning.

Supervised learning

Example - linear regression (data from (Murphy 2012))

lm(y ~ x) 

lm(y ~ poly(x, 2)) 

Linear regression minimizes the squared distance from the points to the regression line.

Unsupervised learning

Example - clustering (data from (Murphy 2012))

Height and weight

Kmeans with K=2 clusters

Goals

1. estimate distribution over number of clusters \(K\)
2. determine which cluster a point belongs to

Harder problem - we don’t know the answer beforehand. How to choose \(K\) wisely? Model selection.

Unsupervised learning

Example - dimensionality reduction with pca (data from (Murphy 2012)))

Question: how do we interpret the so-called principal components?

PCA

Developed by Pearson (1901), finally formalized by Hotelling (1931)

What is PCA?

The central idea of principal component analysis (PCA) is to reduce the dimensionality of a data set consisting of a large number of interrelated variables, while retaining as much as possible of the variation present in the data set

(Principal Component Analysis I.T. Jolliffe Springer 2002)

Geometric interpretation

What is an orthogonal projection?

Orthogonal projections and maximal variance

An animation says more than a million words…

PCA and singular value decomposition

Denote the data by the \(n \times p\) matrix \(\mathbf{X}\). There exists a factorization called singular value decomposition of the form

where

• \(\mathbf{U}\): \(n \times p\) orthogonal matrix whose columns \(\mathbf{u_j}\) are called left singular vectors
• \(\mathbf{V}\): \(p \times p\) orthogonal matrix whose columns \(\mathbf{v_j}\) are called right singular vectors
• \(\mathbf{D}\): \(p\times p\) diagonal matrix with diagonal elements \(d_1 \geq d_2 \geq \cdots \geq d_p \geq 0\).

The columns of \(\mathbf{UD}\) are score vectors of the principal components, columns of \(\mathbf{V}\) are loadings. The first principal component can also be written \(\mathbf{z_1}=\mathbf{X}v_1=\mathbf{U}\mathbf{D}\mathbf{V^T}v_1=\mathbf{u_1}d_1\) where \(\mathsf{Var(}{\mathbf{z_1}}\mathsf{)}=\frac{d^2_1}{n}\) and so on.

Projection points

Back to our PCA plot.

For a given point \(x_i\), the projection point on PC1 is given by \(u_{i1}d_1v_1\). Here, \(x_i\) has a higher PC1 score than \(x_j\).

Interpretation of PCs - loadings

Interpretation of PCs - overlaying metadata

Height and weight

Coloring by metadata also facilitates interpretation.

Standardization

Since PCA aligns PCs along directions of maximum variance, it is sensitive to differences in scale, e.g. caused by different measurement units. Therefore it is standard practice to scale variables to unit variance when running PCA.

PCA as data compression algorithm and the PCA model

\(\mathbf{X}\) can be approximated by using only the first \(q\) principal components, i.e.:

\[ X \approx U D_q V_q^T \]

In essence, PCA can be seen as a compression algorithm in which we throw away the smallest terms (smallest variances) in the diagonal matrix and approximate \(\mathbf{X}\). Any given PC \(j\) explains a proportion \(\pi_j\) of the variance which is defined as \(\pi_j = d_j/\sum_{i=1}^pd_i\).

By introducing an error term, the PCA model can be written as

\[ \mathbf{X} = U D_q V_q^T + E \]

What PCA is revisited

Unfortunately […] there is a lot of confusion about what PCA is really about

(Pachter 2014)

PCA is

• not a statistical method according to one textbook
• a statistical method according to another
• more useful as a visualization technique then as an analytical method
• Principal components is a mathematical algorithm

Other dimensionality reduction techniques

• principal components analysis
• multidimensional scaling
• partial least squares to latent structure analysis (PLS)
• orthogonal partial least squares to latent structures analysis (OPLS)
• PLS-DA (discriminant analysis)
• OPLS-DA
• biplot analysis
• k-means clustering
• hierarchical clustering
• self-organizing maps
• t-SNE

t-distributed stochastic neigbor embedding (t-SNE)

In contrast to PCA, t-SNE (Maaten and Hinton 2008) is a non-linear dimensionality reduction technique

It models high-dimensional points in such a way that points close in high-dimensional space end up close low-dimensional space

The algorithm is non-linear and adapts to the underlying data density, using different transformations for different regions

The perplexity input parameter approximates the number of neigbors a point has. Typical value are between 5 and 50.

In addition, the algorithm is stochastic, meaning successive runs with identical parameters may produce different results

The perplex perplexity effect

Choosing the number of iterations

Cluster size means nothing

Inbetween cluster distance might not mean anything

Example data set

The following examples are based on a study by (Zeisel et al. 2015). The data can be obtained from (Amit Zeisel and Linnarsson 2015). The analyses have been performed in R using utility functions from the scater package (“Scater,” n.d.). The analysis is based on the simpleSingleCell bioconductor tutorial (Aaron T. L. Lun and Marioni 2017).

The authors sequenced RNA from \(>\) 3000 single cells to “perform a molecular census of the primary somatosensory cortex (S1) and the hippocampal CA1 region”.

Quality control

PCA can be used for quality control by basing it on the quality metrics (coldata) for each cell, e.g.:

• the total number of reads
• the total number of features (e.g. expressed genes)
• proportion of mitochondrial or spike-in reads
• batch / lane

or by using (filtered) expression values to get an overview of clustering patterns.

Outliers may be indicative of low-quality cells. Use quality metrics to aid interpretation of principal components.

Quality control

# data has been loaded and preprocessed and stored in a SingleCellExperiment object
sce <- calculateQCMetrics(sce, feature_controls=list(Spike=is.spike, Mt=is.mito))  

plotPCA(sce, pca_data_input="coldata") 

Interpretation? No clear outlier or cluster. Add metadata information.

Quality control

Add metadata information (tissue and total.mRNA.mol) to facilitate interpretation.

plotPCA(sce, colour_by="tissue", size_by="total.mRNA.mol",
        pca_data_input="coldata", legend="none") 

First PC seems to separate by number of mRNA molecules in cell.

Quality control

Size by number of expressed genes (calculated by CalculateQCMetrics)

plotPCA(sce, colour_by="tissue", size_by="total_features",
        pca_data_input="coldata", legend="none") 

Even clearer trend - cells with few features at the right part of the plot.

Quality control

Look at higher principal components.

plotPCA(sce, colour_by="tissue", size_by="total_features",
        pca_data_input="pdata", legend="none", ncomponents=3) 

Quality control

PCA done on expression values (rowdata)

plotPCA(sce, pca_data_input="rowdata", colour_by="tissue",
        ncomponents=2, exprs_values="counts", 
        size_by=as.data.frame(rep(2, length(sce$tissue)))) 

Quality control

Size by library size and number of genes (calculated by calculateQCMetrics)

plotPCA(sce, colour_by="tissue", exprs_values="counts",
        size_by=as.data.frame(sce$total_counts/1e6),
        legend="none", ncomponents=2) 

plotPCA(sce, colour_by="tissue",
        exprs_values="counts",
        size_by=as.data.frame(sce$total_features),
        legend="none", ncomponents=2) 

Data processing and filtering

There are many more quality control steps to perform. Here, we have only discussed the use of PCA.

Before proceeding with data exploration with t-SNE/PCA and clustering we need to process and filter the data via a number of steps:

• outlier removal: remove outliers based on distribution of QC metrics, such as total number of features.
• remove uninteresting genes: remove genes with low average expression
• normalization: account for library size
• check for technical factors: block for uninteresting factors, e.g. sex
• identify correlated HVGs: highly variable genes
• remove batch effects: remove effects accounted for by blocking factor

Data exploration of correlated highly variable genes

tsne1 <- plotTSNE(sce, exprs_values="norm_exprs", colour_by=top.hvg,
    perplexity=10, rand_seed=100, feature_set=chosen) + fontsize
tsne2 <- plotTSNE(sce, exprs_values="norm_exprs", colour_by="Mog",
    perplexity=10, rand_seed=100, feature_set=chosen) + fontsize
multiplot(tsne1, tsne2, cols=2) 

pca1 <- plotPCA(sce, exprs_values="norm_exprs", colour_by=top.hvg) + fontsize
pca2 <- plotPCA(sce, exprs_values="norm_exprs", colour_by="Mog") + fontsize
multiplot(pca1, pca2, cols=2) 

Clustering cells into putative subpopulations

Perform hierarchical clustering on Euclidean distances between cells and use a dynamic tree cut to define clusters. The block-like pattern consistent with dimensionality reduction shown previosly.

Detecting marker genes between subpopulations

The goal: identify marker genes for each subpopulation. Many markers are not uniquely up- or down-regulated in a cluster; rather a combination of markers is used to characterize a subpopulation

Summary

Dimensionality reduction techniques, with a focus on PCA (and t-SNE)

Using PCA:

• find latent structures
• loadings, marker data to aid interpretation of components
• center data!
• scaling matters

t-SNE:

• perplexity
• distances and cluster sizes don’t necessarily mean anything

References