单细胞分析（一）——seurat包单个样本处理

大家好，欢迎来到IT知识分享网。

10X genomics的基本原理

Read10X() 函数是针对于整理好的10X Genomics 数据，如果手头的不是类似文件，可以将其进行转换，成为格式一致的文件。

接下来使用计数矩阵创建一个 Seurat 对象。该对象用作包含单细胞数据集的数据（如计数矩阵）和分析（如 PCA 或聚类结果）的容器。例如，count matrix存储在 pbmc[[“RNA”]]@counts 中。

library(dplyr) library(Seurat) library(patchwork) 

创建对象

加载数据

# Load the PBMC dataset scdata <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/") 

创建 Seurat 对象

 2.创建Seurat对象  counts 输入的是数据，行是基因，列是细胞  project 参数输入的是项目名称,出现在metadata的orig.ident这一列  min.cells 限定的是基因：每个基因在至少多少个细胞中出现  min.features 限定的是细胞: 每个细胞中最少有多少个基因 scobj <- CreateSeuratObject(counts = scdata, project = "pbmc3k", min.cells = 3, min.features = 200) 

count matrix是什么样子?

count矩阵是稀松矩阵，可以减少占用空间

pbmc.data[c("IGF2BP2", "TCL1A", "MS4A1"), 1:30] dense.size <- object.size(as.matrix(pbmc.data)) dense.size sparse.size <- object.size(pbmc.data) sparse.size dense.size/sparse.size 

预处理流程

计算线粒体含量

这是质控的重要步骤，使用PercentageFeatureSet函数

 主要PercentageFeatureSet函数计算线粒体含量  人类使用pattern = "^MT-"，小鼠使用pattern = "^mt-" scobj[["percent.mt"]] <- PercentageFeatureSet(scobj, pattern = "^MT-")  该操作会在metadata数据里面增加一列叫做percent.mt metadata <- scobj@meta.data 

一般情况下，可以认为线粒体含量多，意味着细胞可能趋于死亡，这样的细胞就应该剔除。但是如果本身研究的就是和线粒体相关的内容，例如药物干预会引起线粒体的变化，那么就要根据具体情况来分析是否需要剔除。

质控展示

 质控数据可视化，使用VlnPlot函数  nFeature_RNA, number of Feature, 每个细胞中有多少个基因  nCount_RNA, number of counts, 每个细胞中有多少个counts  percent.mt, 我们自己增加的列, 每个细胞中线粒体基因的比例 VlnPlot(scobj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3) 

大部分情况下，线粒体的含量应该去除较大的，这个可以根据具体的情况进行分析，一般10%以内可以接受。但是依然是要根据相应的干预和特点进行取舍。

也可以展示特征之间的关系

# FeatureScatter is typically used to visualize feature-feature relationships, but can be used # for anything calculated by the object, i.e. columns in object metadata, PC scores etc. plot1 <- FeatureScatter(scobj, feature1 = "nCount_RNA", feature2 = "percent.mt") plot2 <- FeatureScatter(scobj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") plot1 + plot2 

正式筛选

这一步依然是根据metadata的数据进行筛选，这一部分依然是根据相关的具体情况进行筛选。

 正式筛选，筛选的是细胞，最终细胞减少  nFeature_RNA > 200  nFeature_RNA < 2500  percent.mt < 5 scobj <- subset(scobj, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5) 

筛选之后，再次用小提琴图展示，是上面小提琴图的缩减版本。

至于到底应该确定到多少的筛选指标，可以先按照常规指标进行设定，在后续的细胞聚类中，看看是否有与线粒体相关的群聚类出来，然后再更改筛选指标，聚类指标确实消失，那么就可以对筛选指标确定。

数据标准化

每个细胞分别进行检测，所以如果要比较各个细胞，还是要进行标准化

 先除以总数，再乘以系数，然后取log scobj <- NormalizeData(scobj, normalization.method = "LogNormalize", scale.factor = 10000)  默认参数可以省略 scobj <- NormalizeData(scobj) 

特征筛选

如果某些基因的表达再各个细胞之间表达是恒定的，所以要区分各个细胞，就要使用变化差异大的基因。

寻找高变基因

scobj <- FindVariableFeatures(scobj, selection.method = "vst", nfeatures = 2000)  默认参数可以省略 scobj <- FindVariableFeatures(scobj) 

高变基因绘图

 使用VariableFeatures 函数提取高变基因  等同于 scobj@assays$s top10 <- head(VariableFeatures(scobj), 10)  使用VariableFeaturePlot 画图 plot1 <- VariableFeaturePlot(scobj) plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE) plot1 + plot2 

缩放数据

缩放的数据达到的效果是：基因的平均值为0，方差为1。便于不同类型之间的比较。

 降维之前的必备操作 scobj <- ScaleData(scobj, features = rownames(scobj))  如果不限定参数，只会缩放高变基因  scobj <- ScaleData(scobj)  缩放后的数据存放在scobj@assays$,会很大 scale.data <- scobj@assays$RNA@scale.data 

这一步缩放数据会花费比较长的时间，原因是这一步的操作是在所有基因上展开。但是，其实我们最需要的是变异大的基因，也就是可能是标志基因的数据。所以，这个函数默认是可以使用前2000的基因进行缩放的，可以省略features参数的内容。

scobj <- ScaleData(scobj) 

可以这样操作

scobj@assays$RNA@scale.data <-matrix() 

PCA线性降维

PCA操作的对象为缩放过的数据，因此，需要保证前面的缩放数据。结果会保存在对象的reduction中，可以从中提取数据。

为什么要进行主成分分析，主要还是用少量的数据和维度，尽可能的保留原来数据特点。

降维

scobj <- RunPCA(scobj, features = VariableFeatures(object = scobj)) DimPlot(scobj, reduction = "pca")  PCA降维数据存放在scobj@reductions$pca中 pcadata = as.data.frame(scobj@reductions$pca@cell.embeddings) ggplot(pcadata,aes(PC_1,PC_2,color="red"))+ geom_point()  选择合适的PCA维度 ElbowPlot(scobj) 

降维展示

# Examine and visualize PCA results a few different ways print(scobj[["pca"]], dims = 1:5, nfeatures = 5) 

VizDimLoadings(scobj, dims = 1:2, reduction = "pca") 

使用热图进行展示

DimHeatmap(scobj, dims = 1, cells = 500, balanced = TRUE) 

dims 是确定需要展示维度数，可以尝试

UMAP非线性降维

 依赖PCA的结果  dims = 1:10 由上一步确定 scobj <- RunUMAP(scobj, dims = 1:10) DimPlot(scobj, reduction = "umap") 

 UMAP降维数据存放在scobj@reductions$umap中 umapdata = as.data.frame(scobj@reductions$umap@cell.embeddings) ggplot(umapdata,aes(UMAP_1,UMAP_2,color="red"))+ geom_point() 

聚类分群

二维空间中看到了细胞分群，那么在这个基础上，怎么进行亚群分类，这就是这一步的意义。

 找紧邻,dims = 1:10 跟UMAP相同 scobj <- FindNeighbors(scobj, dims = 1:10)  分群 scobj <- FindClusters(scobj, resolution = 0.5)  会在metadata中增加两列数据"RNA_snn_res.0.5" "seurat_clusters" metadata <- scobj@meta.data 

resolution是这个主要参数，也就是亚群切割的分辨率，这个具体怎么挑选，可以按照以下的方法可视化查看。

 设置多个resolution选择合适的resolution scobj <- FindClusters(scobj, resolution = seq(0.2,1.2,0.1)) metadata <- scobj@meta.data library(clustree) clustree(scobj) 

根据上面的可视化结果，可以大致确定一个分辨率

 选择特定分辨率得到的分群此处为RNA_snn_res.0.5 scobj@meta.data$seurat_clusters <- scobj@meta.data$RNA_snn_res.0.5 Idents(scobj) <- "seurat_clusters" DimPlot(scobj, reduction = "umap", label = T) 

大致得到了相应的分群，后续就需要对这些分群进行注释。注释是一个与个人知识背景很相关的分析过程，后续补充。

参考文章

Seurat – Guided Clustering Tutorial

Visualization of gene expression with Nebulosa (in Seurat)

Use regularized negative binomial regression to normalize UMI count data

Tutorial: Integrating stimulated vs. control PBMC datasets to learn cell-type specific responses

Understanding UMAP