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")
