1 of 57

Introduction to �single cell data science and tidyomics

Stephanie Hicks, PhD | Associate Professor

Biostatistics, Bloomberg School of Public Health

Biomedical Engineering, Whiting School of Engineering

Faculty member, Malone Center for Engineering in Healthcare��stephaniehicks.com/talks �Twitter: @stephaniehicks

2 of 57

How are genomic data generated?

3 of 57

ACTGACCTAGATCAGTCGATCGATCGTATACGATTACAAAATCATCGGCAT

DNA

data generation

4 of 57

data generation

GATCGATCGTATACGAT

Fragments

ACTGACCTAGATCAGTC

TACAAAATCATCGGCAT

ACTGACCTAGATCAGTCGATCGATCGTATACGATTACAAAATCATCGGCAT

DNA

5 of 57

data generation

GATCGATCGTATACGAT

Fragments

ACTGACCTAGATCAGTC

TACAAAATCATCGGCAT

Reads

TACAAAATCA

AGATCAGTC

GATCGATCG

ACTGACCTAGATCAGTCGATCGATCGTATACGATTACAAAATCATCGGCAT

DNA

6 of 57

@22:16362385-16362561W:ENST00000440999:2:177:-40:244:S/2

CCAGCCCACCTGAGGCTTCTTTTTCCTTCCCAAGCCACATCACCATCCTGGTGGAACTCTCCTGTGAGGACAGCCA

+

GGFF<BB=>GBGIIIIIIIIIIIIIIEGEHGHHIIIIIIIIHFHBB2/:=??EGGGEGFHHIHHEDBD?@@DDHHD

@22:16362385-16362561W:ENST00000440999:3:177:-56:294:S/2

GCGTGAGCCACAGGGCCCAGCCCACCTGAGGCTTCTTTTTCCTTCCCAAGCCACATCACCATCCTGGTGGAACTCT

+

@=ABBBBIIIIIIIIHHGGGGIIDBDIIIIIIGIIIIHIIIIHFDD@BBDBGGFIDEE8DCC/29>BGFCGHHHGF

@22:16362385-16362561W:ENST00000440999:4:177:137:254:S/1

TCACCATCCTGGTGGAACTCTCCTGTGAGGACAGCCAAGGCCTGAACTACCTGCaGTGGGGAGCACCTCAGGGTTT

+

DDGBBCGGGIGGGBDDDHIIGGDGD77=BDIIIIIIIIFHHHHIIIHEFFHGGDD8A>DEGHHIFDDHH8@BEDDI

@22:16362385-16362561W:ENST00000440999:5:177:68:251:S/2

AGGGTTTGCCCAGGCAACCAGCCAGCCCTGGTCCAAGGCATCCTGGAGCGAGTTGTGGATGGCAAAAAGACNCGCC

+

HIGHIHFHEGE4111:.;8@?@HDIIIIIIIEGGIHHHIIGA?=:FIIIDD8.02506A8=AC#############

@22:16362385-16362561W:ENST00000440999:6:177:348:453:S/1

AAGGCCTGAACTACCTGCGGTGGGGAGCACCTCAGGGTTTGCCCAGGCAACCAGCCAGCCCTGGTCCAAGGCATCC

+

B9?@8=42:E@GDEDIIIIIGGHIIIFBEEAGIIDIIDHHGGHIIEGEIIIIIHIHFHFFEEFGGGGGB88>:DGH

@22:51205934-51222090C:ENST00000464740:132:612:223:359:S/2

GGAAGTATGATGCTGATGACAACGTGAAGATCATCTGCCTGGGAGACAGCGCAGTGGGCAAATCCAAACTCATGGA

+

IIEHHHHHIIIIIIIHGGDGHHEDDG8=;?==19;<<>D@@GGGIIHIIHGGDDHGBA=ABEG@@DFCCAA<:=>8

@22:51205934-51222090C:ENST00000464740:125:612:-1:185:S/1

TGGAGTGCGCTGCGGCGCGAGCTGGGCCGGCGGGCGTGGTTCGAGAGCGCGCAGAGTCCAGACTGGCGGCAGGGCC

+

GGFF<BB=>GBGIIIIIIIIIIIIIIEGEHGHHIIIIIIIIHFHBB2/:=??EGGGEGFHHIHHEDBD?@@DDHHD

@22:16362385-16362561W:ENST00000440999:3:177:-56:294:S/2

GCGTGAGCCACAGGGCCCAGCCCACCTGAGGCTTCTTTTTCCTTCCCAAGCCACATCACCATCCTGGTGGAACTCT

+

GGFF<BB=>GBGIIIIIIIIIIIIIIEGEHGHHIIIIIIIIHFHBB2/:=??EGGGEGFHHIHHEDBD?@@DDHHD

@22:16362385-16362561W:ENST00000440999:3:177:-56:294:S/2

GCGTGAGCCACAGGGCCCAGCCCACCTGAGGCTTCTTTTTCCTTCCCAAGCCACATCACCATCCTGGTGGAACTCT

+

billions more

7 of 57

Let’s talk about one type of genomic data:

�Gene expression (or RNA)

8 of 57

Reads��Need to figure out what part of the DNA they came from →

Reference genome (DNA)

ACTGACCTAGATCAGTCGATCGATCGTATACGA

Gene XYZ

Sum reads

23 reads from �Gene XYZ

9 of 57

23

0 4 4 6 1 2

9 2 1 0 10 3

Gene 1

Gene 2

Gene G

.

Sample 1

Sample 6

Sample 2

.

Chr Start End ��chr1 10 500��chr1 400 882��chr22 12 100

Important �bookkeeping

Data we typically analyze

10 15 3 0 1

10 of 57

How can we analyze this genomic data?

11 of 57

Open-source, open development software project�
Began in 2001�
Big priorities: reproducible research and high-quality documentation

Vignettes
Diverse community support
Workflows (super helpful for new users)
Teaching resources

(and open development)

http://www.bioconductor.org

You can think of Bioconductor as a cousin to CRAN, but specific for the analysis of biological data.

Bioconductor is an open source and open development software project for computational biology.

Open source means anyone can read and modify the underlying code.
Open development means that anyone can contribute to the project.

Began in 2001 by Robert Gentleman and others (now includes 12 full-time core developers to build and maintain core Bioc infrastructure)

Big priorities (focused around reproducible research addressing domain-specific challenges with high-quality documentation)

Each software package must contain a vignette
Big community support forum with Qs on common problems answered by experts in the field
Slack workspace (open to anyone)
Many workflows exist to help users get started analyzing biological data
Lots of great teaching resources
Strongly encourages use of existing data and software infrastructure from Bioc packages to enable interoperability

12 of 57

Standard Bioconductor data structure

Lee et al. (2019) Genome Biology

GenomicRanges (GRanges)

13 of 57

Create GRanges object

14 of 57

Things you can do with GRanges objects

15 of 57

Tidy data" is organized for programming

16 of 57

The pipe

command | command | command > output.txt

"Pipes rank alongside the hierarchical file system and regular expressions as one of the most powerful yet elegant features of Unix-like operating systems." (http://www.linfo.org/pipe)��(In R we use '|>' instead of '|' to chain operations)

17 of 57

Verb-based operations

In the R package dplyr:

mutate() adds new variables that are functions of existing variables.
select() picks variables based on their names.
filter() picks cases based on their values.
slice() picks cases based on their position.
summarize() reduces multiple values down to a single summary.
arrange() changes the ordering of the rows.
group_by() perform any operation by group.

(https://dplyr.tidyverse.org)

18 of 57

Summarize after grouping

A useful paradigm is to group data and then summarize:

19 of 57

Summarize after grouping

A useful paradigm is to group data and then summarize:

20 of 57

And genomic range data are often already tidy

21 of 57

Tidy advantages

Many already familiar with tidyverse packages (e.g. dplyr and ggplot2)
Avoid intermediate variables, e.g.:

```{r eval=FALSE}� dat3 <- dat2[dat2$signal > x]� ```

Aim is for human readable code

22 of 57

Why "tidy analysis" for genomics?

Encourages exploration
Encourages efficiency: fewer calls out of R
Generalizes from simple to complex cases
Developer side: modularity is easier to maintain

23 of 57

Genomic verbs / action + �tidy data = plyranges

Michael Lawrence

Di Cook

Stuart Lee

Goal: Write human readable analysis workflows�
Idea: Define an API (i.e. extend dplyr) that maps relational genomic algebra to “verbs” that act on “tidy” genomic data

Another great idea: Borrow dplyr’s syntax and design principles�
And another great idea: Compose verbs together with pipe operator from magrittr

24 of 57

Recall non-tidy code with GRanges object

25 of 57

26 of 57

27 of 57

Piping directly into plots!

Lee et al. (2019) Genome Biology

28 of 57

(source)

29 of 57

Single-cell transcriptomics (gene expression)

30 of 57

A “typical” single-cell analysis workflow

Amezquita et al. 2020, Nature Methods

31 of 57

Amezquita et al. 2020, Nature Methods

32 of 57

Amezquita et al. 2020, Nature Methods

33 of 57

Standardized ways to store single-cell data

Amezquita et al. 2020, Nature Methods

34 of 57

Demo

https://www.stephaniehicks.com/singlecelldemo

35 of 57

Extra slides

36 of 57

Quality control

37 of 57

Motivation for Quality Control

Low-quality scRNA-seq libraries can come from many places:

Cell damage during dissociation
Failure in library preparation (e.g., inefficient reverse transcription or PCR amplification)

These manifest as “cells” with

low total counts
few expressed genes
high mitochondrial or spike-in proportions, etc

Often form their own distinct cluster(s) (complicating interpretation of results)

38 of 57

Example QC metric: library size

Library size (= total number of reads/UMIs across relevant features per cell)

Typically, we will consider the relevant features to be the endogenous genes
Cells with small library sizes are of low quality as the RNA has been lost at some point during library preparation, either due to cell lysis or inefficient cDNA capture and amplification

39 of 57

Common choices for QC metrics

Can calculate these QC metrics from the scater R/Bioc package using

perCellQCMetrics() function

40 of 57

Common choices for QC metrics

Can calculate these QC metrics from the scater R/Bioc package using

addPerCellQCMetrics() function

41 of 57

Normalization

42 of 57

Normalization and log-transformation

Now that we've removed problematic cells, need to put counts on a comparable scale before/during downstream analysis

Main factors to account for:

library size / sequencing depth
batch effects (known or unknown)

Most common approaches:

Scaling normalization (either cell- OR cell- and gene-specific) followed by
Log-transformation (and variance stabilization broadly see this post)

These systematic differences typically arise from technical differences in cDNA capture or PCR amplification efficiency across cells, attributable to the difficulty of achieving consistent library preparation with minimal starting material
Normalization aims to remove these differences such that they do not interfere with comparisons of the expression profiles between cells.
This ensures that any observed heterogeneity or differential expression within the cell population are driven by biology and not technical biases.

Scaling normalization, which is the simplest and most commonly used class of normalization strategies.

This involves dividing all counts for each cell by a cell-specific scaling factor, often called a “size factor” (Anders and Huber 2010).

The assumption here is that any cell-specific bias (e.g., in capture or amplification efficiency) affects all genes equally via scaling of the expected mean count for that cell (definitely not always true)
The size factor for each cell represents the estimate of the relative bias in that cell, so division of its counts by its size factor should remove that bias.
The resulting “normalized expression values” can then be used for downstream analyses such as clustering and dimensionality reduction.

The log-transformation is useful as differences in the log-values represent log-fold changes in expression.
This is important in downstream procedures based on Euclidean distances, which includes many forms of clustering and dimensionality reduction.
By operating on log-transformed data, we ensure that these procedures are measuring distances between cells based on log-fold changes in expression. Or in other words, which is more interesting - a gene that is expressed at an average count of 50 in cell type A and 10 in cell type B, or a gene that is expressed at an average count of 1100 in A and 1000 in B?
Log-transformation focuses on the former by promoting contributions from genes with strong relative differences.

43 of 57

Library size normalization

Library size = the total sum of counts across all genes for each cell�

“library size factor” for each cell is directly proportional to its library size where the proportionality constant is defined such that the mean size factor across all cells is equal to 1

44 of 57

Feature selection

45 of 57

Motivation for feature selection

Often a primary goal of scRNA-seq analysis is to characterize heterogeneity across cells

Dim reduction & clustering combine per-gene diffs into a similarity metric btw cells

The choice of which genes to use can heavily impact downstream analyses

Keep: genes with information about the biology of the system
Remove: genes that contain random noise

Feature selection aims to keep #1 and remove #2 to reduce the size of the data��Idea: select most variable genes based on expression across the population

46 of 57

Mean-Variance of log-normalized expression values

Simplest approach: Quantify per-gene variation by computing the variance of the log-normalized expression values (i.e. “log-counts”) for each gene across all cells

Feature selection and downstream analysis both on same log-scale

Genes with the largest variances in log-values will contribute most to the Euclidean distances between cells during procedures like clustering and dimensionality reduction.
By using log-values here, we ensure that our quantitative definition of heterogeneity is consistent throughout the entire analysis.

Idea: just rank genes by the empirical variance (most to least)

But this doesn’t capture the mean-variance relationship in seq data

The simplest approach to quantifying per-gene variation is to compute the variance of the log-normalized expression values (i.e., “log-counts” ) for each gene across all cells
The advantage of this approach is that the feature selection is based on the same log-values that are used for later downstream steps.
In particular, genes with the largest variances in log-values will contribute most to the Euclidean distances between cells during procedures like clustering and dimensionality reduction.
By using log-values here, we ensure that our quantitative definition of heterogeneity is consistent throughout the entire analysis.

…

Calculation of the per-gene variance is simple but feature selection requires modelling of the mean-variance relationship.

The log-transformation is not a variance stabilizing transformation in most cases, which means that the total variance of a gene is driven more by its abundance than its underlying biological heterogeneity.

To account for this effect, we use the modelGeneVar() function to fit a trend to the variance with respect to abundance across all genes

47 of 57

Mean-Variance of log-normalized expression values

Total var = biological var (interesting) + � technical var (uninteresting)

48 of 57

Mean-Variance of log-normalized expression values

49 of 57

Selecting highly variable genes (HVGs)

We use the getTopHVGs() function to subset of HVGs to be used in downstream analyses
Typically select top n genes with the largest biological components
Nice because a user can directly control the number of genes retained (and hence the computational complexity of downstream analyses)

50 of 57

Dimensionality reduction

51 of 57

Motivation for dimensionality reduction

A common task in scRNA-seq analysis is to compare cells based on their expression values across multiple genes.

For example, clustering aims to identify cells with similar transcriptomic profiles by computing Euclidean distances across genes.

In these applications, each individual gene represents a dimension of the data.

More intuitively, if we had a scRNA-seq data set with two genes, we could make a two-dimensional plot where each axis represents the expression of one gene and each point in the plot represents a cell.

This concept can be extended to data sets with thousands of genes where each cell’s expression profile defines its location in the high-dimensional expression space.

52 of 57

Motivation for dimensionality reduction

Dimensionality reduction aims to reduce the number of dimensions in the data.

This is possible because different genes are correlated if they are affected by the same biological process (i.e. don’t need to store separate info for individual genes)

Instead can compress multiple features into a single dimension, e.g., an “eigengene”

This reduces computational work in downstream analyses like clustering, as

calculations only need to be performed for a few dimensions rather than thousands of genes
reduces noise by averaging across multiple genes to obtain a more precise representation of the patterns in the data
enables effective plotting of the data, for those of us who are not capable of visualizing more than 3 dimensions

53 of 57

Principal Components Analysis (PCA)

Principal components analysis (PCA)

Finds axes in high-dimensional space that capture the largest amount of variation
The first axis (or “principal component”, PC) is chosen such that it maximizes this variance
The next PC is chosen such that it is orthogonal to the first and captures the greatest remaining amount of variation, and so on.

Novembre et al. 2008

54 of 57

How do we apply PCA to scRNA-seq data?

Apply PCA on the the log-normalized expression values using the fixedPCA() function from the scran R/Bioconductor package.

Default: computes first 50 PCs and stores in reducedDims() slot of a SCE object

Also can compute approximate PCA methods!

55 of 57

Things slow down with large numbers!

Remedies:

fast approximate PCA e.g. irlba - implicitly restarted Lanczos bidiagonalization algorithm
Work with sparse matrices or HDF5 files

56 of 57

Visualizing the PCs

57 of 57

Non-linear methods for visualization

t-stochastic neighbor embedding
Uniform manifold approximation and projection (UMAP