胃癌单细胞数据集GSE163558复现(二)：Seurat V5标准流程

生信技能树jimmy

发布于 2024-06-13 18:24:19

1820

发布于 2024-06-13 18:24:19

文章被收录于专栏：单细胞天地单细胞天地

分享是一种态度

前言

Hello小伙伴们大家好，我是生信技能树的小学徒”我才不吃蛋黄“。今天是胃癌单细胞数据集GSE163558复现系列第二期。第一期我们进行了数据的下载与读取并成功构建Seurat对象。本期，我们将在第一期基础上走Seurat V5标准流程，对Seurat对象进行QC质控、并利用harmony整合去批次、并按标准流程进行降维聚类分群。

1质控

质控（quality control, QC）的目的是为了去除质量较差细胞，低质量的细胞会形成独特的亚群，使分群结果变得复杂；在主成分分析过程中，前几个主要成分将捕获质量差异而不是生物学差异，从而降低降维效果。因此，低质量的细胞可能会导致下游分析出现误导性结果。为了避免上述情况的发生，我们需要在下游分析开始前排除掉这些低质量细胞。

QC主要的指标有nCount_RNA(每个细胞的UMI数目)和nFeature_RNA（每个细胞中检测到的基因数量）以及"percent_mito"（表示细胞中线粒体基因的比例）这三个指标。此外，还可以纳入”percent_ribo“（核糖体基因比例）和”percent_hb“（红血细胞基因比例）两个指标。

细胞的UMI分子数和基因数反映细胞的质量，数量太低可能是细胞碎片，太高则可能是两个或多个细胞结团的情况；线粒体基因高表达的细胞，可能是处于凋亡状态或者裂解状态的细胞；核糖体基因高表达的细胞，细胞内出现RNA降解时；血红蛋白基因高表达的细胞通常为红细胞，而红细胞本身是不含有细胞核的，且携带的基因少，其携带的基因与疾病以及生物发育等过程没有太大的联系，所以可以直接剔除掉。

在这里，我们要利用PercentageFeatureSet函数分别计算每个细胞的"percent_mito"，”percent_ribo“和”percent_hb“。首先加载R包，导入第一期生成的Seurat数据：

getwd()
setwd("")
rm(list=ls())
options(stringsAsFactors = F) 
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(data.table)
library(dplyr)
sce.all <- readRDS("GSE163558.rds")

然后开始计算线粒体基因比例：

#计算线粒体基因比例
mito_genes=rownames(sce.all)[grep("^MT-", rownames(sce.all),ignore.case = T)] 
print(mito_genes) #可能是13个线粒体基因，小鼠数据基因名为小写"^mt-"
#sce.all=PercentageFeatureSet(sce.all, "^MT-", col.name = "percent_mito")
sce.all=PercentageFeatureSet(sce.all, features = mito_genes, col.name = "percent_mito")
fivenum(sce.all@meta.data$percent_mito)

计算核糖体基因比例

#计算核糖体基因比例
ribo_genes=rownames(sce.all)[grep("^Rp[sl]", rownames(sce.all),ignore.case = T)]
print(ribo_genes)
sce.all=PercentageFeatureSet(sce.all,  features = ribo_genes, col.name = "percent_ribo")
fivenum(sce.all@meta.data$percent_ribo)

再计算红血细胞基因比例

#计算红血细胞基因比例
Hb_genes=rownames(sce.all)[grep("^Hb[^(p)]", rownames(sce.all),ignore.case = T)]
print(Hb_genes)
sce.all=PercentageFeatureSet(sce.all,  features = Hb_genes,col.name = "percent_hb")
fivenum(sce.all@meta.data$percent_hb)
head(sce.all@meta.data)

可视化细胞的上述比例情况

feats <- c("nFeature_RNA", "nCount_RNA", "percent_mito",
           "percent_ribo", "percent_hb")
feats <- c("nFeature_RNA", "nCount_RNA")
p1=VlnPlot(sce.all, group.by = "orig.ident", features = feats, pt.size = 0, ncol = 2) + 
  NoLegend()
p1 
w=length(unique(sce.all$orig.ident))/3+5;w
ggsave(filename="Vlnplot1.pdf",plot=p1,width = w,height = 5)
feats <- c("percent_mito", "percent_ribo", "percent_hb")
p2=VlnPlot(sce.all, group.by = "orig.ident", features = feats, pt.size = 0, ncol = 3, same.y.lims=T) + 
  scale_y_continuous(breaks=seq(0, 100, 5)) +
  NoLegend()
p2 
w=length(unique(sce.all$orig.ident))/2+5;w
ggsave(filename="Vlnplot2.pdf",plot=p2,width = w,height = 5)

p3=FeatureScatter(sce.all, "nCount_RNA", "nFeature_RNA", group.by = "orig.ident", pt.size = 0.5)
p3
ggsave(filename="Scatterplot.pdf",plot=p3)

结果如下：

p1：nCount_RNA与nFeature_RNA

p2："percent_mito"，”percent_ribo“和”percent_hb“

p3：nCount_RNA与nFeature_RNA相关性分析

根据上述指标，过滤低质量细胞/基因

过滤指标1:最少表达基因数的细胞 and 最少表达细胞数的基因

一般来说，在CreateSeuratObject的时候已经是进行了这个过滤操作，如果后期看到了自己的单细胞降维聚类分群结果很诡异，就可以回过头来看质量控制环节，先走默认流程即可。

 if(F){
  selected_c <- WhichCells(sce.all, expression = nFeature_RNA > 500)
  selected_f <- rownames(sce.all)[Matrix::rowSums(sce.all@assays$RNA$counts > 0 ) > 3]
  sce.all.filt <- subset(sce.all, features = selected_f, cells = selected_c)
  dim(sce.all) 
  dim(sce.all.filt) 
}
sce.all.filt =  sce.all
# par(mar = c(4, 8, 2, 1))
# 这里的C 这个矩阵，有一点大，可以考虑随抽样 
C=subset(sce.all.filt,downsample=100)@assays$RNA$counts
dim(C)
C=Matrix::t(Matrix::t(C)/Matrix::colSums(C)) * 100

most_expressed <- order(apply(C, 1, median), decreasing = T)[50:1]

pdf("TOP50_most_expressed_gene.pdf",width=14)
boxplot(as.matrix(Matrix::t(C[most_expressed, ])),
        cex = 0.1, las = 1, 
        xlab = "% total count per cell", 
        col = (scales::hue_pal())(50)[50:1], 
        horizontal = TRUE)
dev.off()
rm(C)

过滤指标2:线粒体/核糖体基因比例(根据上面的violin图)

这里，我们筛选的标准是：percent_mito < 25，percent_ribo > 3，percent_hb < 1。不同的数据集，不同的组织，需要根据特定情况来调整筛选阈值。

selected_mito <- WhichCells(sce.all.filt, expression = percent_mito < 25)
selected_ribo <- WhichCells(sce.all.filt, expression = percent_ribo > 3)
selected_hb <- WhichCells(sce.all.filt, expression = percent_hb < 1 )
length(selected_hb)
length(selected_ribo)
length(selected_mito)
sce.all.filt <- subset(sce.all.filt, cells = selected_mito)
sce.all.filt <- subset(sce.all.filt, cells = selected_ribo)
sce.all.filt <- subset(sce.all.filt, cells = selected_hb)
dim(sce.all.filt)
table(sce.all.filt$orig.ident)
length(sce.all.filt$orig.ident)

我们再统计一下QC前后各样本细胞数：

第一期QC前各样本细胞数：

本期QC后各样本细胞数：

可以看到，QC后细胞数量由53093个变成45548个，过滤成功（特别是GSM5004186_O1这个样本，很多细胞线粒体基因的含量较高，未过滤前有5837个细胞，过滤掉低质量的细胞之后，剩下1909个）。

可视化过滤后的情况：

feats <- c("nFeature_RNA", "nCount_RNA")
p1_filtered=VlnPlot(sce.all.filt, group.by = "orig.ident", features = feats, pt.size = 0, ncol = 2) + 
  NoLegend()
w=length(unique(sce.all.filt$orig.ident))/3+5;w 
ggsave(filename="Vlnplot1_filtered.pdf",plot=p1_filtered,width = w,height = 5)

feats <- c("percent_mito", "percent_ribo", "percent_hb")
p2_filtered=VlnPlot(sce.all.filt, group.by = "orig.ident", features = feats, pt.size = 0, ncol = 3) + 
  NoLegend()
w=length(unique(sce.all.filt$orig.ident))/2+5;w 
ggsave(filename="Vlnplot2_filtered.pdf",plot=p2_filtered,width = w,height = 5)

p1_filtered

p2_filtered

2 harmony整合多个单细胞样品

细胞筛选之后，我们需要进行harmony整合。第一期我们在创建总的Seurat对象时，使用了merge函数对多个Seurat进行了简单的合并。merge只是按照行和列进行了合并，并未对数据进行其他处理。

在拿到下游单细胞矩阵前，样本经历了多个实验环节的处理，时间、处理人员、试剂以及技术平台等变量都会成为混杂因素。以上因素混合到一起，就会导致数据产生批次效应（batch effect）。为了尽可能避免混杂因素，我们可以严格把控测序的技术流程，同时也需要在下游分析中进行事后补救（算法去批次）。目前单细胞测序常用的去批次算法有Harmony，CCA，RPCA,FastMNN,scVI等。在这里，我们采用Harmony进行演示。

在Harmony去批次前，我们需要进行以下流程：

数据标准化

数据标准化的方法是通过对原始表达值进行对数转换"LogNormalize"，使其总体更加符合正态分布。经过对数转换之后，不同基因或细胞之间也更具有可比性，从一定程度上消除不同细胞之间的技术差异。

sce.all.filt <- NormalizeData(sce.all.filt, 
                             normalization.method = "LogNormalize",
                             scale.factor = 1e4)

筛选高变基因

高变异基因就是highly variable features（HVGs）是指在某些细胞中高度表达，而在其他细胞中低度表达的基因。默认情况下，此步骤回筛选2000个HVGs用于下游分析。

sce.all.filt <- FindVariableFeatures(sce.all.filt)
p4 <- VariableFeaturePlot(sce.all.filt) 
p4

高变基因

数据归一化

scale归一化：将每个基因在所有细胞中的均值变为0，方差标为1，赋予每个基因在下游分析中同样的权重，不至于使高表达基因占据主导地位

sce.all.filt <- ScaleData(sce.all.filt)

PCA线性降维

单细胞测序是一种高通量测序技术，它产生的数据集在细胞和基因数量上都具有很高的维度。这立即指向了一个事实，即单细胞测序数据受到了"维度诅咒"的困扰（单细胞最好的教程（四）：降维）。

"维度诅咒"这个概念最早是由R. Bellman提出的，它描述的是理论上高维数据包含更多的信息，但在实践中并非如此。更高维度的数据往往包含更多的噪声和冗余，因此增加更多的信息并不利于后续的分析步骤。

为了在数据的处理和可视化更加便捷的同时，保留数据重要的信息，在这里我们需要应用PCA（principal components analysis）即主成分分析技术，降低数据维度（单细胞PCA降维结果理解）。

sce.all.filt <- RunPCA(sce.all.filt, features = VariableFeatures(object = sce.all.filt))
##可视化PCA结果
VizDimLoadings(sce.all.filt, dims = 1:2, reduction = "pca")
DimPlot(sce.all.filt, reduction = "pca") + NoLegend()
DimHeatmap(sce.all.filt, dims = 1:12, cells = 500, balanced = TRUE)

PCA可视化结果

Harmony去批次

seuratObj <- RunHarmony(sce.all.filt, "orig.ident")
names(seuratObj@reductions)
[1] "pca"     "harmony"

我们可以看到，RunHarmony后，Harmony的结果已经在seuratObj的reductions中。

然后使用UMAP/TSNE可视化Harmony去批次效果：

seuratObj <- RunUMAP(seuratObj,  dims = 1:15, 
                       reduction = "harmony")
DimPlot(seuratObj,reduction = "umap",label=F ) 
seuratObj <- RunTSNE(seuratObj, dims = 1:15, 
                       reduction = "harmony")

DimPlot(seuratObj,reduction = "tsne",label=F )

UMAP1

TSNE1

我们再来看一下不用Harmony去批次的效果：

UMAP2

TSNE2

两组图对比我们可以看到，运行Harmony后，各样本散点融合的更好，表明该数据集需要去批次且Harmony去批次效果较好。

3 细胞聚类

Seurat软件使用基于图论的聚类算法对细胞进行聚类和分群，要包括以下步骤：

1.构建细胞间的聚类关系：利用PCA空间中的欧几里得距离构建KNN聚类关系图；
2.优化细胞间聚类关系距离的权重值：利用Jaccard相似性优化任意两个细胞间的边缘权重；
3.聚类和分群：使用Louvain 算法进行细胞群聚类优化。

分群标准的确定使用了两个函数FindNeighbors和FindClusters。

FindNeighbors函数用于计算给定数据集的最近邻图，可以返回包含KNN和SNN的对象列表。还可以通过迭代分组来优化聚类结果，以最大化标准模块度函数。

FindClusters函数是基于共享最近邻（SNN）模块化优化的聚类算法识别细胞簇。该函数的重要参数为分辨率resolution,resolution最小值为0，分为1类；值越大，分群数越多；在0.4-1.2之间通常会对3k 左右的单细胞数据集产生良好的结果。

我们先运行FindNeighbors：

sce.all.filt=seuratObj

sce.all.filt <- FindNeighbors(sce.all.filt, reduction = "harmony",
                             dims = 1:15) 

sce.all.filt.all=sce.all.filt

然后再运行FindClusters，在这里我们使用了for循环，设置了不同的分辨率，观察分群效果。

#设置不同的分辨率，观察分群效果(选择哪一个？)
  for (res in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1)) {
    sce.all.filt.all=FindClusters(sce.all.filt.all, #graph.name = "CCA_snn", 
                               resolution = res, algorithm = 1)
  }
  colnames(sce.all.filt.all@meta.data)
  apply(sce.all.filt.all@meta.data[,grep("RNA_snn",colnames(sce.all.filt.all@meta.data))],2,table)
  
  p1_dim=plot_grid(ncol = 3, DimPlot(sce.all.filt.all, reduction = "umap", group.by = "RNA_snn_res.0.01") + 
                     ggtitle("louvain_0.01"), DimPlot(sce.all.filt.all, reduction = "umap", group.by = "RNA_snn_res.0.1") + 
                     ggtitle("louvain_0.1"), DimPlot(sce.all.filt.all, reduction = "umap", group.by = "RNA_snn_res.0.2") + 
                     ggtitle("louvain_0.2"))
  ggsave(plot=p1_dim, filename="Dimplot_diff_resolution_low.pdf",width = 14)
  
  p1_dim=plot_grid(ncol = 3, DimPlot(sce.all.filt.all, reduction = "umap", group.by = "RNA_snn_res.0.8") + 
                     ggtitle("louvain_0.8"), DimPlot(sce.all.filt.all, reduction = "umap", group.by = "RNA_snn_res.1") + 
                     ggtitle("louvain_1"), DimPlot(sce.all.filt.all, reduction = "umap", group.by = "RNA_snn_res.0.3") + 
                     ggtitle("louvain_0.3"))
  ggsave(plot=p1_dim, filename="Dimplot_diff_resolution_high.pdf",width = 18)
  
  p2_tree=clustree(sce.all.filt.all@meta.data, prefix = "RNA_snn_res.")
  ggsave(plot=p2_tree, filename="Tree_diff_resolution.pdf")
  table(sce.all.filt.all@active.ident)