本文是參考學(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)行Ｖ莸堋！低零！

細(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)行ｌＦ酢！
#獲取基因在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
# 圖片