Normalization and clustering of single-cell RNA-seq data

_{^{Owens (2012) Nature 491, 27–29.}}

_{^{www.fluidigm.com}}

The brain is made up of 100’s if not 1000’s different cell types.
We need a rational way to identify and classify them.

• 285 Layer 4 cells
• 657 Layer 5 cells
• 307 Layer 6 cells

Low sequencing depth: 192 cells per Illumina lane (average 1.2M reads per cell)

Each level of the factor of interest, say layer of origin, is observed in each batch.

We can only isolate one cell type per animal / batch.

_{^{See also in bioRxiv:
Hicks et al. (2015) http://dx.doi.org/10.1101/025528}}

Need to have multiple batches per condition and account for batch effects in the design matrix (nested design).

Note that mouse and C1 run effects are still confounded!

_{^{See also in bioRxiv:
Tung et al. (2016) http://dx.doi.org/10.1101/062919}}

To account for batch $j(i)$ in condition $i$, we can model the log-expression of each sample $k$ as

$$y_{ijk} = \mu + \alpha_i + \beta_{j(i)} + \varepsilon_{ijk},$$

subject to the $n+1$ constraints

$$\sum_{i=1}^n \alpha_i = 0; \\ \sum_{j=1}^{n_i} \beta_{j(i)} = 0.$$

Absolute correlation between PCA of log(TPM+1) and QC scores.

PC1 of the QC score matrix stratified by batch.

RUV can be used to estimate $W \alpha$ using negative control genes.

_{^{Gagnon-Bartsch & Speed (2012) [http://dx.doi.org/10.1093/biostatistics/kxr034]
Risso et al. (2014) [http://dx.doi.org/10.1038/nbt.2931]
RUVSeq package [http://bioconductor.org/packages/RUVSeq/]}}

1. Apply a (combination of) normalization method(s).

   • Global scaling, e.g., DESeq, TMM.
   • Full-quantile (FQ)
   • Unknown factors of unwanted variation, e.g., RUV.
   • Known factors of unwanted variation, e.g., regression on QC measures, C1 batch.

2. Rank the normalizations using a set of performance scores.

_{^{Michael Cole, Nir Yosef, Sandrine Dudoit}}

• BIO_SIL: Average silhouette width by biological condition.
• BATCH_SIL: Average silhouette width by batch.
• PAM_SIL: Maximum average silhouette width for PAM clusterings, for a range of user-supplied numbers of clusters.
• EXP_QC_COR: Maximum squared Spearman correlation between count PCs and QC measures.
• EXP_UV_COR: Maximum squared Spearman correlation between count PCs and factors of unwanted variation (preferably derived from other set of negative control genes than used in RUV).
• EXP_WV_COR: Maximum squared Spearman correlation between count PCs and factors of wanted variation (derived from positive control genes).
• RLE_MED: Mean squared median relative log expression (RLE).
• RLE_IQR: Mean inter-quartile range (IQR) of RLE.

Points color-coded by average score.

Heatmap of the 100 most variable genes

Heatmap of the 100 most variable genes

In the literature, most approaches can be summarized by three steps.

1. Dimensionality reduction (e.g., PCA, t-SNE, most variable genes).
2. Compute a distance matrix between samples in the reduce space.
3. Clustering based on a partitioning method (e.g., PAM, k-means).

For each step there are many tuning parameters. E.g.,

• How many principal components?
• Which distance?
• How many clusters?

Given a base cluster algorithm

• Generate a single candidate clustering using
  • resampling (to find robust clusters)
  • sequential clustering (to find stable clusters)
• Repeat the procedure for different algorithms and tuning parameters to generate a collection of candidate clusterings
• Identify a consensus over the different candidates

Implemented in the R/Bioconductor package clusterExperiment: http://bioconductor.org/packages/clusterExperiment

_{^{Elizabeth Purdom}}

Given an underlying clustering strategy, e.g., k-means or PAM with a particular choice of k, we repeat the following:

• Subsample the data, e.g. 70% of samples.
• Find clusters on the subsample.
• Create a co-clustering matrix D:
  • % of subsamples where samples were in the same cluster.

Our sequential clustering works as follows.

• Range over k in PAM clustering using the subsampling strategy.
• The cluster that remains stable across values of k is identified and removed.
• Repeat until no more stable clusters are found.

_{^{Inspired by the "tight clustering" algorithm
Tseng and Wong (2005) http://dx.doi.org/10.1111/j.0006-341X.2005.031032.x}}

• Find cluster gene expression signatures (marker genes).
• Standard solutions:
  • F-test for any difference between clusters
  • all pairwise comparisons
• Our solution:
  • Create a hierarchy of clusters
  • Select appropriate contrasts that compare sister nodes

Functions to

• Generate a collection of candidate cluster labels.
• Find a consensus clustering.
• Merge clusters.
• Find cluster signatures.
• Visualize clustering results.

Available at http://bioconductor.org/packages/clusterExperiment

_{^{Liam Purvis, Elizabeth Purdom}}

GOAL: High-resolution view of transcriptional changes during differentiation and neurogenesis.

Questions:

• Where does the neuronal lineage branch off?
• Which genes are differentially expressed throughout this process?
• Which genes are driving cell fate decisions?

_{^{Kelly Street}}

• Isolate cells from the OE system by FACS
• Sync cells temporally by conditional knockout of p63, an inhibitor of differentiation.
• Quantify gene expression with single-cell RNA sequencing.

_{^{Russell Fletcher}}

_{^{Kelly Street}}

_{^{Kelly Street}}

_{^{Kelly Street}}

Flexible, supervised branching lineage reconstruction.

Input: scRNA-seq data after normalization, clustering and dimensionality reduction.

R package available at: https://github.com/kstreet13/slingshot

_{^{Kelly Street, Elizabeth Purdom, Sandrine Dudoit}}

_{^{Kelly Street}}

_{^{Kelly Street}}

_{^{Kelly Street, Russell Fletcher, Diya Das}}

_{^{Kelly Street, Russell Fletcher, Diya Das}}