单细胞数据高级分析之消除细胞周期因素 | Removal of cell cycle effect

The normalization method described above aims to reduce the effect of technical factors in scRNA-seq data (primarily, depth) from downstream analyses. However, heterogeneity in cell cycle stage, particularly among mitotic cells transitioning between S and G2/M phases, also can drive substantial transcriptomic variation that can mask biological signal. To mitigate this effect, we use a two-step approach:

1) quantify cell cycle stage for each cell using supervised analyses with known stage-specific markers,

2) regress the effect of cell cycle stage using the same negative binomial regression as outlined above.

For the first step we use a previously published list of cell cycle dependent genes (43S phase genes, 54 G2/M phase genes) for an enrichment analysis similar to that proposed in ref. 11.

For each cell, we compare the sum of phase-specific gene expression (log10 transformed UMIs) to the distribution of 100 random background genes sets, where the number of background genes is identical to the phase gene set, and the background genes are drawn from the same expression bins. Expression bins are defined by 50 non-overlapping windows of the same range based on log10(mean UMI). The phase-specific enrichment score is the expression z-score relative to the mean and standard deviation of the background gene sets. Our final ‘cell cycle score’ (Extended Data Fig. 1) is the difference between S-phase score and G2/M-phase score.

 

For a final normalized dataset with cell cycle effect removed, we perform negative binomial regression with technical factors and cell cycle score as predictors. Although the cell cycle activity was regressed out of the data for downstream analysis, we stored the computed cell cycle score before regression, enabling us to remember the mitotic phase of each individual cell. Notably, our regression strategy is tailored to mitigate the effect of transcriptional heterogeneity within mitotic cells in different phases, and should not affect global differences between mitotic and non-mitotic cells that may be biologically relevant.

 

get.cc.score <- function(cm, N=100, seed=42) {
  set.seed(seed)
  cat('get.cc.score, ')
  cat('number of random background gene sets set to', N, '\n')
  
  min.cells <- 5
  
  cells.mols <- apply(cm, 2, sum)
  gene.cells <- apply(cm>0, 1, sum)
  cm <- cm[gene.cells >= min.cells, ]
  
  gene.mean <- apply(cm, 1, mean)
  
  breaks <- unique(quantile(log10(gene.mean), probs = seq(0,1, length.out = 50)))
  gene.bin <- cut(log10(gene.mean), breaks = breaks, labels = FALSE)
  names(gene.bin) <- rownames(cm)
  gene.bin[is.na(gene.bin)] <- 0
  
  regev.s.genes <- read.table(file='./annotation/s_genes.txt', header=FALSE, stringsAsFactors=FALSE)$V1
  regev.g2m.genes <- read.table(file='./annotation/g2m_genes.txt', header=FALSE, stringsAsFactors=FALSE)$V1
  
  goi.lst <- list('S'=rownames(cm)[!is.na(match(toupper(rownames(cm)), regev.s.genes))],
                  'G2M'=rownames(cm)[!is.na(match(toupper(rownames(cm)), regev.g2m.genes))])
  
  n <- min(40, min(sapply(goi.lst, length)))
  goi.lst <- lapply(goi.lst, function(x) x[order(gene.mean[x], decreasing = TRUE)[1:n]])
  
  bg.lst <- list('S'=get.bg.lists(goi.lst[['S']], N, gene.bin),
                 'G2M'=get.bg.lists(goi.lst[['G2M']], N, gene.bin))
  
  all.genes <- sort(unique(c(unlist(goi.lst, use.names=FALSE), unlist(bg.lst, use.names=FALSE))))
  
  expr <- log10(cm[all.genes, ]+1)
  
  s.score <- enr.score(expr, goi.lst[['S']], bg.lst[['S']])
  g2m.score <- enr.score(expr, goi.lst[['G2M']], bg.lst[['G2M']])
  
  phase <- as.numeric(g2m.score > 2 & s.score <= 2)
  phase[g2m.score <= 2 & s.score > 2] <- -1
  
  return(data.frame(score=s.score-g2m.score, s.score, g2m.score, phase))
}