前言
本文是基于seurat v5版本的教程复现,主要依据来源于官网seurat v5建议大家有时间的话可以去官网进行详细阅读学习!
一、Seurat V4和 Seurat V5的比较优势:
更强大的预处理算法：Seurat v5引入了新的预处理算法，如CellBender和SCTransform，它们在去除噪音和批次效应方面表现更好。这些算法可以提高数据的质量和准确性。
更好的集成分析能力：Seurat v5提供了新的集成分析框架，能够更好地将多个数据集进行整合和比较。这种集成分析能够更准确地识别细胞类型和状态的差异。
多模态数据分析：Seurat v5支持多模态数据的分析，如融合RNA和蛋白质表达数据的分析。这使得研究人员能够更全面地研究单细胞数据，揭示不同数据模态之间的相互作用。
更丰富的可视化工具：Seurat v5提供了新的可视化工具，如可视化互作网络图和互动式数据探索工具。这些工具使得结果的可视化更加方便和直观。
劣势:
学习曲线较陡：Seurat v5相较于Seurat v4在功能上有所扩展，因此需要花费一些时间学习和适应新版本的功能和使用方法。
需要更高的计算资源：由于Seurat v5提供了更多的功能和分析选项，因此可能需要更高的计算资源（如计算能力和内存），以处理更大规模的数据集和更复杂的分析任务。
总之,seurat版本更新之后的效果又好有坏,大家可以自行选择使用哪个版本,不过市面上的很多资源都是基于seurat v4的,但是seurat v5也有其优势,下面则是seurat v5的标准化过程,与seurat有很多相似之处,可供借鉴和参考!
二、标准化流程 1.引入库代码如下（示例）：
# 主要用于做单细胞分析的各种函数 library(Seurat) # 主要与批次效应有关 library(harmony) # 广泛使用的数据分析和可视化的工具集 library(tidyverse) # 是一个功能强大的R绘图包 library(ggplot2) # 增强绘图和图形组合的包 library(cowplot) # 用于数据的操作和处理,管道,过滤,排序等 library(dplyr) # 用于组合多个图形 library(patchwork) # 用于文件的处理 library(R.utils) 2.设置默认工作路径代码如下（示例）：
setwd('G:/diversity_intergration') 3.批量读取单细胞数据代码如下（示例）：
# 使用list.files函数列出文件夹中的所有文件 dir_name <- list.files('data/') dir_name # 创建一个空的list对象来存储转化之后的seurat对象 scRNAlist <- list() for (i in 1:length(dir_name)){ counts <- Read10X(data.dir = paste('G:/diversity_intergration/data/',dir_name[i],sep = '')) # 使用createseuratobject来创建seurat对象,使用counts矩阵,设置样本名为目录名 scRNAlist[[i]] <- CreateSeuratObject(counts = counts, project = strsplit(dir_name[i],'_')[[1]][2], min.cells = 3, min.features = 200) } scRNAlist # 查看样本信息 4.批量计算每个样本的线粒体和红细胞的比例(质控 — 线粒体,红细胞,核糖体)代码如下（示例）：
# 注意如果是小鼠为 '^mt-' 人为 '^MT-' # 如果是小鼠的话,将49-53行改为sc[['HB_percent']] <- PercentageFeatureSet(sc,pattern = "^Hb[^(p)]") for (i in 1:length(scRNAlist)){ sc <- scRNAlist[[i]] # 获取scRNAlist中的第i个seurat对象 # 计算线粒体的比例 sc[['mt_percent']] <- PercentageFeatureSet(sc,pattern = '^MT-') # 计算红细胞的比例 HB_genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ") # 定义红细胞的基因列表 HB_m <- match(HB_genes,rownames(sc@assays$RNA)) # 在seurat对象的RNA数据中有红细胞基因的索引位置 HB_genes <- rownames(sc@assays$RNA)[HB_m] # 获取匹配到的红细胞基因的行名 HB_genes <- HB_genes[!is.na(HB_genes)] # 删除na值(未匹配到的基因) sc[['HB_percent']] <- PercentageFeatureSet(sc,features = HB_genes) # 计算红细胞基因的比例,并将其储存在新列HB_percent中 # 将sc的值赋给scRNAlist scRNAlist[[i]] <- sc # 删除sc rm(sc) } 5.批量绘制质控前的小提琴图代码如下（示例）：
for (i in 1:length(scRNAlist)){ violin_before[[i]] <- VlnPlot(scRNAlist[[i]], features = c('nFeature_RNA','nCount_RNA','mt_percent','HB_percent'), pt.size = 0.01, ncol = 4) } violin_before # 输出质控前的图片,有多少个样本就会输出多少张图片 violin_before[[2]] 6.批量过滤细胞,MT,HB基因代码如下（示例）：
scRNAlist <- lapply(scRNAlist, FUN = function(x){ x <- subset(x, # 怀疑有双细胞的nFeature_RNA可以设置低一点,依据这个样本的小提琴图最低设置4000 subset=nFeature_RNA > 200 & nFeature_RNA < 3500 & mt_percent < 7 & HB_percent < 3 & # 红细胞基本没有,随便过滤1-5%都可以 nCount_RNA < quantile(nCount_RNA,0.97) & # 过滤最高的 3% 的nCount_RNA nCount_RNA > 1000) }) # 一般默认线粒体的含量要少于20% ,红细胞的数据至少要少于5% view(scRNAlist[[1]]@meta.data) 7.merge合并样本代码如下（示例）：
# 手动设定一个样本名称 dir_name <- c('26','27','28') scRNAlist_merge <- merge(x = scRNAlist[[1]],y = scRNAlist[-1],add.cell.ids = dir_name) # 质控后的数据(小提琴图) violin_after <- VlnPlot(scRNAlist_merge, features = c('nFeature_RNA','nCount_RNA','mt_percent','HB_percent'), split.by = 'orig.ident', layer = 'counts', pt.size = 0.01, ncol = 4) # 统计细胞数量 table(scRNAlist_merge$orig.ident) # tips:seuratv5中币seuratv4中多了一个layers,并且layers层下的每个样本都是分开的,这样做的主要目的就是为了优化内存 8.meatadata增加分组信息(可选步骤)代码如下（示例）：
# 方法很多,也可以不增加分组,因为orig.ident就已经是一个分组了 # 这里简单介绍依照分组,但是这个案例中的数据不会使用,只有normal组 group <- str_split(colnames(scRNAlist_merge@assays$RNA),'-',simplify = T)[,1] # 这里的simplify表示将结果简化为一个矩阵形式,并且取出第一列 group <- ifelse(str_detect(group,"N$"),'normal','tumor') table(group) # 增加一列group scRNAlist_merge$group <- group 9.标准化&高变基因&归一化&PCA代码如下（示例）：
scRNAlist_merge <- NormalizeData(scRNAlist_merge) scRNAlist_merge <- FindVariableFeatures(scRNAlist_merge,nFeature_RNA = 2000) # 这里指定最后想要得到的基因数量 scRNAlist_merge <- ScaleData(scRNAlist_merge,vars.to.regress = c('mt_percent')) # 这里去除线粒体的影响,当然这里如果有细胞周期的影响,我们要把细胞周期也回归一下,去除影响 scRNAlist_merge <- RunPCA(scRNAlist_merge) # 可以SCT来替换此步骤(初学者不建议使用) --- 109-111 # scRNAlist_merge <- SCTransform(scRNAlist_merge,vars.to.regress = 'mt_percent') 10.选择harmony,CCA,scVI,RPCA,FastMNN等整合去批次代码如下（示例）：
# 以下为官网的五种去批次的方法,运算速率相比seuratv4快了很多,特别是CCA # Anchor-based CCA integration(method = CCAIntegration) # Anchor-based RPCA integration(method = RPCAIntegration) # harmony (method = HarmonyIntegration) # FastMNN (method = FastMNNIntegration) # scVI (method = scVIInteration) scVI得创建python环境,新手不建议使用 # Integrated with CCA scRNA_cca <- IntegrateLayers(object = scRNAlist_merge, method = CCAIntegration, # 需要什么方法改变这里即可 orig.reduction = 'pca', # 这里降维必须选择pca new.reduction = 'integrated.cca') # 储存在新的降维结果integrated.cca中 # Integrated with harmony scRNA_harmony <- IntegrateLayers(object = scRNAlist_merge, method = HarmonyIntegration, # 需要什么方法改变这里即可 orig.reduction = 'pca', # 这里降维必须选择pca new.reduction = 'harmony') # 储存在新的降维结果integrated.cca中 ############## 9.joinLayers合并样本counts和data ################# scRNA_harmony[['RNA']] <- JoinLayers(scRNA_harmony[['RNA']]) 11.降维聚类umap代码如下（示例）：
ElbowPlot(scRNA_harmony,ndims = 50) scRNA_harmony <- FindNeighbors(scRNA_harmony,reduction = 'harmony',dims = 1:20) # 这里的分辨率resolution选的是0.1 - 1 ,之后根据样本直接选择一个就行 scRNA_harmony <- FindClusters(scRNA_harmony,resolution = seq(from = 0.1,to = 1.0, by = 0.1)) scRNA_harmony <- RunUMAP(scRNA_harmony,dims = 1:20,reduction = 'harmony') scRNA_harmony <- RunTSNE(scRNA_harmony,dims = 1:20,reduction = 'harmony') # 保存降维聚类之后的数据 save(scRNA_harmony,file = 'scRNA_harmony.Rdata') 11.1 使用聚类树来确定分辨率代码如下（示例）：
library(clustree) clustree(scRNA_harmony) # 如何选择? 首先选择中间的,之后选择没有连线的,连线越多,被打乱就越多,所以我们选择0.3 11.2 用umap图查看harmony的整合情况代码如下（示例）：
# 我们这三个样本可以看出融合的还是比较好的 # 如果融合的情况不好的话,比如说肿瘤细胞有很强的抑制性,很难与其他细胞融合,这样我们也可以将它单独拿出来,研究它与其他细胞的特异性等等 DimPlot(scRNA_harmony,reduction = 'umap',group.by = 'orig.ident') + ggtitle('harmony') # 将分辨率改为0.3,分成12个cluster scRNA_harmony$RNA_snn_res.0.3 Idents(scRNA_harmony) <- 'RNA_snn_res.0.3' DimPlot(scRNA_harmony,reduction = 'umap') # 绘图 umap_integrated_1 <- DimPlot(scRNA_harmony,reduction = 'umap',group.by = 'orig.ident',label = T) umap_integrated_2 <- DimPlot(scRNA_harmony,reduction = 'tsne', label = T) umap_integrated_3 <- DimPlot(scRNA_harmony,reduction = 'umap', label = T) # 合并图片 integrated_plot <- CombinePlots(list(umap_integrated_1,umap_integrated_2,umap_integrated_3)) # 输出到画板 integrated_plot # 保存 ggsave('integrated_plot.png',integrated_plot,width = 10,height = 10) # 保存数据 save(scRNA_harmony,file = 'scRNA_harmony_2.Rdata') load('scRNA_harmony_2.Rdata') 三、一些常见问题的讨论 1.关于双细胞是否去除在单细胞RNA测序数据分析中，双细胞（双核细胞或细胞聚类中的双细胞）是指在一个细胞中包含两个或多个细胞核的情况。双细胞可能是由于细胞不完全溶解、细胞间粘连或细胞核分裂等原因引起的。处理双细胞的方法因数据特点和研究问题的不同而有所不同。一般来说，可以考虑以下几种处理双细胞的方法：
去除双细胞：将包含双细胞的细胞从数据集中去除。这种方法可以确保单细胞分析的准确性和可靠性，避免将两个不同的细胞混合在同一个细胞群集中。
将双细胞拆分为多个细胞：对于含有多个细胞核的双细胞，可以通过图像分析或细胞核识别算法将其拆分为多个细胞。然后，这些拆分后的细胞可以作为单独的细胞进行后续分析。
保留双细胞并进行分析：有些研究可能对双细胞感兴趣，或者双细胞可能提供了有关细胞类型和状态的额外信息。在这种情况下，可以选择保留双细胞并进行细胞类型鉴定、差异表达分析或其他相关分析。
最终的选择应基于具体研究的目的和数据特点，在进行单细胞分析之前，仔细检查和处理双细胞以确保数据质量和解读的准确性。
2.关于细胞周期的影响细胞周期是指细胞在不同阶段进行生长、DNA复制和分裂的过程。在单细胞RNA测序数据分析中，细胞周期的影响可以表现在以下几个方面：
基因表达变化：不同细胞周期阶段的细胞在基因表达水平上会有明显的差异。例如，在细胞进入有丝分裂阶段时，某些细胞周期相关基因的表达会显著上调。这种差异可能会掩盖其他细胞类型相关的差异，因此需要考虑在分析中对细胞周期进行校正或考虑。
细胞类型鉴定：细胞周期的变化可能会模糊细胞类型的辨识。同一种细胞类型在不同细胞周期阶段的基因表达模式可能会有所不同，导致在细胞类型鉴定时出现困惑。因此，在进行细胞类型鉴定时，需要注意是否将细胞周期纳入考虑，例如通过选择性排除或考虑细胞周期相关的基因。
数据解释：对于研究细胞功能和调控的研究，细胞周期的影响需要考虑在内。细胞周期的变化可能与特定功能相关的基因表达模式有关。因此，在解读单细胞RNA测序数据时，需要根据研究问题的不同，考虑细胞周期的影响并对结果进行解读。
是否需要去除细胞周期的影响取决于具体的研究目的和分析方法。在单细胞RNA测序数据分析中，有些研究关注的是细胞周期的变化，而有些研究则需要将细胞周期的影响排除或校正。如果研究关注细胞周期的影响，比如探究细胞周期与特定基因表达模式之间的关系，或者研究细胞周期对细胞功能和调控的影响，那么去除细胞周期的影响可能不是必需的。相反，细胞周期的变化本身可能是研究的重要部分，可以通过分析细胞周期相关基因的表达来获得相关结论。
然而，如果研究关注的是细胞类型鉴定、差异表达分析或对细胞功能的整体理解，那么去除细胞周期的影响可能是必要的。细胞周期的变化可能会掩盖其他细胞类型相关差异或基因表达的差异，导致结果的误解。在这种情况下，可以通过校正细胞周期相关基因的表达或选择性排除细胞周期相关的基因（如细胞周期调控基因）来消除细胞周期的影响。因此，是否需要去除细胞周期的影响取决于研究问题的具体需求和分析策略。在进行分析之前，需要明确研究的目的和分析的方向，以确定是否需要去除细胞周期的影响。
在了解背景的基础之上,我们也可以使用umap图等方式来可视化,从而得到细胞周期对分析的影响
代码如下（示例）：
s.genes <- cc.genes$s.genes g2m.genes <- cc.genes$g2m.genes scRNA_harmony <- CellCycleScoring(scRNA_harmony,s.features = s.genes,g2m.features = g2m.genes,set.ident = T) # plot DimPlot(scRNA_harmony,group.by = 'Phase') # 这里可以看出细胞周期的影响很小,几乎都是均匀分布的 # 这里是正常的,不用回归细胞周期 DimPlot(scRNA_harmony,group.by = 'Phase',reduction = 'pca') scRNA_harmony <- ScaleData(scRNA_harmony,vars.to.regress = c('S.Score','G2M.Score'), features = rownames(sc_example)) 3.提取count data代码如下（示例）：
# data 通常包含经过一些预处理（如标准化、缩放等）后的数据，它反映了经过一定转换后每个细胞中基因的表达量 # counts 则是原始的基因计数数据，也就是直接统计得到的每个细胞中各个基因的表达计数 # 1.counts # function1 scRNA_counts <- LayerData(scRNA_harmony,assay = 'RNA',layer = 'counts') scRNA_counts <- as.data.frame(scRNA_counts) # function2 scRNA_counts <- GetAssayData(scRNA_harmony,assay = 'RNA',layer = 'counts') scRNA_counts <- as.data.frame(scRNA_counts) # 2.data scRNA_data <- LayerData(scRNA_harmony,assay = 'RNA',layer = 'data') scRNA_data <- as.data.frame(scRNA_data) # 想要保存的话用write.table()函数 4.求每个cluster的平均表达量代码如下（示例）：
av_seurat <- AverageExpression(scRNA_harmony,assays = 'RNA',group.by = 'RNA_snn_res.0.3') av_seurat <- as.data.frame(av_seurat$RNA) av_seurat 5.伪bulk数据AggregateExpression代码如下（示例）：
av_seurat <- AggregateExpression(scRNA_harmony,assays = 'RNA',group.by = 'RNA_snn_res.0.3') av_seurat <- as.data.frame(av_seurat) av_seurat <- AggregateExpression(scRNA_harmony,assays = 'RNA',group.by = 'orig.ident') av_seurat <- as.data.frame(av_seurat) av_seurat <- AggregateExpression(scRNA_harmony,assays = 'RNA',group.by = c('RNA_snn_res.0.3','orig.ident')) av_seurat <- as.data.frame(av_seurat)伪bulk数据是一种将单细胞RNA测序数据聚合成表观上的批次数据的方法。它的作用包括：
细胞类型鉴定：通过将单细胞数据聚合成伪bulk数据，可以更好地识别不同细胞类型的表达模式，并确定哪些基因在不同细胞类型之间具有差异表达。
基因差异表达分析：由于单细胞数据通常具有较高的噪音和稀疏性，因此对于检测基因之间的差异表达可能存在挑战。通过将单细胞数据聚合成伪bulk数据，可以增加样本的数量从而提高差异表达分析的可靠性。
细胞群体分析：在某些情况下，我们可能更关注整个细胞群体的表达特征，而不是单个细胞。通过将单细胞数据聚合成伪bulk数据，可以更好地获取细胞群体的平均表达模式。
数据集集成：单细胞RNA测序数据可能来自不同的实验批次或研究组，因此可能存在一些技术差异。通过将单细胞数据聚合成伪bulk数据，可以减少技术差异的影响，使得不同数据集之间更容易进行比较和综合分析。
总之，伪bulk数据提供了一种将单细胞数据转化为表观上的批次数据的方法，它可以用于各种单细胞数据分析，包括细胞类型鉴定、差异表达分析、细胞群体分析和数据集集成等。
总结(代码资料领取处)本节内容在之前seurat v4的教程都很多都有提到,文章所用到的数据及代码已经给各位放在这里了,有需要自取!
链接：请点击这里哦
提取码：bui1
