斯坦福大學(xué)Satija lab的 Seurat v3.1 guidline于近日更新啦魄懂！其中包括許多個(gè)性化的模塊失乾，其中我個(gè)人比較感興趣的是Cell-Cycle Scoring and Regression模塊,因?yàn)樵跅l件干預(yù)的情況下扰她，部分細(xì)胞處于非穩(wěn)定狀態(tài)下蔫骂，如增殖類細(xì)胞出現(xiàn)由于細(xì)胞周期相關(guān)基因的不同導(dǎo)致細(xì)胞聚類發(fā)生一定的偏移忧便。

其實(shí)在許多已經(jīng)發(fā)表的文獻(xiàn)中，我們也可以看到由于細(xì)胞周期不同所導(dǎo)致的分類偏移急灭。如在Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing 中 vCAFs和cCAFs中只有cell cycle genes的表達(dá)差異姐浮。

當(dāng)作者使用SC3重新進(jìn)行聚類時(shí)，兩種亞型又重新聚在了一起 (重點(diǎn)關(guān)注圖中的cluster1和cluster3葬馋，不要被顏色誤導(dǎo)了卖鲤。同一個(gè)cluster，不同的配色方案畴嘶，也是服了蛋逾，還是看看史上最全的圖表色彩運(yùn)用原理吧)。

image

系統(tǒng)介紹

通過計(jì)算各個(gè)細(xì)胞可能所處的細(xì)胞周期階段的得分 (G1,S,G2期）窗悯，并在預(yù)處理過程 (主要是ScaleData步驟)中將細(xì)胞周期得分作為混雜因素移除换怖，從而排除單細(xì)胞所處細(xì)胞周期不同對基因表達(dá)和細(xì)胞分型的影響。作者在小鼠造血祖細(xì)胞的數(shù)據(jù)集上證明了該觀點(diǎn) （Nestorowa et al. Blood 2016.）蟀瞧，其實(shí)只要是Seurat v3對象，自己的數(shù)據(jù)都是可以跑得通的条摸。

細(xì)胞周期相關(guān)基因集使用的是人的基因悦污，用小鼠進(jìn)行測試，說明該細(xì)胞周期相關(guān)基因數(shù)據(jù)集適合人和小鼠钉蒲；如果是其它物種切端，準(zhǔn)備方法見 https://github.com/satijalab/seurat/issues/462。

下面是操作代碼：

library(Seurat)

# 讀取表達(dá)矩陣顷啼， The first row is a header row, the first column is rownames
exp.mat <- read.table(file = "../data/nestorawa_forcellcycle_expressionMatrix.txt",
                      header = TRUE, as.is = TRUE, row.names = 1)

# 一系列的細(xì)胞周期相關(guān)的markers,其中包括處于S期的43個(gè)細(xì)胞周期相關(guān)基因踏枣，54個(gè)G2M期的細(xì)胞周期相關(guān)基因昌屉，
# from Tirosh et al, 2015, is loaded with Seurat.  We can
# segregate this list into markers of G2/M phase and markers of S phase
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes

# 創(chuàng)建Seurat對象并進(jìn)行標(biāo)準(zhǔn)化；marrow <- CreateSeuratObject(counts = exp.mat)
marrow <- NormalizeData(marrow)
marrow <- FindVariableFeatures(marrow, selection.method = "vst")
marrow <- ScaleData(marrow, features = rownames(marrow))

如果我們在Seurat對象上進(jìn)行PCA分析茵瀑，使用FindVariableFeatures中找到高可變基因间驮，進(jìn)行PCA分析，并展示對各個(gè)主成分貢獻(xiàn)最大的基因马昨。

PCA分析有不少需要注意的竞帽，具體見用了這么多年的PCA可視化竟然是錯(cuò)的！：枧酢屹篓！ 和 PCA主成分分析實(shí)戰(zhàn)和可視化 | 附R代碼和測試數(shù)據(jù)。

marrow <- RunPCA(marrow, features = VariableFeatures(marrow), ndims.print = 6:10, nfeatures.print = 10)

我們看到對PC8和PC10貢獻(xiàn)最大的基因中一部分是細(xì)胞周期相關(guān)基因匙奴，包括TOP2A和MKI67堆巧。

image

從下面的熱圖也可以看出來。

DimHeatmap(marrow, dims = c(8, 10))#熱圖表示

image

我們將嘗試從數(shù)據(jù)中去除該組分泼菌，從而確保細(xì)胞周期異質(zhì)性對PCA或下游分析沒有貢獻(xiàn)谍肤。

分配細(xì)胞周期分?jǐn)?shù)

首先，根據(jù)其G2/M和S期標(biāo)記基因的表達(dá)為每個(gè)細(xì)胞分配一個(gè)所處周期的分?jǐn)?shù)灶轰。這些標(biāo)記基因的表達(dá)水平應(yīng)該是反相關(guān)的谣沸，而不表達(dá)這些標(biāo)記基因的細(xì)胞可能處于G1期。

我們用CellCycleScoring函數(shù)計(jì)算細(xì)胞周期分?jǐn)?shù)笋颤，并在metadata中存儲(chǔ)S和G2/M分?jǐn)?shù)乳附，以及G2M，S或G1階段中每個(gè)細(xì)胞的預(yù)測分類伴澄。如果指定set.ident=T赋除，則CellCycleScoring將Seurat對象中每個(gè)細(xì)胞的分組信息設(shè)置為其所處的細(xì)胞周期階段。

marrow <- CellCycleScoring(marrow, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)

# view cell cycle scores and phase assignments
head(marrow[[]])

image

觀察細(xì)胞周期基因的表達(dá)情況 （評估計(jì)算的準(zhǔn)確性）

注：R語言可視化學(xué)習(xí)筆記之ggridges包可繪制類似圖形非凌。

# 觀察細(xì)胞周期基因的表達(dá)情況
RidgePlot(marrow, features = c("PCNA", "TOP2A", "MCM6", "MKI67"), ncol = 2)

image

我們用CellCycleScoring函數(shù)計(jì)算細(xì)胞周期分?jǐn)?shù)举农，并在metadata中存儲(chǔ)S和G2/M分?jǐn)?shù)，以及G2M敞嗡，S或G1階段中每個(gè)細(xì)胞的預(yù)測分類颁糟。如果指定set.ident=T，則CellCycleScoring將Seurat對象中每個(gè)細(xì)胞的分組信息設(shè)置為其所處的細(xì)胞周期階段喉悴。利用細(xì)胞周期基因進(jìn)行PCA分析 （這個(gè)例子可以拓展棱貌，使用任意指定的基因集進(jìn)行細(xì)胞周期分析）

PCA分析有不少需要注意的，具體見用了這么多年的PCA可視化竟然是錯(cuò)的；唷；橥选！ 和 PCA主成分分析實(shí)戰(zhàn)和可視化 | 附R代碼和測試數(shù)據(jù)。

# Running a PCA on cell cycle genes reveals, unsurprisingly, that cells separate entirely by
# phase
marrow <- RunPCA(marrow, features = c(s.genes, g2m.genes))
DimPlot(marrow)

image

在數(shù)據(jù)標(biāo)歸一化時(shí)去除細(xì)胞周期影響

marrow <- ScaleData(marrow, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(marrow))

再次做PCA時(shí)障贸，就看不到細(xì)胞周期相關(guān)基因?qū)χ鞒煞值呢暙I(xiàn)了错森。

# 我們可以看到組分中不再出現(xiàn)細(xì)胞周期相關(guān)的基因
marrow <- RunPCA(marrow, features = VariableFeatures(marrow), nfeatures.print = 10)

image

再次根據(jù)細(xì)胞周期相關(guān)基因進(jìn)行PCA分析時(shí)，也不分不出群了篮洁，說明移除細(xì)胞周期影響的效果還是比較好的涩维。

# When running a PCA on only cell cycle genes, cells no longer separate by cell-cycle phase
marrow <- RunPCA(marrow, features = c(s.genes, g2m.genes))
DimPlot(marrow)

image

如果細(xì)胞周期不合適時(shí)怎么辦？

上述過程去除了與細(xì)胞周期相關(guān)的所有信息嘀粱。在某些情況下激挪，我們發(fā)現(xiàn)這會(huì)對下游分析產(chǎn)生負(fù)面影響，特別是在分化過程（如小鼠血細(xì)胞分化生成過程中 hematopoiesis）中锋叨，干細(xì)胞處于靜止?fàn)顟B(tài)垄分，分化細(xì)胞正在增殖（反之亦然）。在這種情況下娃磺，消除所有細(xì)胞周期效應(yīng)也會(huì)模糊干細(xì)胞和前體細(xì)胞之間的區(qū)別薄湿。

作為替代方案，我們建議消除G2M和S期分?jǐn)?shù)之間的差異偷卧。這意味著將保持非周期細(xì)胞和周期細(xì)胞的組分差異豺瘤，但是增殖細(xì)胞之間的細(xì)胞周期階段的差異將從數(shù)據(jù)中去除。

# 計(jì)算并移除分?jǐn)?shù)差異
marrow$CC.Difference <- marrow$S.Score - marrow$G2M.Score
marrow <- ScaleData(marrow, vars.to.regress = "CC.Difference", features = rownames(marrow))

細(xì)胞周期相關(guān)的基因?qū)Ω鱾€(gè)主成分貢獻(xiàn)減小

# cell cycle effects strongly mitigated in PCA
# 在PCA中不再出現(xiàn)大量的細(xì)胞周期相關(guān)的基因
marrow <- RunPCA(marrow, features = VariableFeatures(marrow), nfeatures.print = 10)

image

G1期區(qū)分開听诸，G2/M和S期聚在一起坐求。

# 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(marrow, features = c(s.genes, g2m.genes))
DimPlot(marrow)

image

單細(xì)胞是目前很火的領(lǐng)域，分析工具很多晌梨，而且也還在不斷發(fā)展中桥嗤。Seurat是其中一個(gè)，雖然好用仔蝌，卻不一定是最好的泛领。而且運(yùn)用好工具，需要對原理有很好的理解敛惊，易生信8月份的單細(xì)胞課程邀請來在單細(xì)胞算法開發(fā)上很有經(jīng)驗(yàn)的中科院博士開課渊鞋，深入淺出的講述了單細(xì)胞分析的方法和注意事項(xiàng)，我個(gè)人認(rèn)為課程講的特別好瞧挤，很多學(xué)員也反映特別好锡宋，在此強(qiáng)烈推薦。