斯坦福大學(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)用原理吧)。

系統(tǒng)介紹

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

細(xì)胞周期相關(guān)基因集使用的是人的基因，用小鼠進(jìn)行測(cè)試，說(shuō)明該細(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對(duì)象并進(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))

如果我們?cè)赟eurat對(duì)象上進(jìn)行PCA分析，使用FindVariableFeatures中找到高可變基因，進(jìn)行PCA分析，并展示對(duì)各個(gè)主成分貢獻(xiàn)最大的基因。

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

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

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

從下面的熱圖也可以看出來(lái)。

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

我們將嘗試從數(shù)據(jù)中去除該組分，從而確保細(xì)胞周期異質(zhì)性對(duì)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ù)測(cè)分類。如果指定set.ident=T，則CellCycleScoring將Seurat對(duì)象中每個(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[[]])

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

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

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

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

PCA分析有不少需要注意的，具體見用了這么多年的PCA可視化竟然是錯(cuò)的！！！?和?PCA主成分分析實(shí)戰(zhàn)和可視化 | 附R代碼和測(cè)試數(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)

在數(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)

再次根據(jù)細(xì)胞周期相關(guān)基因進(jìn)行PCA分析時(shí)，也不分不出群了，說(shuō)明移除細(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)

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

上述過(guò)程去除了與細(xì)胞周期相關(guān)的所有信息。在某些情況下，我們發(fā)現(xiàn)這會(huì)對(duì)下游分析產(chǎn)生負(fù)面影響，特別是在分化過(guò)程（如小鼠血細(xì)胞分化生成過(guò)程中 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)

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)

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

總結(jié)

以上是生活随笔為你收集整理的Seurat亮点之细胞周期评分和回归的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。