Cell-Cycle Scoring and Regression
Compiled: 2018-03-24
We demonstrate how to mitigate the effects of cell cycle heterogeneity in scRNA-seq data by calculating cell cycle phase scores based on canonical markers, and regressing these out of the data during pre-processing. We demonstrate this on a dataset of murine hematopoietic progenitors (Nestorowa et al., Blood 2016.You can download the files needed to run this vignette here.
library(Seurat)
# Read in the expression matrix The first row is a header row, the first
# column is rownames
exp.mat <- read.table(file = "~/Downloads/seurat_resources/nestorawa_forcellcycle_expressionMatrix.txt",
header = TRUE, as.is = TRUE, row.names = 1)
# Also read in a list of cell cycle markers, from Tirosh et al, 2015
cc.genes <- readLines(con = "~/Downloads/seurat_resources/regev_lab_cell_cycle_genes.txt")
# We can segregate this list into markers of G2/M phase and markers of S
# phase
s.genes <- cc.genes[1:43]
g2m.genes <- cc.genes[44:97]
# Create our Seurat object and complete the initalization steps
marrow <- CreateSeuratObject(raw.data = exp.mat)
marrow <- NormalizeData(object = marrow)
marrow <- FindVariableGenes(object = marrow, do.plot = FALSE, display.progress = FALSE)
marrow <- ScaleData(object = marrow, display.progress = FALSE)If we run a PCA on our object, using the variable genes we found in FindVariableGenes above, we see that while most of the variance can be explained by lineage, PC4 is split on cell-cycle genes including TOP2A and MKI67. We will attempt to regress this signal from the data, so that cell-cycle heterogneity does not contribute to PCA or downstream analysis.
marrow <- RunPCA(object = marrow, pc.genes = marrow@var.genes, pcs.print = 1:4,
genes.print = 10)## [1] "PC1"
## [1] "CORO1A" "FXYD5" "PLAC8" "H2AFY" "LAPTM5" "PKM"
## [7] "CD34" "TMEM176B" "LCP1" "LIMD2"
## [1] ""
## [1] "BLVRB" "KLF1" "ERMAP" "FAM132A" "CES2G" "RHD" "AQP1"
## [8] "SPHK1" "CAR2" "CAR1"
## [1] ""
## [1] ""
## [1] "PC2"
## [1] "APOE" "GATA2" "RAB37" "TXNIP" "VAMP5" "MFSD2B" "S100A10"
## [8] "CTLA2A" "F2R" "H2-Q6"
## [1] ""
## [1] "CTSG" "CLEC12A" "ELANE" "ANXA3" "MPO" "HP" "LY6C2"
## [8] "PRTN3" "IGSF6" "TYROBP"
## [1] ""
## [1] ""
## [1] "PC3"
## [1] "CDK6" "APOE" "FAM46A" "MS4A3" "LGALS1" "FCGR3" "RAB44"
## [8] "ELANE" "NKG7" "HDC"
## [1] ""
## [1] "FLT3" "MYL10" "LSP1" "GIMAP6"
## [5] "9030619P08RIK" "WFDC17" "DNTT" "FTH1"
## [9] "CD24A" "GCNT2"
## [1] ""
## [1] ""
## [1] "PC4"
## [1] "TRIB2" "ATP1B2" "YPEL3" "CAR1" "TCP11L2"
## [6] "LPIN1" "AQP9" "IFNGR1" "PLD4" "HIST1H2BC"
## [1] ""
## [1] "NUSAP1" "TOP2A" "FAM64A" "KIF11" "MIS18BP1" "UBE2C"
## [7] "MKI67" "CDK1" "PRC1" "SPC25"
## [1] ""
## [1] ""PCHeatmap(object = marrow, pc.use = 4, do.balanced = TRUE, label.columns = FALSE,
remove.key = TRUE)
Assign Cell-Cycle Scores
First, we assign each cell a score, based on its expression of G2/M and S phase markers. These marker sets should be anticorrelated in their expression levels, and cells expressing neither are likely not cycling and in G1 phase.
We assign scores in the CellCycleScoring function, which stores S and G2/M scores in object@meta.data, along with the predicted classification of each cell in either G2M, S or G1 phase. CellCycleScoring can also set the identity of the Seurat object to the cell-cycle phase by passing set.ident = TRUE (the original identities are stored as old.ident). Please note that Seurat does not use the discrete classifications (G2M/G1/S) in downstream cell cycle regression. Instead, it uses the quantitative scores for G2M and S phase. However, we provide our predicted classifications in case they are of interest.
marrow <- CellCycleScoring(object = marrow, s.genes = s.genes, g2m.genes = g2m.genes,
set.ident = TRUE)
# view cell cycle scores and phase assignments
head(x = marrow@meta.data)## nGene nUMI orig.ident S.Score G2M.Score Phase old.ident
## Prog_013 10211 2563089 Prog -0.12278698 -0.4370075 G1 Prog
## Prog_019 9991 3030620 Prog -0.16508027 0.5857611 G2M Prog
## Prog_031 10192 1293487 Prog -0.35666253 -0.3447439 G1 Prog
## Prog_037 9599 1357987 Prog -0.43064955 0.6771199 G2M Prog
## Prog_008 10540 4079891 Prog 0.56611382 0.1376925 S Prog
## Prog_014 10788 2569783 Prog 0.07574505 0.3079849 G2M Prog# Visualize the distribution of cell cycle markers across
RidgePlot(object = marrow, features.plot = c("PCNA", "TOP2A", "MCM6", "MKI67"),
nCol = 2)
# Running a PCA on cell cycle genes reveals, unsurprisingly, that cells
# separate entirely by phase
marrow <- RunPCA(object = marrow, pc.genes = c(s.genes, g2m.genes), do.print = FALSE)
PCAPlot(object = marrow)
We score single cells based on the scoring strategy described in Tirosh et al. 2016. See ?AddModuleScore in Seurat for more information, this function can be used to calculate supervised module scores for any gene list.
Regress out cell cycle scores during data scaling
We now attempt to subtract (‘regress out’) this source of heterogeneity from the data. For users of Seurat v1.4, this was implemented in RegressOut. However, as the results of this procedure are stored in object@scale.data (therefore overwriting the output of ScaleData), we now merge this functionality into the ScaleData function itself.
For each gene, Seurat models the relationship between gene expression and the S and G2M cell cycle scores. The scaled residuals of this model represent a ‘corrected’ expression matrix, that can be used downstream for dimensional reduction.
marrow <- ScaleData(object = marrow, vars.to.regress = c("S.Score", "G2M.Score"),
display.progress = FALSE)
# Now, a PCA on the variable genes no longer returns components associated
# with cell cycle
marrow <- RunPCA(object = marrow, pc.genes = marrow@var.genes, genes.print = 10)## [1] "PC1"
## [1] "FXYD5" "H2AFY" "CORO1A" "PKM" "CD34" "EMB" "LCP1"
## [8] "PRTN3" "PLAC8" "MPO"
## [1] ""
## [1] "BLVRB" "KLF1" "CAR2" "AQP1" "CAR1" "CES2G" "ERMAP"
## [8] "ATP1B2" "FAM132A" "TSPO2"
## [1] ""
## [1] ""
## [1] "PC2"
## [1] "CD24A" "LY6C2" "ELANE" "CLEC12A" "ANXA3" "HP" "CTSG"
## [8] "IGSF6" "ALAS1" "RGCC"
## [1] ""
## [1] "APOE" "GATA2" "CTLA2A" "RAB37" "ANGPT1" "F2R" "MUC13"
## [8] "S100A10" "ITGA2B" "TXNIP"
## [1] ""
## [1] ""
## [1] "PC3"
## [1] "GIMAP6" "LSP1" "MYL10" "FLT3"
## [5] "GCNT2" "9030619P08RIK" "WFDC17" "HLF"
## [9] "LY6A" "FTH1"
## [1] ""
## [1] "CDK6" "MS4A3" "ELANE" "FAM46A" "LGALS1" "FCGR3" "RAB44"
## [8] "APOE" "HDC" "IGSF6"
## [1] ""
## [1] ""
## [1] "PC4"
## [1] "SERPINA3G" "CD9" "MPL" "PF4" "CD63"
## [6] "ITGA2B" "CD63-PS" "PPIC" "GP1BB" "SDPR"
## [1] ""
## [1] "DNTT" "S100A10" "IRF8" "FLT3" "MPEG1" "TRIB2" "IGHM"
## [8] "CTSS" "SATB1" "HDAC7"
## [1] ""
## [1] ""
## [1] "PC5"
## [1] "CPA3" "LMO4" "HDC" "CD63-PS" "CD63" "MS4A2" "IKZF2"
## [8] "FUT8" "LTB4R1" "LPCAT2"
## [1] ""
## [1] "SDPR" "GP1BB" "PF4" "F2RL2" "TREML1" "F2R" "SLC14A1"
## [8] "MEF2C" "RAB27B" "LY86"
## [1] ""
## [1] ""# When running a PCA on only cell cycle genes, cells no longer separate by
# cell-cycle phase
marrow <- RunPCA(object = marrow, pc.genes = c(s.genes, g2m.genes), do.print = FALSE)
PCAPlot(object = marrow)
As the best cell cycle markers are extremely well conserved across tissues and species, we have found this procedure to work robustly and reliably on diverse datasets.
Alternate Workflow
The procedure above removes all signal associated with cell cycle. In some cases, we’ve found that this can negatively impact downstream analysis, particularly in differentiating processes (like murine hematopoiesis), where stem cells are quiescent and differentiated cells are proliferating (or vice versa). In this case, regressing out all cell cycle effects can blur the distinction between stem and progenitor cells as well.
As an alternative, we suggest regressing out the difference between the G2M and S phase scores. This means that signals separating non-cycling cells and cycling cells will be maintained, but differences in cell cycle phase amongst proliferating cells (which are often uninteresting), will be regressed out of the data
marrow@meta.data$CC.Difference <- marrow@meta.data$S.Score - marrow@meta.data$G2M.Score
marrow <- ScaleData(object = marrow, vars.to.regress = "CC.Difference", display.progress = FALSE)
# cell cycle effects strongly mitigated in PCA
marrow <- RunPCA(object = marrow, pc.genes = marrow@var.genes, genes.print = 10)## [1] "PC1"
## [1] "CORO1A" "FXYD5" "PLAC8" "LAPTM5" "H2AFY" "PKM"
## [7] "CD34" "TMEM176B" "LCP1" "LIMD2"
## [1] ""
## [1] "BLVRB" "KLF1" "ERMAP" "FAM132A" "CES2G" "RHD" "SPHK1"
## [8] "AQP1" "CAR2" "SLC38A5"
## [1] ""
## [1] ""
## [1] "PC2"
## [1] "CTSG" "CLEC12A" "ANXA3" "MPO" "ELANE" "PRTN3" "HP"
## [8] "LY6C2" "TYROBP" "IGSF6"
## [1] ""
## [1] "APOE" "GATA2" "TXNIP" "VAMP5" "RAB37" "MFSD2B" "S100A10"
## [8] "CTLA2A" "F2R" "H2-Q6"
## [1] ""
## [1] ""
## [1] "PC3"
## [1] "FLT3" "MYL10" "LSP1" "GIMAP6"
## [5] "WFDC17" "9030619P08RIK" "FTH1" "DNTT"
## [9] "CD24A" "GCNT2"
## [1] ""
## [1] "CDK6" "APOE" "FAM46A" "MS4A3" "RAB44" "FCGR3" "LGALS1"
## [8] "ELANE" "NKG7" "HDC"
## [1] ""
## [1] ""
## [1] "PC4"
## [1] "PF4" "ITGA2B" "GP1BB" "TUBA8" "SDPR" "NUSAP1" "PBX1"
## [8] "TOP2A" "MPL" "RAB27B"
## [1] ""
## [1] "TRIB2" "YPEL3" "IFNGR1" "ATP1B2"
## [5] "LPIN1" "CAR1" "TSPO2" "SCIN"
## [9] "HDAC7" "D930028M14RIK"
## [1] ""
## [1] ""
## [1] "PC5"
## [1] "RGCC" "CD63-PS" "CD63" "HDC" "FCGR3" "MS4A3" "PGLYRP1"
## [8] "GSTM1" "ELANE" "CPA3"
## [1] ""
## [1] "IRF8" "TUBB5" "LY86" "MPEG1" "DNTT" "TOP2A" "CTSS" "IGHM"
## [9] "TPM4" "IL12A"
## [1] ""
## [1] ""# when running a PCA on cell cycle genes, actively proliferating cells
# remain distinct from G1 cells however, within actively proliferating
# cells, G2M and S phase cells group together
marrow <- RunPCA(object = marrow, pc.genes = c(s.genes, g2m.genes), do.print = FALSE)
PCAPlot(object = marrow)