About Install Get Started Frequently Asked Questions Frequently Requested Vignettes Contact
Cell-Cycle Scoring and Regression

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)