本文是參考學(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)行lF酢!
#獲取基因在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
# 圖片