寫在前面
在scRNAseq數(shù)據(jù)中，不同細(xì)胞處在不同細(xì)胞周期耸三，如果不進(jìn)行異質(zhì)性校正的話乱陡，會(huì)對(duì)結(jié)果有很大的影響。??
常用的方法是根據(jù)經(jīng)典的細(xì)胞周期基因進(jìn)行評(píng)分仪壮，即cell cycle scores, 然后在預(yù)處理時(shí)將score納入回歸中憨颠。??
用到的包
rm(list = ls())
library(Seurat)
library(tidyverse)
library(ggsci)
示例數(shù)據(jù)
exp.mat <- read.table(file = "./nestorawa_forcellcycle_expressionMatrix.txt",
header = TRUE, as.is = TRUE, row.names = 1)
Seurat標(biāo)準(zhǔn)流程處理
4.1 創(chuàng)建Seurat對(duì)象
marrow <- CreateSeuratObject(counts = exp.mat)
marrow
4.2 歸一化（Normalization）
marrow <- NormalizeData(marrow)
4.3 尋找高變基因
marrow <- FindVariableFeatures(marrow, selection.method = "vst")
4.4 標(biāo)準(zhǔn)化
marrow <- ScaleData(marrow, features = rownames(marrow))
查看并提取細(xì)胞周期基因
Seurat包內(nèi)置了細(xì)胞周期的相關(guān)基因集，我們來看一下不同周期里的基因吧积锅。??
這里提供的是人的細(xì)胞周期的相關(guān)基因爽彤，如果你做的是小鼠，你需要做一下轉(zhuǎn)換缚陷，解決方案如下：??
https://www.r-bloggers.com/2016/10/converting-mouse-to-human-gene-names-with-biomart-package/ （如果大家需要講解的人多适篙，以后可以專門寫一期這個(gè)方面的東西，歡迎留言蹬跃。）
5.1 查看細(xì)胞周期基因集
Note! 這里演示我們用的是老版本匙瘪，這里補(bǔ)充一下新版本，cc.genes.updated.2019蝶缀。 ??
cc.genes
5.2 提取細(xì)胞周期基因
這里我們將不同周期的基因提取出來丹喻，即S期,G2期和M期。??
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes
主成分分析
這里我們可以看到一些細(xì)胞周期基因在PC8和PC10是有顯著差別的, 如TOP2A和MKI67等翁都。??
marrow <- RunPCA(marrow, features = VariableFeatures(marrow),
ndims.print = 1:10, nfeatures.print = 10)
DimHeatmap(marrow, dims = c(8, 10),
fast = T)
計(jì)算細(xì)胞周期評(píng)分
7.1 計(jì)算評(píng)分
現(xiàn)在我們可以根據(jù)這些細(xì)胞周期基因開始計(jì)算評(píng)分了碍论。??
marrow <- CellCycleScoring(marrow,
s.features = s.genes,
g2m.features = g2m.genes,
set.ident = TRUE)
head(marrow[[]])
7.2 可視化-ridgeplot
看一下幾個(gè)細(xì)胞周期基因的分布情況。??
RidgePlot(marrow,
features = c("PCNA", "TOP2A", "MCM6", "MKI67"),
cols = pal_npg("nrc", alpha = 0.7)(3),
ncol = 2)
7.3 可視化-PCA
我們用細(xì)胞周期基因做一下PCA柄慰，有明顯的分別鳍悠。??
marrow <- RunPCA(marrow, features = c(s.genes, g2m.genes))
DimPlot(marrow,
cols = pal_npg("nrc", alpha = 0.7)(3))
排除細(xì)胞周期異質(zhì)性的影響
計(jì)算好細(xì)胞周期的評(píng)分以后税娜，我們就可以在標(biāo)準(zhǔn)化的時(shí)候加入這個(gè)變量了，去除它的影響藏研。??
8.1 開始去除
marrow <- ScaleData(marrow,
vars.to.regress = c("S.Score", "G2M.Score"),
features = rownames(marrow))
8.2 可視化-PCA
這個(gè)時(shí)候我們?cè)儆?strong>細(xì)胞周期基因做一下PCA看看結(jié)果敬矩，成功聚在一起，沒有明顯的區(qū)分啦蠢挡。
marrow <- RunPCA(marrow, features = VariableFeatures(marrow), nfeatures.print = 10)
marrow <- RunPCA(marrow, features = c(s.genes, g2m.genes))
DimPlot(marrow,
cols = pal_npg("nrc", alpha = 0.7)(3))
可選步驟
以上講述的方法刪除了所有與細(xì)胞周期相關(guān)的信號(hào)值弧岳。??
在某些情況下，如分化過程中（如小鼠造血）业踏，干細(xì)胞處于靜止狀態(tài)禽炬，而分化細(xì)胞正處于增殖狀態(tài)（反之亦然）。??
在這種情況下勤家，回歸所有細(xì)胞周期效應(yīng)腹尖，會(huì)影響干細(xì)胞和祖細(xì)胞的區(qū)分。??
所以伐脖，在這里我們采用G2M期和S期得分之間的差值進(jìn)行回歸热幔。
marrow$CC.Difference <- marrow$S.Score - marrow$G2M.Score
marrow <- ScaleData(marrow, vars.to.regress = "CC.Difference", features = rownames(marrow))
可視化一下吧！~??
這里雖然細(xì)胞群聚在一起讼庇，但G1期和G2M期/S期是可以區(qū)分開的断凶。
marrow <- RunPCA(marrow, features = VariableFeatures(marrow), nfeatures.print = 10)
marrow <- RunPCA(marrow, features = c(s.genes, g2m.genes))
DimPlot(marrow,
cols = pal_npg("nrc", alpha = 0.7)(3))
<img src="https://upload-images.jianshu.io/upload_images/24475539-53b3ce362eeeccd8.png" alt="鮮肉" style="zoom:25%;" />
<center>最后祝大家早日不卷!~</center>
需要示例數(shù)據(jù)的小伙伴，在公眾號(hào)回復(fù)
Cellcycle獲取吧巫俺！點(diǎn)個(gè)在看吧各位~ ?.???? ??? ?
本文由mdnice多平臺(tái)發(fā)布
