单细胞转录组基础分析三：降维与聚类

本文是参考学习单细胞转录组基础分析三：降维与聚类
的学习笔记。可能根据学习情况有所改动。

寻找高变基因Seurat负责筛选高变基因的函数是FindVariableFeatures()，它并不删除scRNA对象中的非高变基因。找出的结果可以通过VariableFeatures()函数获取。

library(Seurat)

选择的高变基因可通过cluster/VariableFeatures.png文件查看，横坐标是某基因在所有细胞中的平均表达值，纵坐标是此基因的方差。红色的点是被选中的高变基因，黑色的点是未被选中的基因，变异程度最高的10个基因在如图中标注了基因名称。

图片

在进行PCA降维之前还有要对数据进行中心化（必选）和细胞周期回归分析（可选），它们可以使用ScaleData()函数一步完成。这两项分析虽然可以使用一个函数完成，但是我觉得有必要分开讲一下它们的原理与作用。

数据中心化数据中心化是PCA分析和一些可视化步骤的必须要求，它使数据具有统一的尺度（scale），其运行效果示意图如下所示：

图片

左边的图是原始数据画的散点图，右边的图是用scale之后的数据画的。细心的同学可能发现了，这两张图点的相对位置没变，但是坐标轴的尺度改变了。Scaling代码如下：

##如果内存足够最好对所有基因进行中心化

scRNA对象中原始表达矩阵经过标准化和中心化之后，已经产生了三套基因表达数据，可以通过以下命令获得：

#原始表达矩阵

细胞周期回归
上一步找到的高变基因，常常会包含一些细胞周期相关基因。它们会导致细胞聚类发生一定的偏移，即相同类型的细胞在聚类时会因为细胞周期的不同而分开。如在Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing 中 vCAFs和cCAFs中只有cell cycle genes的表达差异。查看我们选择的高变基因中有哪些细胞周期相关基因：

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

图片

细胞周期评分

g2m_genes = cc.genes$g2m.genes

以上代码运行之后会在scRNA@meta.data中添加S.Score、G2M.Score和Phase三列有关细胞周期的信息。

查看细胞周期基因对细胞聚类的影响

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

图片

右图是seurat细胞周期回归分析教程的案例，提示这种情况需要剔除细胞周期对聚类的影响。左图是我们这次分析的pbmc的结果，看情况没必要做细胞周期回归分析。

##如果需要消除细胞周期的影响

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

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

运行之后会在cluster文件夹下面得到pca图，如下所示：

图片

左图是根据主成分1和2的值将细胞在平面上展示出来，右图展示前20个主成分的解释量（pca中等同于标准差）。后续分析要根据右图选择提取的pc轴数量，一般选择斜率平滑的点之前的所有pc轴，此图我的建议是选择前18个pc轴。

pc.num=1:18

**获取PCA结果 **

此部分代码为分析非必须代码，不建议运行！！！

细胞聚类此步利用 细胞-PC值矩阵计算细胞之间的距离，然后利用距离矩阵来聚类。其中有两个参数需要人工选择，第一个是FindNeighbors()函数中的dims参数，需要指定哪些pc轴用于分析，选择依据是之前介绍的cluster/pca.png文件中的右图。第二个是FindClusters()函数中的resolution参数，需要指定0.1-0.9之间的一个数值，用于决定clusters的相对数量，数值越大cluters越多。

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

运行之后会在rstudio中显示每个cluster的细胞数量，cluter文件夹中也会得到细胞-cluter的映射文件cell_cluster.csv

图片

非线性降维此步也是利用 细胞-PC值矩阵作为输入数据

#tSNE

运行之后会得到TSNE图和UMAP图，以及细胞在这两张图上的坐标值文件embed_tsne.csv和embed_umap.csv

图片

# 寻找高变基因
# Seurat负责筛选高变基因的函数是FindVariableFeatures()，它并不删除scRNA对象中的非高变基因。找出的结果可以通过VariableFeatures()函数获取。
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)
# 选择的高变基因可通过cluster/VariableFeatures.png文件查看，横坐标是某基因在所有细胞中的平均表达值，纵坐标是此基因的方差。红色的点是被选中的高变基因，黑色的点是未被选中的基因，变异程度最高的10个基因在如图中标注了基因名称。
# 
# 
# 图片
# 
# 
# 在进行PCA降维之前还有要对数据进行中心化（必选）和细胞周期回归分析（可选），它们可以使用ScaleData()函数一步完成。这两项分析虽然可以使用一个函数完成，但是我觉得有必要分开讲一下它们的原理与作用。
# 
# 
# 数据中心化
# 数据中心化是PCA分析和一些可视化步骤的必须要求，它使数据具有统一的尺度（scale），其运行效果示意图如下所示：
# 图片
# 左边的图是原始数据画的散点图，右边的图是用scale之后的数据画的。细心的同学可能发现了，这两张图点的相对位置没变，但是坐标轴的尺度改变了。
# Scaling代码如下：
##如果内存足够最好对所有基因进行中心化
scale.genes <-  rownames(scRNA)
scRNA <- ScaleData(scRNA, features = scale.genes)
##如果内存不够，可以只对高变基因进行标准化
scale.genes <-  VariableFeatures(scRNA)
scRNA <- ScaleData(scRNA, features = scale.genes)
# scRNA对象中原始表达矩阵经过标准化和中心化之后，已经产生了三套基因表达数据，可以通过以下命令获得：
#原始表达矩阵
GetAssayData(scRNA,slot="counts",assay="RNA") 
#标准化之后的表达矩阵              
GetAssayData(scRNA,slot="data",assay="RNA")
#中心化之后的表达矩阵 
GetAssayData(scRNA,slot="scale.data",assay="RNA") 

# 细胞周期回归
# 上一步找到的高变基因，常常会包含一些细胞周期相关基因。它们会导致细胞聚类发生一定的偏移，即相同类型的细胞在聚类时会因为细胞周期的不同而分开。如在Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing 中 vCAFs和cCAFs中只有cell cycle genes的表达差异。
# 查看我们选择的高变基因中有哪些细胞周期相关基因：
CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(scRNA))
# 图片


# 细胞周期评分
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)
# 以上代码运行之后会在scRNA@meta.data中添加S.Score、G2M.Score和Phase三列有关细胞周期的信息。



# 查看细胞周期基因对细胞聚类的影响

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细胞周期回归分析教程的案例，提示这种情况需要剔除细胞周期对聚类的影响。左图是我们这次分析的pbmc的结果，看情况没必要做细胞周期回归分析。
##如果需要消除细胞周期的影响
#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)
# 运行之后会在cluster文件夹下面得到pca图，如下所示：
# 图片
# 
# 左图是根据主成分1和2的值将细胞在平面上展示出来，右图展示前20个主成分的解释量（pca中等同于标准差）。后续分析要根据右图选择提取的pc轴数量，一般选择斜率平滑的点之前的所有pc轴，此图我的建议是选择前18个pc轴。
pc.num=1:18


# 获取PCA结果 
# 此部分代码为分析非必须代码，不建议运行！！！
#获取基因在pc轴上的投射值
Loadings(object = scRNA[["pca"]])
#获取各个细胞的pc值
Embeddings(object = scRNA[["pca"]])
#获取各pc轴解释量方差
Stdev(scRNA)
#查看决定pc值的top10基因， 此例查看pc1-pc5轴
print(scRNA[["pca"]], dims = 1:5, nfeatures = 10) 
#查看决定pc值的top10基因在500个细胞中的热图，此例查看pc1-pc9轴
DimHeatmap(scRNA, dims = 1:9, nfeatures=10, cells = 500, balanced = TRUE)

# 细胞聚类
# 此步利用 细胞-PC值 矩阵计算细胞之间的距离，然后利用距离矩阵来聚类
# 。其中有两个参数需要人工选择，
# 第一个是FindNeighbors()函数中的dims参数，需要指定哪些pc轴用于分析，
# 选择依据是之前介绍的cluster/pca.png文件中的右图。
# 第二个是FindClusters()函数中的resolution参数，
# 需要指定0.1-0.9之间的一个数值，用于决定clusters的相对数量，数值越大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)
# 运行之后会在rstudio中显示每个cluster的细胞数量，cluter文件夹中也会得到细胞-cluter的映射文件cell_cluster.csv
# 图片
# 
# 
# 非线性降维
# 此步也是利用 细胞-PC值 矩阵作为输入数据
#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)
##保存数据
saveRDS(scRNA, file="scRNA.rds")
# 运行之后会得到TSNE图和UMAP图，以及细胞在这两张图上的坐标值文件embed_tsne.csv和embed_umap.csv
# 图片