Seurat

\(10X\) 数据读取、标准化数据、特征选择、主成分分析、细胞聚类、非线性降维、簇生物标志物选取、集群注释

整体代码

setwd("E:\\AAA大创\\代码注释版\\暑假练习\\Seurat")
# BiocManager::install("Seurat")
library(Seurat)
library(ggplot2)
library(dplyr)
library(patchwork)

#数据导入
# 导入数据
pbmc.data<-Read10X(data.dir = "./filtered_gene_bc_matrices/hg19")
# 初始化seurat数据
pbmc<-CreateSeuratObject(counts = pbmc.data,project = "pbmc3k",min.cells = 3,min.features = 200)
# pbmc
# 基因过滤：去除在少于min.cells个细胞中表达的基因。
# 细胞过滤：去除检测到基因数（即特征数）少于min.features的细胞。

#数据预处理
#使用PercentageFeatureSet函数计算线粒体QC指标
pbmc[['percent.mt']]<-PercentageFeatureSet(pbmc,pattern = "^MT-")
#nFeature_RNA：代表每个细胞中检测到的基因数量。
#nCount_RNA：代表每个细胞中 RNA 分子的总数。高值可能表示高表达的细胞，
#但过高也可能是细胞双重性（doublets）或低质量细胞的标志。
#percent.mt：代表每个细胞中线粒体基因的表达百分比。
#一般认为线粒体含量增高的细胞为裂解的细胞（死细胞），
#而活细胞中检测到的线粒体含量应低于10%。但如果是高代谢的细胞就需要谨慎过滤该指标。
#percent.rp：代表每个细胞中核糖体蛋白基因的表达百分比。
#percent.hb：代表血红蛋白基因的表达百分比。一般不研究这个可以过滤
# 使用violin plot可视化  QC指标，并使用这些指标过滤单元格
VlnPlot(pbmc,features = c("nFeature_RNA","nCount_RNA","percent.mt"),ncol = 3)
# FeatureScatter 通常用于可视化两个特征之间的关系
plot1<-FeatureScatter(pbmc,feature1 = "nCount_RNA",feature2 = "percent.mt")
plot2<-FeatureScatter(pbmc,feature1 = "nCount_RNA",feature2 = "nFeature_RNA")
plot1+plot2
# 将QC指标可视化，并使用这些指标过滤单元格
# 过滤具有2500或少于200的独特特征计数的单元格，过滤线粒体计数>5%的细胞
pbmc<-subset(pbmc,subset=nFeature_RNA>200 & nFeature_RNA<2500 & percent.mt<5)
# 数据规范化
#采用全局缩放归一化方法（LogNormalize）,该方法将每个单元格的特征表达式
#测量值按总表达进行归一化，将其乘以比例因子（默认为10000），并对结果进行对数转换。
pbmc<-NormalizeData(pbmc,normalization.method = "LogNormalize",scale.factor = 10000)

#识别高度可变基因
# 识别高度可变的特征（特征选择）
#vst表示方差稳定变换法
#保留前2000个变异度最高的基因
pbmc<-FindVariableFeatures(pbmc,selection.method = "vst",nfeatures = 2000)
# 确定高表达的前十个基因
top10<-head(VariableFeatures(pbmc),10)
plot1<-VariableFeaturePlot(pbmc)
plot2<-LabelPoints(plot = plot1,points = top10,repel = TRUE)
#标记前十个基因标签，repel避免标签和点重叠
pdf("Var.pdf", width = 15, height = 6, family = 'GB1')
plot1+plot2
#pdf无法输出是因为dev.off()没有关闭pdf输出，多关几次就好
#数据缩放
#如PCA的降维计数之前的标准预处理步骤，使得每个基因平均表达为0， 方差为1
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

#PCA
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
print(pbmc[['pca']], dims = 1:5, nfeatures = 5)
#可视化
# Seurat提供可视化细胞和定义PCA,包括功能的几种有用的方法
#ViziDimReduction(),DimPlot()和DimHeatmap()
VizDimLoadings(pbmc,dims = 1:2,reduction = "pca")
#使用pca降维结果
DimPlot(pbmc,reduction = "pca")
#默认显示前两个pca
DimHeatmap(pbmc,dims = 1:2,cells = 500,balanced = TRUE)
#仅显示PC1，使用500个细胞，平衡选择高低得分细胞
# 确定数据集的维度
pbmc<-JackStraw(pbmc,num.replicate = 100)
#数字为执行置换检验的次数，数值越大越精确越慢，常用100
pbmc<-ScoreJackStraw(pbmc,dims = 1:20)
#评估前20个PC的显著性
JackStrawPlot(pbmc, dims = 1 : 15)
#可视化，观察PC在哪里急剧下降
ElbowPlot(pbmc)
#生成“肘图”通过观察肘部位置来寻找显著的PC
# 聚类细胞
# 建立KNN图，并基于其局部领域中的共享重叠细化任意两个单元之间的边权重
#使用前十个PC
pbmc<-FindNeighbors(pbmc,dims = 1:10)
#对细胞进行聚类，应用模块化优化技术，Louvain算法（默认）或SLM
# 以迭代方式将细胞分组在一起，目标是优化标准模块化函数
#参数控制聚类粒度，常用0.4-1.2
pbmc<-FindClusters(pbmc,resolution = 0.5)
head(Idents(pbmc),5)

#非线性聚类 UMAP/tSNE
pbmc <- RunUMAP(pbmc, dims = 1:10)
#运行非线性降维
DimPlot(pbmc, reduction = 'umap')
# 寻找差异表达的特征（簇生物标志物）
# findmarkers为所有集群自动执行此过程，也可以测试集群组之间的对比，或针对所有单元格进行测试
# 默认情况下，ident.1与所有其他细胞相比，他识别单个簇的阳性和阴性标记。
# min.pct参数要求在两组细胞中的任何一组中以最小百分比检测到一个特征
# 而 thresh.test 参数要求一个特征在两组之间差异表达（平均）一定量
# 寻找cluster2的所有markers
cluster2.markers<-FindMarkers(pbmc,ident.1 = 2,min.pct = 0.25)
#ident.1 = 2：目标簇（簇2）
#min.pct = 0.25：基因必须在簇2中≥25%的细胞表达，或在对照细胞中≥25%表达
head(cluster2.markers,n=5)

# 寻找cluster5中与cluster0和cluster3n不同的所有markers
cluster5.markers<-FindMarkers(pbmc,ident.1 = 5,ident.2=c(0,3),min.pct = 0.25)
#找出簇5相对于簇0和簇3组合的特异性基因
head(cluster5.markers,n=5)

# 找出每个细胞簇的标记物，与所有剩余的细胞进行比较，只报告阳性细胞 
pbmc.markers<-FindAllMarkers(pbmc,only.pos = TRUE,min.pct = 0.25,logfc.threshold =0.25 )
#只保留上调基因， 表达比例阈值， 最小logFC
pbmc.markers %>%
  group_by(cluster) %>%
  slice_max(n=2,order_by = avg_log2FC)
#按簇分组，按表达倍数排序，每簇取top2基因
#小提琴图可视化
# 还有用于可视化标记表达的工具，Vlnplot(显示跨集群的表达概率分布)和FeaturePlot()(在tSNE或PCA图上可视化特征表达)是最常用的可视化，建议探索RidgePlot(),CellScatter（），和DotPlot()作为查看数据集的其他方法
VlnPlot(pbmc,features = c("MS4A1","CD79A"))
# 使用原始未归一化计数数据
VlnPlot(pbmc,features = c("NKG7","PF4"),slot = "counts",log = TRUE)
#观察个别基因分布
FeaturePlot(pbmc,features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",
                              "CD8A"))
#DoHeatmap()为给定的细胞和特征生成一个表达热图。在这种情况下，
#我们为每个集群绘制前 10 个标记（或所有标记，如果小于 10）。
pbmc.markers%>%
  group_by(cluster)%>%
  top_n(n=10,wt=avg_log2FC)->top10
DoHeatmap(pbmc,features = top10$gene)+NoLegend()
#NoLegend去除右侧标签

#注释分群
# 将细胞类型标识分配给集群
new.cluster.ids<-c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono",
                   "NK", "DC", "Platelet")
names(new.cluster.ids)<-levels(pbmc)
#levels(pbmc) 返回当前 pbmc 对象中 cluster 的名称（比如 "0", "1", ..., "8"）。
#此步骤建立了旧编号和新名称之间的映射关系。
pbmc<-RenameIdents(pbmc,new.cluster.ids)
DimPlot(pbmc,reduction = "umap",label = TRUE,pt.size = 0.5)+NoLegend()
#使用 UMAP 降维结果作图，展示细胞的空间分布。
#label = TRUE：在图中直接标注每个 cluster 的细胞类型名称。
#pt.size = 0.5：每个点（细胞）大小设置为 0.5。
#+ NoLegend()：去除右侧图例（因为标签已经直接标注在图上）。
saveRDS(pbmc,file = "./pbmc3k_final.rds")
#将处理好的 pbmc Seurat 对象保存为 RDS 文件，以便后续载入使用，不需要重复处理数据。
#文件名为 "pbmc3k_final.rds"，保存在当前目录的子目录 ./ 下

一、\(10X\) 数据读取

分析数据来源

Show Code

#数据导入
# 导入数据
pbmc.data<-Read10X(data.dir = "./filtered_gene_bc_matrices/hg19")
# 初始化seurat数据
pbmc<-CreateSeuratObject(counts = pbmc.data,project = "pbmc3k",min.cells = 3,min.features = 200)
# pbmc
# 基因过滤：去除在少于min.cells个细胞中表达的基因。
# 细胞过滤：去除检测到基因数（即特征数）少于min.features的细胞。

二、标准化数据

Show Code

#数据预处理
#使用PercentageFeatureSet函数计算线粒体QC指标
pbmc[['percent.mt']]<-PercentageFeatureSet(pbmc,pattern = "^MT-")
#nFeature_RNA：代表每个细胞中检测到的基因数量。
#nCount_RNA：代表每个细胞中 RNA 分子的总数。高值可能表示高表达的细胞，
#但过高也可能是细胞双重性（doublets）或低质量细胞的标志。
#percent.mt：代表每个细胞中线粒体基因的表达百分比。
#一般认为线粒体含量增高的细胞为裂解的细胞（死细胞），
#而活细胞中检测到的线粒体含量应低于10%。但如果是高代谢的细胞就需要谨慎过滤该指标。
#percent.rp：代表每个细胞中核糖体蛋白基因的表达百分比。
#percent.hb：代表血红蛋白基因的表达百分比。一般不研究这个可以过滤
# 使用violin plot可视化  QC指标，并使用这些指标过滤单元格
VlnPlot(pbmc,features = c("nFeature_RNA","nCount_RNA","percent.mt"),ncol = 3)
# FeatureScatter 通常用于可视化两个特征之间的关系
plot1<-FeatureScatter(pbmc,feature1 = "nCount_RNA",feature2 = "percent.mt")
plot2<-FeatureScatter(pbmc,feature1 = "nCount_RNA",feature2 = "nFeature_RNA")
plot1+plot2
# 将QC指标可视化，并使用这些指标过滤单元格
# 过滤具有2500或少于200的独特特征计数的单元格，过滤线粒体计数>5%的细胞
pbmc<-subset(pbmc,subset=nFeature_RNA>200 & nFeature_RNA<2500 & percent.mt<5)
# 数据规范化
#采用全局缩放归一化方法（LogNormalize）,该方法将每个单元格的特征表达式
#测量值按总表达进行归一化，将其乘以比例因子（默认为10000），并对结果进行对数转换。
pbmc<-NormalizeData(pbmc,normalization.method = "LogNormalize",scale.factor = 10000)

结果

三、特征选择

Show Code

#识别高度可变基因
# 识别高度可变的特征（特征选择）
#vst表示方差稳定变换法
#保留前2000个变异度最高的基因
pbmc<-FindVariableFeatures(pbmc,selection.method = "vst",nfeatures = 2000)
# 确定高表达的前十个基因
top10<-head(VariableFeatures(pbmc),10)
plot1<-VariableFeaturePlot(pbmc)
plot2<-LabelPoints(plot = plot1,points = top10,repel = TRUE)
#标记前十个基因标签，repel避免标签和点重叠
pdf("Var.pdf", width = 15, height = 6, family = 'GB1')
plot1+plot2
#pdf无法输出是因为dev.off()没有关闭pdf输出，多关几次就好
#数据缩放
#如PCA的降维计数之前的标准预处理步骤，使得每个基因平均表达为0， 方差为1
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

结果

四、主成分分析（PCA）

Show Code

#PCA
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
print(pbmc[['pca']], dims = 1:5, nfeatures = 5)
#可视化
# Seurat提供可视化细胞和定义PCA,包括功能的几种有用的方法
#ViziDimReduction(),DimPlot()和DimHeatmap()
VizDimLoadings(pbmc,dims = 1:2,reduction = "pca")
#使用pca降维结果
DimPlot(pbmc,reduction = "pca")
#默认显示前两个pca
DimHeatmap(pbmc,dims = 1:2,cells = 500,balanced = TRUE)
#仅显示PC1，使用500个细胞，平衡选择高低得分细胞
# 确定数据集的维度
pbmc<-JackStraw(pbmc,num.replicate = 100)
#数字为执行置换检验的次数，数值越大越精确越慢，常用100
pbmc<-ScoreJackStraw(pbmc,dims = 1:20)
#评估前20个PC的显著性
JackStrawPlot(pbmc, dims = 1 : 15)
#可视化，观察PC在哪里急剧下降
ElbowPlot(pbmc)
#生成“肘图”通过观察肘部位置来寻找显著的PC

结果

六、细胞聚类

Show Code

# 聚类细胞
# 建立KNN图，并基于其局部领域中的共享重叠细化任意两个单元之间的边权重
#使用前十个PC
pbmc<-FindNeighbors(pbmc,dims = 1:10)
#对细胞进行聚类，应用模块化优化技术，Louvain算法（默认）或SLM
# 以迭代方式将细胞分组在一起，目标是优化标准模块化函数
#参数控制聚类粒度，常用0.4-1.2
pbmc<-FindClusters(pbmc,resolution = 0.5)
head(Idents(pbmc),5)

结果

七、非线性降维

Show Code

#非线性聚类 UMAP/tSNE
pbmc <- RunUMAP(pbmc, dims = 1:10)
#运行非线性降维
DimPlot(pbmc, reduction = 'umap')

结果

八、簇生物标志物选取

Show Code

# 寻找差异表达的特征（簇生物标志物）
# findmarkers为所有集群自动执行此过程，也可以测试集群组之间的对比，或针对所有单元格进行测试
# 默认情况下，ident.1与所有其他细胞相比，他识别单个簇的阳性和阴性标记。
# min.pct参数要求在两组细胞中的任何一组中以最小百分比检测到一个特征
# 而 thresh.test 参数要求一个特征在两组之间差异表达（平均）一定量
# 寻找cluster2的所有markers
cluster2.markers<-FindMarkers(pbmc,ident.1 = 2,min.pct = 0.25)
#ident.1 = 2：目标簇（簇2）
#min.pct = 0.25：基因必须在簇2中≥25%的细胞表达，或在对照细胞中≥25%表达
head(cluster2.markers,n=5)

# 寻找cluster5中与cluster0和cluster3n不同的所有markers
cluster5.markers<-FindMarkers(pbmc,ident.1 = 5,ident.2=c(0,3),min.pct = 0.25)
#找出簇5相对于簇0和簇3组合的特异性基因
head(cluster5.markers,n=5)

# 找出每个细胞簇的标记物，与所有剩余的细胞进行比较，只报告阳性细胞 
pbmc.markers<-FindAllMarkers(pbmc,only.pos = TRUE,min.pct = 0.25,logfc.threshold =0.25 )
#只保留上调基因， 表达比例阈值， 最小logFC
pbmc.markers %>%
  group_by(cluster) %>%
  slice_max(n=2,order_by = avg_log2FC)
#按簇分组，按表达倍数排序，每簇取top2基因
#小提琴图可视化
# 还有用于可视化标记表达的工具，Vlnplot(显示跨集群的表达概率分布)和FeaturePlot()(在tSNE或PCA图上可视化特征表达)是最常用的可视化，建议探索RidgePlot(),CellScatter（），和DotPlot()作为查看数据集的其他方法
VlnPlot(pbmc,features = c("MS4A1","CD79A"))
# 使用原始未归一化计数数据
#VlnPlot(pbmc,features = c("NKG7","PF4"),slot = "counts",log = TRUE)
#观察个别基因分布
FeaturePlot(pbmc,features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",
                              "CD8A"))
#DoHeatmap()为给定的细胞和特征生成一个表达热图。在这种情况下，
#我们为每个集群绘制前 10 个标记（或所有标记，如果小于 10）。
pbmc.markers%>%
  group_by(cluster)%>%
  top_n(n=10,wt=avg_log2FC)->top10
DoHeatmap(pbmc,features = top10$gene)+NoLegend()
#NoLegend去除右侧标签

结果

九、集群注释

Show Code

#注释分群
# 将细胞类型标识分配给集群
new.cluster.ids<-c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono",
                   "NK", "DC", "Platelet")
names(new.cluster.ids)<-levels(pbmc)
#levels(pbmc) 返回当前 pbmc 对象中 cluster 的名称（比如 "0", "1", ..., "8"）。
#此步骤建立了旧编号和新名称之间的映射关系。
pbmc<-RenameIdents(pbmc,new.cluster.ids)
DimPlot(pbmc,reduction = "umap",label = TRUE,pt.size = 0.5)+NoLegend()
#使用 UMAP 降维结果作图，展示细胞的空间分布。
#label = TRUE：在图中直接标注每个 cluster 的细胞类型名称。
#pt.size = 0.5：每个点（细胞）大小设置为 0.5。
#+ NoLegend()：去除右侧图例（因为标签已经直接标注在图上）。
saveRDS(pbmc,file = "./pbmc3k_final.rds")
#将处理好的 pbmc Seurat 对象保存为 RDS 文件，以便后续载入使用，不需要重复处理数据。
#文件名为 "pbmc3k_final.rds"，保存在当前目录的子目录 ./ 下

结果