scRNA基础分析-2：降维与聚类

scRNA基础分析-1:安装包、导入数据、过滤质控 - 简书
scRNA基础分析-2：降维与聚类 - 简书
scRNA基础分析-3：鉴定细胞类型 - 简书
scRNA基础分析-4：细胞亚类再聚类、注释 - 简书
scRNA基础分析-5：伪时间分析 - 简书
scRNA基础分析-6：富集分析 - 简书

library(Seurat)
library(tidyverse)
library(patchwork)
rm(list=ls())
scRNA <- readRDS("scRNA.rds")

1、首先选择2000个（默认）表达值变化大的基因代表细胞转录谱

• Seurat包负责筛选高变基因的函数是FindVariableFeatures()，它并不删除scRNA对象中的非高变基因。

scRNA <- FindVariableFeatures(scRNA, selection.method = "vst", nfeatures = 2000) 

• 找出的结果可以通过VariableFeatures()函数获取

top10 <- head(VariableFeatures(scRNA), 10) 
top10
# [1] "GNLY"   "IGLC2"  "S100A9" "IGLC3"  "FCGR3A" "S100A8" "CDKN1C"
# [8] "GZMB"   "ITM2C"  "LYZ" 

• 可视化

plot1 <- VariableFeaturePlot(scRNA) 
LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) 
#如下图，横坐标是某基因在所有细胞中的平均表达值，纵坐标是此基因的方差;红点即为高变基因（2000个）

1-1

补充：细胞周期相关基因

• 上一步找到的高变基因，可能会包含一些细胞周期相关基因；
• 它们会导致细胞聚类发生一定的偏移，即相同类型的细胞在聚类时会因为细胞周期的不同而分开。
• 因此有必要查看是否有细胞周期相关基因的存在；若有，则剔除

#细胞周期有关基因
head(c(cc.genes$s.genes,cc.genes$g2m.genes))
# [1] "MCM5" "PCNA" "TYMS" "FEN1" "MCM2" "MCM4"

#查看我们选择的高变基因中有哪些细胞周期相关基因
CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(scRNA))

• 在scRNA@meta.data中添加S.Score、G2M.Score和Phase三列有关细胞周期的信息。

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)
head(scRNA@meta.data)

1-2

• 观察细胞周期相关基因是否影响聚类

scRNAa <- RunPCA(scRNA, features = c(s_genes, g2m_genes))
DimPlot(scRNA, reduction = "pca", group.by = "Phase")
#如下图结果，细胞周期基因对细胞聚类的影响不大，不需要去除。

# 如需去除，代码如下
# scRNAb <- ScaleData(scRNA, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(scRNA))

1-3

2、PCA降维（线性降维）

• 使用主成分分析将2000维的信息投射到50个维度，并提取前10-20个维度（人工选择）的信息代表细胞的转录特征

scRNA <- RunPCA(scRNA, features = VariableFeatures(scRNA)) 
plot1 <- DimPlot(scRNA, reduction = "pca", group.by="orig.ident")
#（左图）根据主成分1和2的值将细胞在平面上展示出来
plot2 <- ElbowPlot(scRNA, ndims=20, reduction="pca") 
#（右图）展示前20个主成分的解释量
plot1+plot2

• 重点关注下右图，后续分析要根据右图选择提取的pc轴数量，一般选择斜率平滑的点之前的所有pc轴；此图，作者的建议是选择前18个pc轴。

  
  2-1

3、非线性降维（tSNE或UMAP）

• 最后使用非线性降维方法（tSNE或UMAP）将这10-20个PC值降维到二维空间。
• 经过上述操作，转录谱的特征信息会损失一些，但是大部分转录特征会在二维空间呈现出来。

先聚类cluster

pc.num=1:18
scRNA <- FindNeighbors(scRNA, dims = pc.num) 
# dims参数，需要指定哪些pc轴用于分析；这里利用上面的分析，选择18
scRNA <- FindClusters(scRNA, resolution = 0.5)
# resolution参数，需要指定0.1-0.9之间的一个数值，用于决定clusters的相对数量；
#数值越大，cluters越多。
table(scRNA@meta.data$seurat_clusters) #分成了0-9，共10个cluster
#  0   1   2   3   4   5   6   7   8   9 
#296 184 120 119  70  55  53  50  44  22 

tSNE

scRNA = RunTSNE(scRNA, dims = pc.num)
embed_tsne <- Embeddings(scRNA, 'tsne')
plot1 = DimPlot(scRNA, reduction = "tsne") 

3-1

UMAP

scRNA <- RunUMAP(scRNA, dims = pc.num)
embed_umap <- Embeddings(scRNA, 'umap')
plot2 = DimPlot(scRNA, reduction = "umap")

3-2

plotc <- plot1+plot2+ plot_layout(guides = 'collect')
saveRDS(scRNA, file="scRNA.rds")