1. 拿到数据后先别急着跑代码搞懂文件是第一步拿到单细胞转录组数据尤其是第一次接触的时候很多人会有点懵。公司或者合作方给过来的往往是一个压缩包解压后里面一堆文件夹和文件名字还都挺奇怪。别慌这第一步其实很简单就是搞清楚你手里到底有什么“食材”才能决定后面怎么做“菜”。我处理过很多项目发现不同测序平台给的数据格式差别挺大。最常见的就是10X Genomics的数据这也是目前市场占有率最高的平台。他们给的标准输出里最关键的就是那个叫filtered_feature_bc_matrix的文件夹。你打开它里面通常就三个压缩文件barcodes.tsv.gz、features.tsv.gz和matrix.mtx.gz。别看就三个文件它们构成了单细胞分析的基石。barcodes.tsv.gz里面是每个细胞的“身份证号”也就是条形码序列features.tsv.gz是基因列表告诉你测了哪些基因而matrix.mtx.gz这个文件最重要它是一个稀疏矩阵存储着每个细胞里每个基因的表达量也就是我们常说的UMI计数。这个矩阵非常“稀疏”因为一个细胞通常只表达一两万个基因中的几千个所以大部分值都是0用稀疏矩阵存储能节省大量空间。除了10X你可能会遇到BDBecton Dickinson平台的数据。他们处理完通常会给你一个CSV文件比如叫RSEC_MolsPerCell.csv。这个文件本质上也是一个表达矩阵行是基因列是细胞里面的数字是每个细胞的分子计数。文件通常很大用Excel是打不开的我们都是在R或者Python里直接读取。还有Smart-seq2技术的数据它更像传统的转录组测序给出的表达量往往是经过标准化的RPKM或TPM值而不是原始的UMI计数这一点在后续分析时要特别注意。所以第一步千万别跳过。花十分钟用几行简单的代码看看数据的基本信息比如用R的dim()看看矩阵有多大多少基因多少细胞用head()看看前几行数据长什么样。确认一下基因名称是不是你熟悉的符号比如是TP53还是ENSG00000141510这种编号细胞数量是不是和实验设计预期相符。这一步能帮你提前发现很多潜在问题比如文件损坏、格式错乱或者用了你不熟悉的基因组版本比如hg19还是hg38避免在分析到一半才发现基础数据有问题那可就白忙活了。2. 搭建你的分析环境R和Seurat的安装与配置工欲善其事必先利其器。单细胞分析的主流工具目前基本都集中在R语言的生态里尤其是Seurat这个包可以说是行业标准。所以第一步就是搭建好你的R分析环境。首先你需要安装R语言本身和RStudio。R是一个免费的统计计算环境而RStudio是一个让R用起来更舒服的集成开发环境IDE强烈建议两者都装。去R官网和RStudio官网下载对应你操作系统Windows、macOS或Linux的安装包一路下一步就行没什么坑。安装好之后打开RStudio我们首先要做的是配置软件包安装镜像。因为很多生物信息学包Bioconductor系列的服务器在国外直接下载可能会很慢甚至失败。在R控制台里运行下面两行代码换成国内的镜像源速度会快很多options(BioC_mirrorhttps://mirrors.tuna.tsinghua.edu.cn/bioconductor) options(repos c(CRANhttps://mirrors.tuna.tsinghua.edu.cn/CRAN/))接下来就是安装核心的分析包了。我习惯一次性把常用的都装上避免后面用到时再手忙脚乱。基础的数据处理和可视化包是必须的install.packages(Seurat) # 单细胞分析核心包 install.packages(dplyr) # 数据操作神器 install.packages(ggplot2) # 画图包 install.packages(stringr) # 字符串处理 install.packages(future) # 并行计算加快速度安装Seurat的时候可能会自动安装很多依赖包耐心等待即可。然后我们需要从Bioconductor安装一些专门用于生物信息学分析的包比如做差异表达、通路富集分析的install.packages(BiocManager) # Bioconductor的安装管理器 library(BiocManager) BiocManager::install(limma) # 差异表达分析 BiocManager::install(clusterProfiler) # GO/KEGG富集分析 BiocManager::install(org.Hs.eg.db) # 人类基因注释数据库 BiocManager::install(SingleR) # 自动细胞注释工具 BiocManager::install(celldex) # SingleR需要的参考数据集有时候安装会问你“Update all/some/none? [a/s/n]:”直接输入a然后回车。如果问“Do you want to install from sources the packages which need compilation? (Yes/no/cancel)”输入no回车。这样能避免从源码编译节省时间。最后还有一些好用的工具包在GitHub上比如检测双细胞Doublet的DoubletFinder做细胞间通讯分析的CellChat以及用于批次效应校正的harmony。我们可以用devtools来安装它们install.packages(devtools) library(devtools) devtools::install_github(chris-mcginnis-ucsf/DoubletFinder) devtools::install_github(sqjin/CellChat) devtools::install_github(immunogenomics/harmony)环境搭好之后我建议你新建一个R项目Project为你的单细胞分析单独创建一个工作目录把数据、脚本和结果都放在里面这样管理起来非常清晰也方便复现。3. 数据读入与初筛把原始数据变成Seurat对象环境准备好了数据也看明白了现在就要开始真正的分析了。第一步是把原始数据读入R并转换成Seurat对象。Seurat对象是一个容器它把表达矩阵、细胞元数据比如样本来源、线粒体基因比例等、分析结果如PCA坐标、聚类信息都打包在一起后续所有操作都围绕这个对象进行非常方便。对于10X Genomics的数据Seurat提供了非常方便的函数Read10X()和CreateSeuratObject()。假设你的数据文件夹路径是./data/filtered_feature_bc_matrix/操作如下library(Seurat) library(dplyr) # 设置工作目录到你的数据所在位置 setwd(/path/to/your/data/) # 读取10X数据 data_matrix - Read10X(data.dir ./filtered_feature_bc_matrix/) # 查看一下数据维度基因数 x 细胞数 dim(data_matrix) # 看看前5个细胞和前5个基因的表达情况 data_matrix[1:5, 1:5] # 创建Seurat对象 seurat_obj - CreateSeuratObject(counts data_matrix, project MyProject, # 给你的项目起个名 min.cells 3, # 一个基因至少在3个细胞中表达才保留 min.features 200) # 一个细胞至少检测到200个基因才保留这里的min.cells和min.features是两个重要的初筛参数。min.cells 3会过滤掉那些在绝大多数细胞里都不表达的基因这些基因可能是噪音对下游分析没帮助反而增加计算量。min.features 200则会过滤掉那些检测到的基因数太少的细胞这些细胞可能是空液滴没有细胞或者死细胞质量很差。这个阈值可以根据你的数据灵活调整比如对于高质量的数据可以设高一点如500对于捕获效率不高的数据可以设低一点。创建好对象后我们立刻要做一个关键的质量控制QC步骤计算每个细胞的线粒体基因比例。线粒体基因是细胞能量的来源但当细胞开始死亡或处于应激状态时细胞膜破裂线粒体RNA会大量泄露出来导致测序数据中线粒体基因的占比异常升高。因此线粒体基因比例是衡量细胞活性的一个非常重要的指标。# 计算线粒体基因百分比人类基因通常以MT-开头 seurat_obj[[percent.mt]] - PercentageFeatureSet(seurat_obj, pattern ^MT-) # 查看元数据现在多了三列nFeature_RNA基因数nCount_RNA总UMI数percent.mt线粒体基因百分比 head(seurat_objmeta.data) # 画个小提琴图直观感受一下数据质量分布 VlnPlot(seurat_obj, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3)通过这个图你可以快速判断数据的质量。健康的细胞其nFeature_RNA每个细胞检测到的基因数和nCount_RNA每个细胞的总UMI数应该在一个合理的范围内并且percent.mt线粒体基因百分比较低。如果看到有些细胞基因数极少比如小于100或者线粒体基因比例极高比如大于20%这些很可能就是需要过滤掉的低质量细胞或死细胞。4. 质量控制与细胞过滤把“坏”细胞踢出去看完QC图心里就有数了接下来就是动手过滤。这一步的目标是剔除低质量细胞保留高质量的数据用于下游分析。过滤的标准主要基于三个指标每个细胞检测到的基因数量nFeature_RNA、每个细胞的总分子数nCount_RNA以及线粒体基因比例percent.mt。怎么定过滤阈值呢没有绝对的金标准需要结合你的数据和生物学背景来判断。我通常的经验是基因数下限一般不低于200-300。低于这个数细胞包含的生物学信息太少。基因数上限通常不高于6000-7000。过高的基因数可能意味着这个“细胞”实际上是由两个或多个细胞粘在一起形成的“双细胞”或多细胞聚集体。线粒体基因比例一般设定在10%-20%以下。超过这个阈值细胞很可能正在死亡或已经死亡。我们可以用subset()函数来执行过滤# 执行过滤保留基因数在200-7000之间且线粒体基因比例低于20%的细胞 seurat_obj_filtered - subset(seurat_obj, subset nFeature_RNA 200 nFeature_RNA 7000 percent.mt 20) # 过滤前后对比一下细胞数量 cat(过滤前细胞数, ncol(seurat_obj), \n) cat(过滤后细胞数, ncol(seurat_obj_filtered), \n) cat(过滤掉了, ncol(seurat_obj) - ncol(seurat_obj_filtered), 个低质量细胞。\n)过滤掉多少细胞算正常这取决于样本类型和实验质量。对于组织解离做得好的样本可能只过滤掉5%-10%对于难以解离的组织或者新鲜度较差的样本过滤掉20%-30%也是常见的。如果过滤掉超过一半的细胞那就要回头检查实验环节是不是出问题了。过滤之后我强烈建议再画一次QC图确认一下过滤效果VlnPlot(seurat_obj_filtered, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3)你应该会看到小提琴图的“尾巴”代表异常值被砍掉了数据的分布更加集中和“健康”。此外还可以画一个相关性散点图检查基因数和总分子数之间是否呈良好的线性关系这通常是高质量数据的特征plot - FeatureScatter(seurat_obj_filtered, feature1 nCount_RNA, feature2 nFeature_RNA) print(plot)一个干净、高质量的数据集是后续所有分析可靠性的基础所以这一步千万不能马虎。花点时间调整阈值多看看数据分布确保你留下的细胞是真正能代表你生物学问题的“好”细胞。5. 数据标准化与特征基因选择让细胞之间可比过滤完细胞我们手里的表达矩阵还是原始的UMI计数。这种数据有两个特点一是不同细胞测序深度即总分子数不同直接比较不公平二是其方差会随着均值增大而增大。所以我们需要进行标准化目的就是消除技术噪音如测序深度差异让细胞之间的基因表达水平具有可比性。Seurat使用LogNormalize方法进行标准化这是最常用的一种。它的原理很简单首先将每个细胞的表达量除以该细胞的总分子数scale.factor默认为10000得到一个“每万分子计数”CP10K然后对这个值加1避免对0取对数最后取自然对数。# 数据标准化 seurat_obj_filtered - NormalizeData(seurat_obj_filtered, normalization.method LogNormalize, scale.factor 10000)标准化之后矩阵里的数值就从整数计数变成了连续的对数转换值。接下来我们要从成千上万个基因中找出那些在细胞间表达变化最大、最能区分不同细胞状态的基因这些基因被称为“高变基因”Highly Variable Features, HVGs。下游的降维和聚类分析主要就依赖于这些高变基因因为它们携带了最多的生物学差异信息。Seurat中找高变基因的函数是FindVariableFeatures()常用的方法是vst方差稳定变换。我们会选择排名前2000个左右的高变基因。# 寻找高变基因 seurat_obj_filtered - FindVariableFeatures(seurat_obj_filtered, selection.method vst, nfeatures 2000) # 查看找到的高变基因 top10_hvgs - head(VariableFeatures(seurat_obj_filtered), 10) print(top10_hvgs) # 可视化高变基因绘制基因的均值-方差关系图 plot - VariableFeaturePlot(seurat_obj_filtered) print(plot)在得到的图上高变基因会被高亮显示。你可以看到这些基因通常分布在散点图的右上部分意味着它们既有较高的平均表达水平又有较高的细胞间变异。选择2000个是一个经验值对于大多数数据集够用了。如果你的细胞数特别多比如超过10万或者想捕捉更细微的差异可以适当增加到3000或5000。找到高变基因后我们还需要进行一步“缩放”Scaling。这一步的目的是将所有基因的表达值进行线性变换使其均值为0方差为1。这样做是为了在后续的主成分分析PCA中让所有基因具有相同的权重避免高表达基因主导分析结果。同时我们还可以在这一步回归掉一些不想要的变异来源比如我们之前计算的线粒体基因比例。# 数据缩放并回归掉线粒体基因比例的影响 all_genes - rownames(seurat_obj_filtered) seurat_obj_filtered - ScaleData(seurat_obj_filtered, features all_genes, vars.to.regress percent.mt)执行完ScaleData后数据就准备好了。我们可以用GetAssayData函数查看缩放后的数据你会发现数值变成了有正有负的Z-score。至此数据预处理的核心步骤就完成了接下来我们将进入降维和可视化的世界。6. 线性降维与细胞分群看清细胞的“社会关系”现在我们有了一个经过标准化、缩放并聚焦于高变基因的数据集。但它的维度仍然很高2000个基因就是2000维我们人脑无法理解。降维的目的就是把数据从高维空间映射到低维空间比如2维或3维同时尽可能保留细胞之间的差异信息。最经典、最常用的线性降维方法就是主成分分析PCA。PCA会找到数据中变异最大的几个方向主成分每个主成分都是原始基因的线性组合。我们通常取前几十个主成分来进行后续分析。# 运行PCA基于之前找到的高变基因 seurat_obj_filtered - RunPCA(seurat_obj_filtered, features VariableFeatures(object seurat_obj_filtered)) # 查看前5个主成分中正负贡献最大的基因有助于理解PC的生物学意义 print(seurat_obj_filtered[[pca]], dims 1:5, nfeatures 5) # 可视化每个主成分的基因载荷Loading VizDimLoadings(seurat_obj_filtered, dims 1:2, reduction pca) # 绘制碎石图Elbow Plot帮助决定使用多少个主成分 ElbowPlot(seurat_obj_filtered, ndims 50)ElbowPlot生成的碎石图非常重要。横坐标是主成分序号纵坐标是每个主成分解释的方差百分比。图形通常会在某个点出现一个“拐点”elbow拐点之后的主成分解释的方差增量变得很小。一个常见的经验法则是选择拐点处的主成分数量。例如如果拐点在PC15附近那么我们可以选择前15-30个主成分用于下游分析。选择太少会丢失信息选择太多会引入噪音。确定好使用的主成分数量比如前30个PC后我们就可以基于这些PC来对细胞进行聚类了。聚类的目标是找到在基因表达模式上相似的细胞群体。Seurat使用一种基于图的方法首先根据细胞在PCA空间中的距离构建一个K近邻图FindNeighbors然后使用模块化优化算法如Louvain算法来识别图中的社区也就是细胞簇FindClusters。# 基于PCA结果构建细胞邻接图 seurat_obj_filtered - FindNeighbors(seurat_obj_filtered, dims 1:30) # 进行细胞聚类resolution参数控制分群的粗细程度 seurat_obj_filtered - FindClusters(seurat_obj_filtered, resolution 0.8)这里的resolution参数是关键它越大得到的细胞簇越多、越细越小则簇越少、越粗。对于初次分析我建议尝试一组分辨率如0.4, 0.6, 0.8, 1.0然后结合后续的标记基因和生物学知识来判断哪个分辨率最合理。你可以通过Idents(seurat_obj_filtered)来查看每个细胞被分配到了哪个簇。聚类完成后我们得到了细胞的分群信息但PCA坐标对于可视化来说还不够直观。我们需要进一步使用非线性降维方法将细胞投射到2维平面上最常用的就是t-SNE和UMAP。# 运行UMAP降维目前更流行能更好地保留全局结构 seurat_obj_filtered - RunUMAP(seurat_obj_filtered, dims 1:30) # 运行t-SNE降维传统方法局部结构保持好 seurat_obj_filtered - RunTSNE(seurat_obj_filtered, dims 1:30) # 可视化聚类结果 p1 - DimPlot(seurat_obj_filtered, reduction umap, label TRUE) p2 - DimPlot(seurat_obj_filtered, reduction tsne, label TRUE) p1 | p2 # 并排显示现在你就能看到一张漂亮的点图了每个点代表一个细胞颜色代表它所属的簇在UMAP/t-SNE空间上位置相近的细胞其基因表达谱也相似。这是你第一次“看到”自己数据中的细胞群体结构是不是很激动但这只是开始我们还需要知道这些簇到底代表什么类型的细胞。7. 寻找标记基因与手动细胞注释给细胞群“贴标签”细胞聚类给了我们一群一群的细胞但这些都是计算机根据数学相似性划分的我们需要给这些群赋予生物学意义这就是细胞注释。注释的第一步是为每个细胞簇寻找其特异的“标记基因”Marker Genes。标记基因是在某个细胞簇中高表达而在其他簇中低表达或不表达的基因。Seurat的FindAllMarkers()函数可以一次性为所有簇寻找标记基因。# 为所有聚类寻找标记基因只找在每个簇中高表达的基因 cluster_markers - FindAllMarkers(seurat_obj_filtered, only.pos TRUE, # 只保留高表达的 min.pct 0.25, # 基因至少在25%的该簇细胞中表达 logfc.threshold 0.25) # 对数倍变化阈值 # 查看结果按平均对数倍变化avg_log2FC降序排列 top_markers - cluster_markers %% group_by(cluster) %% slice_max(n 5, order_by avg_log2FC) print(top_markers, n50)结果是一个数据框包含了基因名、所属簇、平均表达量、在簇内和簇外的表达比例、p值等信息。avg_log2FC是关键的指标它越大说明这个基因在该簇中相对于其他簇的表达特异性越强。拿到标记基因列表后我们就需要像侦探一样根据已知的生物学知识来推断每个细胞簇的身份。这是一项需要经验的工作。下面我列一个常见细胞类型及其经典标记基因的速查表你可以对照着看细胞类型经典标记基因人类功能提示上皮细胞EPCAM, CDH1, KRT8, KRT18构成组织和器官的屏障癌症研究中的主角。T细胞CD3D, CD3E, CD8A (细胞毒性T), CD4 (辅助T)免疫系统的核心负责细胞免疫。B细胞CD19, MS4A1 (CD20), CD79A产生抗体负责体液免疫。浆细胞SDC1 (CD138), MZB1, JCHAIN, XBP1分泌抗体的终末分化B细胞。自然杀伤细胞NKG7, GNLY, PRF1, KLRD1天然免疫细胞能快速杀伤异常细胞。髓系细胞CD14 (单核/巨噬), FCGR3A (CD16非经典单核/巨噬), LYZ (溶菌酶)包括单核细胞、巨噬细胞、树突状细胞等是免疫系统的“清道夫”和“哨兵”。内皮细胞PECAM1 (CD31), VWF, CDH5 (VE-cadherin)构成血管内壁。成纤维细胞COL1A1, COL3A1, DCN, LUM, PDGFRA产生细胞外基质维持组织结构。周细胞RGS5, CSPG4 (NG2), PDGFRB包裹在微血管外壁的细胞。肥大细胞TPSAB1, CPA3, HDC参与过敏反应和炎症。有了这个表你就可以检查每个簇的top标记基因。比如如果一个簇高表达EPCAM和KRT8而不表达免疫细胞标记那它很可能是上皮细胞。如果一个簇高表达CD3D和CD8A那它很可能是CD8 T细胞。我们还可以用FeaturePlot()函数在UMAP/t-SNE图上直观地看某个标记基因的表达分布# 可视化几个关键标记基因的表达 FeaturePlot(seurat_obj_filtered, features c(EPCAM, CD3D, CD19, PECAM1, COL1A1), reduction umap, ncol 3)结合标记基因列表和可视化你可以开始手动给细胞簇重命名# 根据标记基因分析结果手动重命名细胞簇 new_cluster_names - c(0 CD8_T_Cells, 1 Fibroblasts, 2 Endothelial, 3 Epithelial, 4 Monocytes, 5 B_Cells, 6 NK_Cells, 7 CD4_T_Cells) # 应用新的命名 seurat_obj_filtered - RenameIdents(seurat_obj_filtered, new_cluster_names) # 查看重命名后的结果 DimPlot(seurat_obj_filtered, reduction umap, label TRUE)手动注释是金标准但非常依赖研究者的经验而且对于新的、未知的细胞类型可能无能为力。这时候自动注释工具就能帮上大忙了。8. 利用自动注释工具辅助判断自动细胞注释工具顾名思义就是利用已有的、经过人工精心注释的参考数据集通过算法来预测你数据中每个细胞的类型。这对于新手或者处理复杂组织如大脑、肿瘤微环境时特别有帮助可以提供一个强大的参考。这里我介绍两个最常用的R包SingleR和scHCL。SingleR的原理是将你的每个细胞与参考数据集的每个细胞进行比对找出最相似的“邻居”然后根据邻居的标签来投票决定你的细胞的类型。它内置了多个高质量的参考数据集比如HumanPrimaryCellAtlasData。library(SingleR) library(celldex) # 将Seurat对象的表达矩阵提取出来注意要用counts数据 test_data - GetAssayData(seurat_obj_filtered, slot counts) # 加载人类原代细胞图谱参考数据集 ref_data - HumanPrimaryCellAtlasData() # 运行SingleR进行预测 pred_results - SingleR(test test_data, ref ref_data, labels ref_data$label.main) # 使用主要的细胞类型标签 # 查看预测结果 table(pred_results$labels) # 将预测结果添加回Seurat对象的元数据中 seurat_obj_filtered$singleR_labels - pred_results$labels # 在UMAP图上用预测的标签着色与我们的聚类结果对比 DimPlot(seurat_obj_filtered, reduction umap, group.by singleR_labels, label TRUE)另一个工具是scHCL它基于人类细胞图谱Human Cell Landscape项目涵盖了更广泛的成人组织细胞类型。# 安装和加载scHCL # devtools::install_github(ggjlab/scHCL) library(scHCL) # 注意scHCL需要输入log归一化后的数据 test_data_log - GetAssayData(seurat_obj_filtered, slot data) # 运行scHCL预测这一步可能比较耗时 hcl_result - scHCL(scdata test_data_log, numbers_plot 1) # 查看预测的细胞类型 head(hcl_result$scHCL) # 同样将结果添加回元数据 seurat_obj_filtered$scHCL_labels - hcl_result$scHCL自动注释的结果并非100%准确尤其是当你的数据来自非常规组织或疾病状态时。但它能提供极其有价值的线索。我的做法是将自动注释结果、手动标记基因分析以及UMAP/t-SNE图的形态结合起来进行综合判断。比如SingleR预测某个簇是“T cells”scHCL也预测是“T cells”并且该簇高表达CD3D、CD8A等基因在UMAP图上和已知的免疫细胞区域挨在一起那么我们就可以比较有信心地将其注释为CD8 T细胞。自动注释和手动注释是相辅相成的。对于明确的、经典的细胞类型自动工具很准对于新的、过渡态的或疾病特异的细胞亚群则需要我们结合更多的标记基因和生物学背景进行深入挖掘。最终你会得到一份带有细胞类型标签的、可以用于下游深入分析如差异表达、通路分析、细胞通讯等的干净数据集。记得用saveRDS()函数保存你的劳动成果这个Seurat对象包含了从原始数据到注释结果的全部信息是后续所有分析的起点。