單細(xì)胞轉(zhuǎn)錄組基礎(chǔ)分析三:降維與聚類(lèi)

本文是參考學(xué)習(xí)單細(xì)胞轉(zhuǎn)錄組基礎(chǔ)分析三:降維與聚類(lèi)
的學(xué)習(xí)筆記拐邪。可能根據(jù)學(xué)習(xí)情況有所改動(dòng)。

尋找高變基因Seurat負(fù)責(zé)篩選高變基因的函數(shù)是FindVariableFeatures()拢操,它并不刪除scRNA對(duì)象中的非高變基因忽匈。找出的結(jié)果可以通過(guò)VariableFeatures()函數(shù)獲取恐仑。

library(Seurat)

選擇的高變基因可通過(guò)cluster/VariableFeatures.png文件查看嗤攻,橫坐標(biāo)是某基因在所有細(xì)胞中的平均表達(dá)值集索,縱坐標(biāo)是此基因的方差泳赋。紅色的點(diǎn)是被選中的高變基因雌桑,黑色的點(diǎn)是未被選中的基因,變異程度最高的10個(gè)基因在如圖中標(biāo)注了基因名稱(chēng)祖今。

圖片

在進(jìn)行PCA降維之前還有要對(duì)數(shù)據(jù)進(jìn)行中心化(必選)和細(xì)胞周期回歸分析(可選)校坑,它們可以使用ScaleData()函數(shù)一步完成。這兩項(xiàng)分析雖然可以使用一個(gè)函數(shù)完成千诬,但是我覺(jué)得有必要分開(kāi)講一下它們的原理與作用耍目。

數(shù)據(jù)中心化數(shù)據(jù)中心化是PCA分析和一些可視化步驟的必須要求,它使數(shù)據(jù)具有統(tǒng)一的尺度(scale)徐绑,其運(yùn)行效果示意圖如下所示:

圖片

左邊的圖是原始數(shù)據(jù)畫(huà)的散點(diǎn)圖邪驮,右邊的圖是用scale之后的數(shù)據(jù)畫(huà)的。細(xì)心的同學(xué)可能發(fā)現(xiàn)了傲茄,這兩張圖點(diǎn)的相對(duì)位置沒(méi)變毅访,但是坐標(biāo)軸的尺度改變了。Scaling代碼如下:

##如果內(nèi)存足夠最好對(duì)所有基因進(jìn)行中心化

scRNA對(duì)象中原始表達(dá)矩陣經(jīng)過(guò)標(biāo)準(zhǔn)化和中心化之后盘榨,已經(jīng)產(chǎn)生了三套基因表達(dá)數(shù)據(jù)喻粹,可以通過(guò)以下命令獲得:

#原始表達(dá)矩陣

細(xì)胞周期回歸
上一步找到的高變基因,常常會(huì)包含一些細(xì)胞周期相關(guān)基因草巡。它們會(huì)導(dǎo)致細(xì)胞聚類(lèi)發(fā)生一定的偏移守呜,即相同類(lèi)型的細(xì)胞在聚類(lèi)時(shí)會(huì)因?yàn)榧?xì)胞周期的不同而分開(kāi)。如在Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing 中 vCAFs和cCAFs中只有cell cycle genes的表達(dá)差異。查看我們選擇的高變基因中有哪些細(xì)胞周期相關(guān)基因:

CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(scRNA))
圖片

細(xì)胞周期評(píng)分

g2m_genes = cc.genes$g2m.genes

以上代碼運(yùn)行之后會(huì)在scRNA@meta.data中添加S.Score查乒、G2M.Score和Phase三列有關(guān)細(xì)胞周期的信息弥喉。

查看細(xì)胞周期基因?qū)?xì)胞聚類(lèi)的影響

scRNAa <- RunPCA(scRNA, features = c(s_genes, g2m_genes))
圖片

右圖是seurat細(xì)胞周期回歸分析教程的案例,提示這種情況需要剔除細(xì)胞周期對(duì)聚類(lèi)的影響玛迄。左圖是我們這次分析的pbmc的結(jié)果档桃,看情況沒(méi)必要做細(xì)胞周期回歸分析。

##如果需要消除細(xì)胞周期的影響

PCA降維并提取主成分****PCA降維

scRNA <- RunPCA(scRNA, features = VariableFeatures(scRNA)) 

運(yùn)行之后會(huì)在cluster文件夾下面得到pca圖憔晒,如下所示:

圖片

左圖是根據(jù)主成分1和2的值將細(xì)胞在平面上展示出來(lái)藻肄,右圖展示前20個(gè)主成分的解釋量(pca中等同于標(biāo)準(zhǔn)差)。后續(xù)分析要根據(jù)右圖選擇提取的pc軸數(shù)量拒担,一般選擇斜率平滑的點(diǎn)之前的所有pc軸嘹屯,此圖我的建議是選擇前18個(gè)pc軸。

pc.num=1:18

**獲取PCA結(jié)果 **

此部分代碼為分析非必須代碼从撼,不建議運(yùn)行V莸堋!低零!

細(xì)胞聚類(lèi)此步利用 細(xì)胞-PC值矩陣計(jì)算細(xì)胞之間的距離婆翔,然后利用距離矩陣來(lái)聚類(lèi)。其中有兩個(gè)參數(shù)需要人工選擇掏婶,第一個(gè)是FindNeighbors()函數(shù)中的dims參數(shù)啃奴,需要指定哪些pc軸用于分析,選擇依據(jù)是之前介紹的cluster/pca.png文件中的右圖雄妥。第二個(gè)是FindClusters()函數(shù)中的resolution參數(shù)最蕾,需要指定0.1-0.9之間的一個(gè)數(shù)值,用于決定clusters的相對(duì)數(shù)量老厌,數(shù)值越大cluters越多瘟则。

scRNA <- FindNeighbors(scRNA, dims = pc.num) 

運(yùn)行之后會(huì)在rstudio中顯示每個(gè)cluster的細(xì)胞數(shù)量,cluter文件夾中也會(huì)得到細(xì)胞-cluter的映射文件cell_cluster.csv

圖片

非線性降維此步也是利用 細(xì)胞-PC值矩陣作為輸入數(shù)據(jù)

#tSNE

運(yùn)行之后會(huì)得到TSNE圖和UMAP圖枝秤,以及細(xì)胞在這兩張圖上的坐標(biāo)值文件embed_tsne.csv和embed_umap.csv

圖片
# 尋找高變基因
# Seurat負(fù)責(zé)篩選高變基因的函數(shù)是FindVariableFeatures()醋拧,它并不刪除scRNA對(duì)象中的非高變基因。找出的結(jié)果可以通過(guò)VariableFeatures()函數(shù)獲取淀弹。
library(Seurat)
library(tidyverse)
library(patchwork)
rm(list=ls())
dir.create("cluster")
scRNA <- readRDS("scRNA.rds")
scRNA <- FindVariableFeatures(scRNA, selection.method = "vst", nfeatures = 2000) 
top10 <- head(VariableFeatures(scRNA), 10) 
plot1 <- VariableFeaturePlot(scRNA) 
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) 
plot <- CombinePlots(plots = list(plot1, plot2),legend="bottom") 
ggsave("cluster/VariableFeatures.pdf", plot = plot, width = 8, height = 6) 
ggsave("cluster/VariableFeatures.png", plot = plot, width = 8, height = 6)
# 選擇的高變基因可通過(guò)cluster/VariableFeatures.png文件查看丹壕,橫坐標(biāo)是某基因在所有細(xì)胞中的平均表達(dá)值,縱坐標(biāo)是此基因的方差垦页。紅色的點(diǎn)是被選中的高變基因雀费,黑色的點(diǎn)是未被選中的基因干奢,變異程度最高的10個(gè)基因在如圖中標(biāo)注了基因名稱(chēng)痊焊。
# 
# 
# 圖片
# 
# 
# 在進(jìn)行PCA降維之前還有要對(duì)數(shù)據(jù)進(jìn)行中心化(必選)和細(xì)胞周期回歸分析(可選),它們可以使用ScaleData()函數(shù)一步完成。這兩項(xiàng)分析雖然可以使用一個(gè)函數(shù)完成薄啥,但是我覺(jué)得有必要分開(kāi)講一下它們的原理與作用辕羽。
# 
# 
# 數(shù)據(jù)中心化
# 數(shù)據(jù)中心化是PCA分析和一些可視化步驟的必須要求,它使數(shù)據(jù)具有統(tǒng)一的尺度(scale)垄惧,其運(yùn)行效果示意圖如下所示:
# 圖片
# 左邊的圖是原始數(shù)據(jù)畫(huà)的散點(diǎn)圖刁愿,右邊的圖是用scale之后的數(shù)據(jù)畫(huà)的。細(xì)心的同學(xué)可能發(fā)現(xiàn)了到逊,這兩張圖點(diǎn)的相對(duì)位置沒(méi)變铣口,但是坐標(biāo)軸的尺度改變了。
# Scaling代碼如下:
##如果內(nèi)存足夠最好對(duì)所有基因進(jìn)行中心化
scale.genes <-  rownames(scRNA)
scRNA <- ScaleData(scRNA, features = scale.genes)
##如果內(nèi)存不夠觉壶,可以只對(duì)高變基因進(jìn)行標(biāo)準(zhǔn)化
scale.genes <-  VariableFeatures(scRNA)
scRNA <- ScaleData(scRNA, features = scale.genes)
# scRNA對(duì)象中原始表達(dá)矩陣經(jīng)過(guò)標(biāo)準(zhǔn)化和中心化之后脑题,已經(jīng)產(chǎn)生了三套基因表達(dá)數(shù)據(jù),可以通過(guò)以下命令獲得:
#原始表達(dá)矩陣
GetAssayData(scRNA,slot="counts",assay="RNA") 
#標(biāo)準(zhǔn)化之后的表達(dá)矩陣              
GetAssayData(scRNA,slot="data",assay="RNA")
#中心化之后的表達(dá)矩陣 
GetAssayData(scRNA,slot="scale.data",assay="RNA") 

# 細(xì)胞周期回歸
# 上一步找到的高變基因铜靶,常常會(huì)包含一些細(xì)胞周期相關(guān)基因叔遂。它們會(huì)導(dǎo)致細(xì)胞聚類(lèi)發(fā)生一定的偏移,即相同類(lèi)型的細(xì)胞在聚類(lèi)時(shí)會(huì)因?yàn)榧?xì)胞周期的不同而分開(kāi)争剿。如在Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing 中 vCAFs和cCAFs中只有cell cycle genes的表達(dá)差異已艰。
# 查看我們選擇的高變基因中有哪些細(xì)胞周期相關(guān)基因:
CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(scRNA))
# 圖片


# 細(xì)胞周期評(píng)分
g2m_genes = cc.genes$g2m.genes
g2m_genes = CaseMatch(search = g2m_genes, match = rownames(scRNA))
s_genes = cc.genes$s.genes
s_genes = CaseMatch(search = s_genes, match = rownames(scRNA))
scRNA <- CellCycleScoring(object=scRNA,  g2m.features=g2m_genes,  s.features=s_genes)
# 以上代碼運(yùn)行之后會(huì)在scRNA@meta.data中添加S.Score、G2M.Score和Phase三列有關(guān)細(xì)胞周期的信息蚕苇。



# 查看細(xì)胞周期基因?qū)?xì)胞聚類(lèi)的影響

scRNAa <- RunPCA(scRNA, features = c(s_genes, g2m_genes))
p <- DimPlot(scRNAa, reduction = "pca", group.by = "Phase")
ggsave("cluster/cellcycle_pca.png", p, width = 8, height = 6)
# 圖片

# 右圖是seurat細(xì)胞周期回歸分析教程的案例哩掺,提示這種情況需要剔除細(xì)胞周期對(duì)聚類(lèi)的影響。左圖是我們這次分析的pbmc的結(jié)果涩笤,看情況沒(méi)必要做細(xì)胞周期回歸分析疮丛。
##如果需要消除細(xì)胞周期的影響
#scRNAb <- ScaleData(scRNA, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(scRNA))

# PCA降維并提取主成分
# PCA降維
scRNA <- RunPCA(scRNA, features = VariableFeatures(scRNA)) 
plot1 <- DimPlot(scRNA, reduction = "pca", group.by="orig.ident") 
plot2 <- ElbowPlot(scRNA, ndims=20, reduction="pca") 
plot1+plot2
plotc <- plot1+plot2
ggsave("cluster/pca.pdf", plot = plotc, width = 8, height = 4) 
ggsave("cluster/pca.png", plot = plotc, width = 8, height = 4)
# 運(yùn)行之后會(huì)在cluster文件夾下面得到pca圖,如下所示:
# 圖片
# 
# 左圖是根據(jù)主成分1和2的值將細(xì)胞在平面上展示出來(lái)辆它,右圖展示前20個(gè)主成分的解釋量(pca中等同于標(biāo)準(zhǔn)差)誊薄。后續(xù)分析要根據(jù)右圖選擇提取的pc軸數(shù)量,一般選擇斜率平滑的點(diǎn)之前的所有pc軸锰茉,此圖我的建議是選擇前18個(gè)pc軸呢蔫。
pc.num=1:18


# 獲取PCA結(jié)果 
# 此部分代碼為分析非必須代碼,不建議運(yùn)行lF酢!
#獲取基因在pc軸上的投射值
Loadings(object = scRNA[["pca"]])
#獲取各個(gè)細(xì)胞的pc值
Embeddings(object = scRNA[["pca"]])
#獲取各pc軸解釋量方差
Stdev(scRNA)
#查看決定pc值的top10基因协屡, 此例查看pc1-pc5軸
print(scRNA[["pca"]], dims = 1:5, nfeatures = 10) 
#查看決定pc值的top10基因在500個(gè)細(xì)胞中的熱圖俏脊,此例查看pc1-pc9軸
DimHeatmap(scRNA, dims = 1:9, nfeatures=10, cells = 500, balanced = TRUE)

# 細(xì)胞聚類(lèi)
# 此步利用 細(xì)胞-PC值 矩陣計(jì)算細(xì)胞之間的距離,然后利用距離矩陣來(lái)聚類(lèi)
# 肤晓。其中有兩個(gè)參數(shù)需要人工選擇爷贫,
# 第一個(gè)是FindNeighbors()函數(shù)中的dims參數(shù)认然,需要指定哪些pc軸用于分析,
# 選擇依據(jù)是之前介紹的cluster/pca.png文件中的右圖漫萄。
# 第二個(gè)是FindClusters()函數(shù)中的resolution參數(shù)卷员,
# 需要指定0.1-0.9之間的一個(gè)數(shù)值,用于決定clusters的相對(duì)數(shù)量腾务,數(shù)值越大cluters越多毕骡。
scRNA <- FindNeighbors(scRNA, dims = pc.num) 


scRNA <- FindClusters(scRNA, resolution = 0.5)
table(scRNA@meta.data$seurat_clusters)
metadata <- scRNA@meta.data
cell_cluster <- data.frame(cell_ID=rownames(metadata), cluster_ID=metadata$seurat_clusters)
write.csv(cell_cluster,'cluster/cell_cluster.csv',row.names = F)
# 運(yùn)行之后會(huì)在rstudio中顯示每個(gè)cluster的細(xì)胞數(shù)量,cluter文件夾中也會(huì)得到細(xì)胞-cluter的映射文件cell_cluster.csv
# 圖片
# 
# 
# 非線性降維
# 此步也是利用 細(xì)胞-PC值 矩陣作為輸入數(shù)據(jù)
#tSNE
scRNA = RunTSNE(scRNA, dims = pc.num)
embed_tsne <- Embeddings(scRNA, 'tsne')
write.csv(embed_tsne,'cluster/embed_tsne.csv')
plot1 = DimPlot(scRNA, reduction = "tsne") 
ggsave("cluster/tSNE.pdf", plot = plot1, width = 8, height = 7)
ggsave("cluster/tSNE.png", plot = plot1, width = 8, height = 7)
#UMAP
scRNA <- RunUMAP(scRNA, dims = pc.num)
embed_umap <- Embeddings(scRNA, 'umap')
write.csv(embed_umap,'cluster/embed_umap.csv') 
plot2 = DimPlot(scRNA, reduction = "umap") 
ggsave("cluster/UMAP.pdf", plot = plot2, width = 8, height = 7)
ggsave("cluster/UMAP.png", plot = plot2, width = 8, height = 7)
#合并tSNE與UMAP
plotc <- plot1+plot2+ plot_layout(guides = 'collect')
ggsave("cluster/tSNE_UMAP.pdf", plot = plotc, width = 10, height = 5)
ggsave("cluster/tSNE_UMAP.png", plot = plotc, width = 10, height = 5)
##保存數(shù)據(jù)
saveRDS(scRNA, file="scRNA.rds")
# 運(yùn)行之后會(huì)得到TSNE圖和UMAP圖岩瘦,以及細(xì)胞在這兩張圖上的坐標(biāo)值文件embed_tsne.csv和embed_umap.csv
# 圖片
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
  • 序言:七十年代末未巫,一起剝皮案震驚了整個(gè)濱河市,隨后出現(xiàn)的幾起案子启昧,更是在濱河造成了極大的恐慌橱赠,老刑警劉巖,帶你破解...
    沈念sama閱讀 222,590評(píng)論 6 517
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件箫津,死亡現(xiàn)場(chǎng)離奇詭異狭姨,居然都是意外死亡,警方通過(guò)查閱死者的電腦和手機(jī)苏遥,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 95,157評(píng)論 3 399
  • 文/潘曉璐 我一進(jìn)店門(mén)饼拍,熙熙樓的掌柜王于貴愁眉苦臉地迎上來(lái),“玉大人田炭,你說(shuō)我怎么就攤上這事师抄。” “怎么了教硫?”我有些...
    開(kāi)封第一講書(shū)人閱讀 169,301評(píng)論 0 362
  • 文/不壞的土叔 我叫張陵叨吮,是天一觀的道長(zhǎng)。 經(jīng)常有香客問(wèn)我瞬矩,道長(zhǎng)茶鉴,這世上最難降的妖魔是什么? 我笑而不...
    開(kāi)封第一講書(shū)人閱讀 60,078評(píng)論 1 300
  • 正文 為了忘掉前任景用,我火速辦了婚禮涵叮,結(jié)果婚禮上,老公的妹妹穿的比我還像新娘伞插。我一直安慰自己割粮,他們只是感情好,可當(dāng)我...
    茶點(diǎn)故事閱讀 69,082評(píng)論 6 398
  • 文/花漫 我一把揭開(kāi)白布媚污。 她就那樣靜靜地躺著舀瓢,像睡著了一般。 火紅的嫁衣襯著肌膚如雪耗美。 梳的紋絲不亂的頭發(fā)上京髓,一...
    開(kāi)封第一講書(shū)人閱讀 52,682評(píng)論 1 312
  • 那天航缀,我揣著相機(jī)與錄音,去河邊找鬼朵锣。 笑死,一個(gè)胖子當(dāng)著我的面吹牛甸私,可吹牛的內(nèi)容都是我干的诚些。 我是一名探鬼主播,決...
    沈念sama閱讀 41,155評(píng)論 3 422
  • 文/蒼蘭香墨 我猛地睜開(kāi)眼皇型,長(zhǎng)吁一口氣:“原來(lái)是場(chǎng)噩夢(mèng)啊……” “哼诬烹!你這毒婦竟也來(lái)了?” 一聲冷哼從身側(cè)響起弃鸦,我...
    開(kāi)封第一講書(shū)人閱讀 40,098評(píng)論 0 277
  • 序言:老撾萬(wàn)榮一對(duì)情侶失蹤绞吁,失蹤者是張志新(化名)和其女友劉穎,沒(méi)想到半個(gè)月后唬格,有當(dāng)?shù)厝嗽跇?shù)林里發(fā)現(xiàn)了一具尸體家破,經(jīng)...
    沈念sama閱讀 46,638評(píng)論 1 319
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡,尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 38,701評(píng)論 3 342
  • 正文 我和宋清朗相戀三年购岗,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了汰聋。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點(diǎn)故事閱讀 40,852評(píng)論 1 353
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡喊积,死狀恐怖烹困,靈堂內(nèi)的尸體忽然破棺而出,到底是詐尸還是另有隱情乾吻,我是刑警寧澤髓梅,帶...
    沈念sama閱讀 36,520評(píng)論 5 351
  • 正文 年R本政府宣布,位于F島的核電站绎签,受9級(jí)特大地震影響枯饿,放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜诡必,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 42,181評(píng)論 3 335
  • 文/蒙蒙 一鸭你、第九天 我趴在偏房一處隱蔽的房頂上張望。 院中可真熱鬧擒权,春花似錦袱巨、人聲如沸。這莊子的主人今日做“春日...
    開(kāi)封第一講書(shū)人閱讀 32,674評(píng)論 0 25
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽(yáng)。三九已至剖效,卻和暖如春嫉入,著一層夾襖步出監(jiān)牢的瞬間焰盗,已是汗流浹背。 一陣腳步聲響...
    開(kāi)封第一講書(shū)人閱讀 33,788評(píng)論 1 274
  • 我被黑心中介騙來(lái)泰國(guó)打工咒林, 沒(méi)想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留熬拒,地道東北人。 一個(gè)月前我還...
    沈念sama閱讀 49,279評(píng)論 3 379
  • 正文 我出身青樓垫竞,卻偏偏與公主長(zhǎng)得像澎粟,于是被迫代替她去往敵國(guó)和親。 傳聞我的和親對(duì)象是個(gè)殘疾皇子欢瞪,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 45,851評(píng)論 2 361

推薦閱讀更多精彩內(nèi)容