1. 从零开始理解单细胞数据与实战准备大家好我是老张在生物信息这个行当里摸爬滚打了十来年经手的单细胞项目少说也有几十个了。今天我想和你聊聊单细胞数据分析里最基础、也最关键的一步数据预处理。这就像是做菜前的备菜环节菜洗得干不干净、切得好不好直接决定了最后这盘菜能不能吃、好不好吃。很多新手朋友拿到数据后看着一堆以.mtx、.tsv.gz结尾的文件常常感到无从下手。别担心这篇文章就是为你准备的。我会用一个真实的乳腺癌单细胞数据集GSE161529作为例子手把手带你走完从原始数据到高质量细胞图谱的全过程。你不需要有高深的数学背景只要会基本的R语言操作跟着我的步骤就能得到一个干净、可靠的Seurat对象为后续的细胞亚群分析、差异表达等高级分析打下坚实基础。在开始敲代码之前我们得先搞清楚手里拿的是什么东西。单细胞测序数据特别是像10X Genomics这种主流平台产生的数据通常以三个核心文件的形式提供matrix.mtx.gz表达矩阵记录每个细胞里每个基因的计数、barcodes.tsv.gz细胞条形码每个条形码对应一个细胞和features.tsv.gz基因特征通常是基因ID和基因名。你可以把它们想象成一个巨大的Excel表格行是基因列是细胞表格里的数字就是每个细胞里每个基因被检测到的分子数UMI计数。我们的目标就是把这个“原始表格”读进R里把它变成一个Seurat能认识、能处理的“智能对象”。工欲善其事必先利其器。实战的第一步永远是搭建环境。你需要确保电脑上安装了R建议4.0以上版本和RStudio。接下来就是安装一系列必要的R包。Seurat是单细胞分析的“瑞士军刀”是绝对的核心。其他像dplyr、Matrix、ggplot2等则是数据处理和可视化的好帮手。我建议你直接在RStudio的控制台里运行下面的安装命令。如果遇到网络问题可以尝试切换镜像源。# 安装CRAN上的包 install.packages(c(Seurat, dplyr, Matrix, ggplot2, cowplot, plyr)) # 安装Bioconductor上的包 if (!require(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(c(scater, scran)) # 安装GitHub上的包如Harmony用于批次校正 install.packages(devtools) devtools::install_github(immunogenomics/harmony)安装过程可能会花点时间喝杯咖啡等待一下就好。安装完成后记得用library()函数把它们都加载进来。我习惯在脚本开头一次性加载所有需要的包并把工作环境清空这样能避免之前运行的变量干扰本次分析。# 加载所有需要的R包 library(Seurat) library(dplyr) library(Matrix) library(ggplot2) library(cowplot) library(scater) library(scran) library(harmony) library(plyr) # 清空当前环境确保分析纯净 rm(list ls())2. 数据读取与整合把原始文件变成Seurat对象环境准备好了数据也下载到了本地假设你解压后放在./GSE161529_RAW文件夹里接下来就是真刀真枪地读数据了。这一步看似简单却最容易出岔子比如文件路径不对、文件名不匹配等等。我的经验是先用list.files()函数看看文件夹里到底有哪些文件做到心中有数。我们面对的数据集往往包含多个样本。以乳腺癌数据为例可能同时有癌组织和癌旁组织或者不同病人的样本。我们需要分别读取每个样本再把它们巧妙地合并起来。下面这段代码是我在实际项目中反复打磨过的逻辑清晰且健壮。它首先识别出所有样本的矩阵和条形码文件然后通过一个循环逐个读取并创建Seurat对象最后使用merge函数将它们整合成一个大的对象。这里的add.cell.ids参数至关重要它会给每个样本的细胞ID加上前缀比如Sample1_AAACAT...这样在后续分析中我们永远能追溯到某个细胞来自哪个样本避免了混淆。# 1. 设置数据路径并列出文件 file.dir - ./GSE161529_RAW all.files - list.files(file.dir) # 2. 区分矩阵文件、条形码文件和特征文件 # 注意文件名可能因数据集而异关键是匹配模式 matrix.files - all.files[grepl(matrix.mtx.gz$, all.files)] barcode.files - all.files[grepl(barcodes.tsv.gz$, all.files)] # 假设只有一个全局的特征文件 feature.file - all.files[grepl(features.tsv.gz$, all.files)] # 3. 提取样本名假设文件名以样本ID开头后接下划线或点 # 例如GSM4901234_matrix.mtx.gz - GSM4901234 sample.names - unique(gsub(_(matrix.mtx.gz|barcodes.tsv.gz).*, , c(matrix.files, barcode.files))) # 4. 读取特征基因信息 features - read.table(gzfile(paste(file.dir, feature.file, sep/)), header FALSE, stringsAsFactors FALSE) # 通常第二列是基因符号gene symbol我们用它作为行名 gene.names - features[, 2] # 5. 定义一个函数用于读取单个样本并创建Seurat对象 create_seurat_for_sample - function(sample.id) { message(正在处理样本: , sample.id) # 构建完整的文件路径 mtx.path - paste(file.dir, grep(paste0(sample.id, .*matrix.mtx.gz), matrix.files, valueTRUE), sep/) barcode.path - paste(file.dir, grep(paste0(sample.id, .*barcodes.tsv.gz), barcode.files, valueTRUE), sep/) # 读取稀疏矩阵和条形码 counts - readMM(gzfile(mtx.path)) barcodes - read.table(gzfile(barcode.path), headerFALSE, stringsAsFactorsFALSE)[,1] # 设置矩阵的行名基因和列名细胞 rownames(counts) - gene.names colnames(counts) - barcodes # 创建Seurat对象project参数标记样本来源 seurat.obj - CreateSeuratObject(counts counts, project sample.id, min.cells 3, min.features 200) return(seurat.obj) } # 6. 循环处理所有选定的样本并存入列表 # 你可以通过修改select.samples来选择部分样本进行分析 select.samples - sample.names # 这里处理所有样本你也可以指定如 c(GSM4901234, GSM4901235) seurat.list - lapply(select.samples, create_seurat_for_sample) names(seurat.list) - select.samples # 7. 合并所有样本的Seurat对象 # 如果只有单个样本这一步可以跳过 if (length(seurat.list) 1) { merged.seurat - merge(seurat.list[[1]], y seurat.list[-1], add.cell.ids select.samples, project BRCA_Project) } else { merged.seurat - seurat.list[[1]] } # 8. 查看合并后的对象基本信息 print(merged.seurat)运行完这段代码后在控制台你会看到类似“An object of class Seurat... 13714 features across 5000 samples”的输出。这表示成功创建了一个包含13714个基因和5000个细胞的Seurat对象。你可以用head(merged.seuratmeta.data)查看前几个细胞的元信息确认样本来源orig.ident列已正确标记。这一步的稳健完成是整个分析流程的坚实起点。3. 核心预处理质控、归一化与特征选择拿到原始的Seurat对象后我们面对的是一个包含所有检测到的细胞和基因的“毛坯房”。里面不可避免地混入了一些“不合格产品”比如死细胞RNA含量极低、破损细胞线粒体基因比例过高以及双细胞或多细胞一个液滴里包裹了多个细胞导致基因数异常高。质量控制QC的目的就是把这些不合格的细胞剔除掉保证我们分析的都是高质量的单个细胞。那么如何判断一个细胞是好是坏呢Seurat对象的元数据meta.data里自动计算了几个关键指标nCount_RNA一个细胞的总UMI数反映测序深度、nFeature_RNA一个细胞中检测到的唯一基因数反映细胞复杂度。我们还需要手动计算两个指标mitoRatio线粒体基因占比和log10GenesPerUMI每个UMI对应的基因数反映数据的“效率”。线粒体基因占比高通常意味着细胞状态不好如凋亡或破损而log10GenesPerUMI过低则可能暗示大量重复序列或技术噪音。我通常先用VlnPlot函数可视化这些指标的分布直观感受一下数据的整体质量和需要设定的阈值。下面这段代码会生成一个组合图帮你做决策。# 计算线粒体基因和核糖体基因比例 merged.seurat[[percent.mt]] - PercentageFeatureSet(merged.seurat, pattern ^MT-) merged.seurat[[percent.rb]] - PercentageFeatureSet(merged.seurat, pattern ^RP[SL]) # 计算log10(基因数/UMI数) merged.seurat$log10GenesPerUMI - log10(merged.seurat$nFeature_RNA) / log10(merged.seurat$nCount_RNA) # 可视化四个QC指标 qc_plots - VlnPlot(merged.seurat, features c(nFeature_RNA, nCount_RNA, percent.mt, log10GenesPerUMI), pt.size 0.1, # 点的大小设为0可隐藏点 ncol 4) theme(axis.text.x element_text(angle 0)) print(qc_plots)看图说话我们就能设定合理的过滤阈值了。对于大多数10X数据我的经验阈值是nFeature_RNA在200-6000之间nCount_RNA在500-30000之间percent.mt小于20%log10GenesPerUMI大于0.8。当然这需要根据你的数据分布灵活调整。设定好阈值后用subset函数进行过滤。# 应用QC阈值过滤细胞 filtered.seurat - subset(merged.seurat, subset nFeature_RNA 200 nFeature_RNA 6000 nCount_RNA 500 nCount_RNA 30000 percent.mt 20 log10GenesPerUMI 0.8) message(paste(过滤前细胞数:, ncol(merged.seurat))) message(paste(过滤后细胞数:, ncol(filtered.seurat))) message(paste(过滤掉了, ncol(merged.seurat)-ncol(filtered.seurat), 个低质量细胞))质控之后就是数据转换的“三部曲”归一化Normalization、特征选择Feature Selection和缩放Scaling。归一化是为了消除不同细胞间测序深度差异带来的影响最常用的是LogNormalize方法。特征选择是为了找出那些在不同细胞间表达变化最大高变异的基因这些基因往往更能驱动细胞类型的区分Seurat默认找2000个这样的基因。缩放过去常叫标准化则是将每个基因的表达值进行线性转换使其均值为0方差为1这样后续的PCA等基于距离的算法才能公平地对待所有基因。# 1. 归一化使用LogNormalize方法默认 filtered.seurat - NormalizeData(filtered.seurat, normalization.method LogNormalize, scale.factor 10000) # 默认值可理解为将每个细胞的UMI总数缩放到10000 # 2. 寻找高变基因选择2000个变异最大的基因用于下游分析 filtered.seurat - FindVariableFeatures(filtered.seurat, selection.method vst, # 方差稳定变换法效果稳定 nfeatures 2000) # 可视化高变基因 top10 - head(VariableFeatures(filtered.seurat), 10) variable_feature_plot - VariableFeaturePlot(filtered.seurat) LabelPoints(plot variable_feature_plot, points top10, repel TRUE) # 3. 缩放数据为所有基因或高变基因进行线性变换 filtered.seurat - ScaleData(filtered.seurat, features rownames(filtered.seurat)) # 对所有基因进行缩放耗时但全面 # 如果数据量大想节省时间可以只缩放高变基因features VariableFeatures(filtered.seurat)这里有个小技巧ScaleData函数有一个强大的vars.to.regress参数可以用来校正一些我们不关心的技术变异来源比如细胞周期的影响、线粒体基因比例、甚至批次效应如果后续不用Harmony等专门方法。例如我们可以通过CellCycleScoring函数为每个细胞分配S期或G2M期分数然后在缩放时将其回归掉。# 加载细胞周期基因集Seurat内置 s.genes - cc.genes$s.genes g2m.genes - cc.genes$g2m.genes # 为每个细胞计算细胞周期分数 filtered.seurat - CellCycleScoring(filtered.seurat, s.features s.genes, g2m.features g2m.genes, set.ident FALSE) # 在缩放时回归掉细胞周期分数和线粒体基因比例的影响 filtered.seurat - ScaleData(filtered.seurat, vars.to.regress c(S.Score, G2M.Score, percent.mt), features rownames(filtered.seurat))4. 降维、聚类与批次校正看清细胞的“群落”经过前面的处理我们的数据已经干净多了但维度依然极高上万个基因。我们需要通过降维Dimensionality Reduction来捕捉数据中最主要的变化模式同时大幅降低计算复杂度。主成分分析PCA是这一步的标准工具。它会把数据投射到一系列新的坐标轴主成分PC上其中前几个PC就包含了数据的大部分信息。# 运行PCA分析默认使用前面找到的高变基因 filtered.seurat - RunPCA(filtered.seurat, features VariableFeatures(object filtered.seurat), npcs 50, # 计算50个主成分通常足够 verbose FALSE) # 可视化PCA结果 # 1. 查看每个主成分对方差的解释度 ElbowPlot(filtered.seurat, ndims 50) # 寻找“拐点”通常拐点后的PC包含的生物学信息变少 # 2. 查看指定基因在PC上的载荷loading VizDimLoadings(filtered.seurat, dims 1:4, reduction pca) # 3. 用前两个PC画散点图 DimPlot(filtered.seurat, reduction pca)ElbowPlot图非常重要它帮助我们决定在后续分析中该使用多少个PC。通常我们选择“拐点”处对应的PC数。比如拐点在PC15附近那么我们可能选择前15或20个PC来进行下游的聚类和UMAP可视化。接下来是聚类Clustering。简单说就是把表达谱相似的细胞归到同一个组里。Seurat基于共享最近邻SNN图算法进行聚类。你需要理解两个关键参数dims使用前多少个PC和resolution分辨率。分辨率越高得到的聚类数越多细胞分得越细分辨率越低聚类数越少细胞被归为更粗的类别。对于初次分析可以从0.4到1.2之间尝试几个值。# 基于PCA结果构建细胞间的邻接关系图 filtered.seurat - FindNeighbors(filtered.seurat, dims 1:20, # 这里假设我们根据肘部图选择了前20个PC reduction pca) # 进行聚类分析resolution是关键参数 filtered.seurat - FindClusters(filtered.seurat, resolution 0.5, # 尝试0.4, 0.6, 0.8等值 algorithm 1, # Louvain算法 random.seed 42) # 设置随机种子保证结果可重复 # 查看聚类结果会新增一个meta.data列叫seurat_clusters table(filtered.seuratmeta.data$seurat_clusters)现在我们有了聚类结果但PCA图二维很难直观展示细胞的群落关系。这时就需要UMAP或t-SNE这类非线性降维方法将高维数据映射到2D或3D空间便于我们肉眼观察。UMAP是目前更主流的选择因为它运行更快且能更好地保留全局数据结构。# 运行UMAP降维基于之前PCA的结果 filtered.seurat - RunUMAP(filtered.seurat, dims 1:20, # 必须与FindNeighbors使用的dims一致 reduction pca, seed.use 42, # 设置随机种子 n.neighbors 30, # 邻居数默认30对于大样本可适当增加 min.dist 0.3) # 控制点分布的紧密程度默认0.3值越小簇越紧 # 可视化UMAP按聚类着色 cluster_umap - DimPlot(filtered.seurat, reduction umap, label TRUE, # 在图上显示聚类编号 label.size 6, pt.size 0.5) ggtitle(UMAP plot colored by Seurat clusters) print(cluster_umap)如果你的数据来自多个样本、多个批次比如不同实验日期、不同测序通道批次效应可能会让相同类型的细胞因为技术原因而分开或者让不同类型的细胞因为批次相似而聚在一起。这时就需要进行批次校正Batch Correction。Harmony是一个整合在Seurat流程中非常好用的工具。使用起来非常简单在运行PCA之后调用RunHarmony函数指定批次变量比如样本IDorig.ident然后用Harmony降维的结果替代PCA的结果进行后续的找邻居、聚类和UMAP。# 假设我们的批次信息存储在orig.ident中即样本ID # 在运行PCA之后执行Harmony整合 filtered.seurat - RunHarmony(filtered.seurat, group.by.vars orig.ident, # 指定批次变量 reduction pca, assay.use RNA, project.dim FALSE) # 使用Harmony降维的结果reduction harmony进行后续分析 filtered.seurat - FindNeighbors(filtered.seurat, reduction harmony, dims 1:20) filtered.seurat - FindClusters(filtered.seurat, resolution 0.5) filtered.seurat - RunUMAP(filtered.seurat, reduction harmony, dims 1:20) # 比较校正前后的UMAP图 p1 - DimPlot(filtered.seurat, reduction umap, group.by orig.ident, pt.size0.5) ggtitle(Before Harmony (colored by batch)) p2 - DimPlot(filtered.seurat, reduction umap, group.by seurat_clusters, pt.size0.5) ggtitle(After Harmony (colored by cluster)) cowplot::plot_grid(p1, p2, ncol2)5. 细胞类型注释为细胞群落贴上“身份标签”得到一个个细胞簇Cluster后最激动人心也最具挑战性的环节来了细胞类型注释Cell Annotation。简单说就是回答“第0簇是什么细胞第1簇又是什么细胞”的问题。注释的准确性直接决定了整个生物学故事能否讲对。新手常犯的错误是只看一两个标记基因就下结论这很容易误判。我推荐结合多种方法进行交叉验证。方法一基于已知标记基因Manual Annotation。这是最经典的方法需要你查阅领域内的文献或权威数据库如CellMarker整理出不同细胞类型的经典标记基因。然后在数据中检查这些基因在各个簇中的表达情况。Seurat提供了FeaturePlot可视化基因在UMAP图上的表达和VlnPlot/DotPlot展示基因在不同簇的表达水平或比例等强大工具。# 定义一组乳腺癌微环境中常见的细胞类型标记基因示例需根据你的研究背景调整 marker_genes - list( T细胞 c(CD3D, CD3E, CD8A, CD4), B细胞 c(CD79A, MS4A1, CD19), 髓系细胞 c(CD68, CD14, FCGR3A), 成纤维细胞 c(DCN, COL1A1, FAP), 内皮细胞 c(PECAM1, VWF, CDH5), 上皮细胞/癌细胞 c(EPCAM, KRT8, KRT18) ) # 1. 用FeaturePlot查看关键标记基因的空间表达 FeaturePlot(filtered.seurat, features unlist(marker_genes)[1:6], # 先看前6个 reduction umap, order TRUE, # 将高表达细胞点画在上面 pt.size 0.5, ncol 3) # 2. 用DotPlot更系统地比较所有标记基因在所有簇中的平均表达和表达比例 DotPlot(filtered.seurat, features marker_genes, cols c(lightgrey, blue), # 颜色梯度 dot.scale 6) RotatedAxis()方法二使用参考数据集进行自动注释Automatic Annotation。随着单细胞数据库的丰富现在有很多工具可以利用已注释好的参考数据集通过算法自动为你注释。比如SingleR、scCATCH、CellAssign等。这里以SingleR为例它需要从celldex包加载一个参考数据集比如人类免疫细胞参考数据集HumanPrimaryCellAtlasData。# 安装并加载必要的包 BiocManager::install(c(SingleR, celldex)) library(SingleR) library(celldex) # 加载人类初级细胞图谱参考数据 ref - HumanPrimaryCellAtlasData() # 提取我们数据的表达矩阵使用归一化后的数据 test.data - GetAssayData(filtered.seurat, slotdata) # log-normalized counts # 运行SingleR进行注释 pred - SingleR(test test.data, ref ref, labels ref$label.main) # 使用粗粒度的细胞类型标签 # 将注释结果添加到Seurat对象的元数据中 filtered.seurat$singleR_labels - pred$labels # 可视化自动注释的结果 DimPlot(filtered.seurat, reduction umap, group.by singleR_labels, label TRUE) ggtitle(Cell annotation by SingleR)方法三差异表达基因DEG辅助注释。有时候某个簇没有非常特异的经典标记基因或者它可能是一个新的、未被很好定义的亚群。这时我们可以找出这个簇相对于其他所有簇特异性高表达的基因即该簇的差异表达基因然后去查这些基因的功能从而推断其身份。FindAllMarkers函数可以一键完成所有簇的差异分析。# 找出每个簇的差异表达基因与其余所有簇相比 # 这步计算量较大可能需要一些时间 cluster_markers - FindAllMarkers(filtered.seurat, only.pos TRUE, # 只保留在该簇中高表达的基因 min.pct 0.25, # 基因至少在25%的该簇细胞中表达 logfc.threshold 0.25) # 表达量差异倍数阈值 # 查看第0簇的前10个标记基因 top10_cluster0 - cluster_markers %% filter(cluster 0) %% top_n(n 10, wt avg_log2FC) print(top10_cluster0) # 我们可以用这些新发现的标记基因结合文献来更精确地注释细胞类型在实际项目中我通常会将这三种方法结合起来。先用SingleR等工具得到一个初步的、自动化的注释结果作为参考。然后仔细查看经典标记基因和差异表达基因的表达模式手动调整和确认每个簇的标签。特别是对于肿瘤样本癌细胞通常表达上皮标记如EPCAM的异质性很大可能需要结合拷贝数变异CNV分析来区分恶性细胞和正常上皮细胞。注释是一个迭代和需要生物学知识判断的过程不要期望一键得到完美结果。6. 流程封装与结果保存打造可复用的分析脚本经过上面五个步骤我们已经完成了一个完整的单细胞数据预处理流程。但在实际科研中我们经常需要分析多个数据集或者用不同的参数反复尝试。每次都从头写一遍代码既低效又容易出错。因此将核心流程封装成函数是提升分析效率和可重复性的关键一步。下面我把我常用的预处理流程封装成一个函数run_sc_preprocess它集成了质控、归一化、降维、聚类和基础可视化。# 单细胞数据预处理一站式函数 # # param seurat_obj 输入的原始Seurat对象 # param project_name 项目名称用于绘图标题和保存文件 # param qc_params 质控参数列表包含nFeature_min, nFeature_max, nCount_min, nCount_max, mt_percent_max, log10GPU_min # param do_harmony 是否进行Harmony批次校正 # param batch_var 进行Harmony校正时指定批次变量的列名在meta.data中 # param pcs_use 用于聚类和UMAP的主成分数量 # param cluster_resolution 聚类分辨率参数 # param seed 随机种子保证结果可重复 # # return 返回一个处理完成的Seurat对象 run_sc_preprocess - function(seurat_obj, project_name MyProject, qc_params list(nFeature_min200, nFeature_max6000, nCount_min500, nCount_max30000, mt_percent_max20, log10GPU_min0.8), do_harmony FALSE, batch_var NULL, pcs_use 20, cluster_resolution 0.5, seed 42) { set.seed(seed) message( 开始预处理流程, project_name) # --- 步骤1: 计算QC指标 --- message(1. 计算质量控制指标...) seurat_obj[[percent.mt]] - PercentageFeatureSet(seurat_obj, pattern ^MT-) seurat_obj$log10GenesPerUMI - log10(seurat_obj$nFeature_RNA) / log10(seurat_obj$nCount_RNA) # --- 步骤2: 质控过滤 --- message(2. 过滤低质量细胞...) n_before - ncol(seurat_obj) seurat_obj - subset(seurat_obj, subset nFeature_RNA qc_params$nFeature_min nFeature_RNA qc_params$nFeature_max nCount_RNA qc_params$nCount_min nCount_RNA qc_params$nCount_max percent.mt qc_params$mt_percent_max log10GenesPerUMI qc_params$log10GPU_min) n_after - ncol(seurat_obj) message(paste( 保留细胞:, n_after, /, n_before, sprintf((过滤掉 %.1f%%), 100*(n_before-n_after)/n_before))) # --- 步骤3: 归一化、找高变基因、缩放 --- message(3. 归一化与特征缩放...) seurat_obj - NormalizeData(seurat_obj, normalization.method LogNormalize, scale.factor 10000) seurat_obj - FindVariableFeatures(seurat_obj, selection.method vst, nfeatures 2000) seurat_obj - ScaleData(seurat_obj, features rownames(seurat_obj)) # --- 步骤4: 降维 (PCA) --- message(4. 线性降维 (PCA)...) seurat_obj - RunPCA(seurat_obj, features VariableFeatures(seurat_obj), npcs 50, verbose FALSE) # --- 步骤5: 批次校正 (可选) --- if (do_harmony !is.null(batch_var)) { message(5. 使用Harmony进行批次校正...) require(harmony) seurat_obj - RunHarmony(seurat_obj, group.by.vars batch_var, reduction pca, assay.use RNA, project.dim FALSE) reduction_use - harmony dims_use - 1:pcs_use } else { reduction_use - pca dims_use - 1:pcs_use } # --- 步骤6: 聚类 --- message(6. 细胞聚类...) seurat_obj - FindNeighbors(seurat_obj, reduction reduction_use, dims dims_use) seurat_obj - FindClusters(seurat_obj, resolution cluster_resolution, random.seed seed) # --- 步骤7: 非线性降维 (UMAP) --- message(7. UMAP可视化...) seurat_obj - RunUMAP(seurat_obj, reduction reduction_use, dims dims_use, seed.use seed, n.neighbors 30, min.dist 0.3) # --- 步骤8: 生成质控报告图 --- message(8. 生成可视化图表...) # 质控指标小提琴图 qc_vln - VlnPlot(seurat_obj, features c(nFeature_RNA, nCount_RNA, percent.mt), pt.size 0, ncol 3) plot_annotation(title paste(project_name, - QC Metrics (Post-filtering)), theme theme(plot.title element_text(hjust 0.5))) # 聚类UMAP图 cluster_umap - DimPlot(seurat_obj, reduction umap, label TRUE, pt.size 0.5) ggtitle(paste(project_name, - Clusters (Resolution, cluster_resolution, ))) # 打印图形 print(qc_vln) print(cluster_umap) message( 预处理流程完成) return(seurat_obj) } # 使用示例假设你的原始对象叫raw_seurat processed_seurat - run_sc_preprocess(seurat_obj raw_seurat, project_name BRCA_GSE161529, qc_params list(nFeature_min300, nFeature_max5000, nCount_min800, nCount_max25000, mt_percent_max15, log10GPU_min0.85), do_harmony TRUE, batch_var orig.ident, pcs_use 20, cluster_resolution 0.6, seed 2024)最后别忘了将辛苦处理好的Seurat对象保存下来。RDS格式是R语言存储单个对象的原生格式它能完整保留对象的所有结构、数据和元数据下次用readRDS函数读取后可以直接继续分析非常方便。# 保存最终的Seurat对象 saveRDS(object processed_seurat, file ./processed_BRCA_GSE161529_seurat.rds) # 同时也建议将关键的元数据如细胞注释信息单独保存为表格方便在其他软件中查看 write.csv(processed_seuratmeta.data, file ./processed_BRCA_GSE161529_metadata.csv, row.names TRUE)预处理后的高质量细胞图谱就像一份精炼过的原材料已经为后续的深入分析做好了准备。你可以用它来探索不同细胞亚群之间的差异表达基因进行细胞通讯分析或者拟时序分析来推断细胞分化轨迹。每一步的代码和参数选择都蕴含着对数据的理解多动手尝试多看看结果的分布你会越来越有感觉。记住预处理没有绝对的金标准最适合你数据的参数往往藏在数据本身的故事里。