【大瓜】颜宁：明年将不再接受国内任何个人或单位的邀请做学术报告；同时非公务邮件、快递一概拒收

非常简单，适合新手学习。作者就是对7个HGSOC样本和5个正常卵巢样本进行了单细胞测序：

🔽①长按下方二维码关注🔽

②输入关键词：7787

OK，我们开始操作，先下载这12个样本的单细胞测序数据，作者说数据已经传到GEO数据库了，号码是184880:

那我们打开生物医学之家（www.swyxzj.com）网站后进入GEO：

这个页面的每一个字都要细读，非常重要，作者怎么提取的单细胞，怎么做的测序，一定要细读，如下下载单细胞测序原始数据：

下载后解压我们可以看到36个文件夹：

每个GSM是一个样本，每个样本是3个文件夹（这是单细胞测序最常见的文件格式，三文件格式）。一共12个样本，所以是36个文件夹。随后我会自己新建文件夹，7个癌症样本和5个正常样本，如下。

每个样本下都是三文件形式：

随后用R语言读取数据，进行单细胞测序数据分析。首先是加载需要的包，将运行设置到最大（灵活的节省时间）：

然后是读取数据，构建seurat单细胞处理对象(后续所有的操作都是基于这个对象进行处理，质控降维等等)，并对正常对照的细胞都加上ctrl前缀分组和癌症组都加上exp前缀分组（这是万能代码，只要你以后的数据也是这样放置，三文件在一个文件夹中，都可以这样读取）：

1. 读取数据

这个时候我们已经过滤基因数量小于300的细胞（原文说是200，但是我试过之几个参数后发现300才能对上原文说的59342个细胞）。随后就是非常常规的单细胞数据分析处理流程，最终进行细胞降维得到细胞群，那具体处理流程如下：

2.1 质量控制（QC）过滤前：检查线粒体基因

为什么要做质量控制（QC）分析？

1. 评估数据质量：

• 线粒体基因（percent.mt）：高比例的线粒体基因表达通常表明细胞可能已经经历了应激或凋亡。

• 核糖体基因（percent.ribo）：核糖体基因表达的异常比例可能反映了特定的生物学特性或技术噪声。

• 血红蛋白基因（percent.hb）：某些样本中，血红蛋白表达高可能是血液污染的标志。

• nFeature_RNA 和 nCount_RNA 的分布可以帮助识别质量不佳的细胞，例如双重细胞（细胞核外有过多基因）、细胞碎片或低捕获效率的细胞。

• 通过 QC 图表，可以直观地观察指标分布，并据此设置合理的过滤阈值，例如线粒体比例阈值或基因数目范围。

结果如何解读？

1. 小提琴图：

• nFeature_RNA 和 nCount_RNA：

• 如果某些样本的这些值明显偏低，可能是由于文库构建效率较低或捕获效率低。

• 如果某些样本的值过高，可能是因为双重细胞或其他异常细胞。

• percent.mt、percent.ribo 和 percent.hb：

• 高比例线粒体基因（percent.mt > 15%-20%）通常表明细胞正在经历凋亡，需要过滤。

• 核糖体基因（percent.ribo）通常有一个正常范围，如果某些细胞偏高或偏低，可以根据实验背景进行解释。

• 血红蛋白基因（percent.hb）的高比例可能提示血液污染，尤其是在非血液来源的样本中。

2. 散点图：

• nCount_RNA vs percent.mt：

• 理想情况下，没有显著的相关性。如果有高相关性，可能是由于捕获效率异常导致线粒体基因的相对比例上升。

• nCount_RNA vs nFeature_RNA：

• 通常有较强的正相关性，表明捕获到的 RNA 数目与检测到的基因数目一致。

• 如果某些点偏离趋势，可能是由于双重细胞或低质量细胞。

当前结果解读：

从结果来看：

1. percent.mt 与 nCount_RNA 散点图：

• 图中显示了一些细胞具有非常高的线粒体比例（超过 50%-75%）。这些细胞可能是坏死或凋亡细胞，建议过滤。

• 散点呈现较好的正相关性，说明大部分细胞的质量较高。但右侧可能存在少量偏离趋势的高表达细胞，可能需要进一步检查。

通过这些图表，你可以根据实验背景设置合理的过滤参数。

2.2 过滤细胞：保留基因数小于10000的细胞，同时线粒体基因转录小于40%的细胞:

过滤后的 QC 结果解读

1. 数据过滤后的目标：

1. 基因数量限制：

• 过滤掉了检测到基因数量（nFeature_RNA）小于300的细胞（低质量或死细胞）和大于10000的细胞（可能是双细胞或异常高捕获）。

• 剩下的细胞应该反映生物学特征，而不是技术噪声。

• 设置了 percent.mt < 40.0% 的阈值（参考原文，一般也是这个参数），排除了高比例线粒体基因表达的细胞（可能是凋亡或应激细胞）。

• 剩下的细胞线粒体基因表达比例处于合理范围内。

2. 过滤后散点图解读：

• 左图（nCount_RNA vs percent.mt）：

• 现象：大部分细胞的线粒体基因比例已经集中在40%以下（过滤条件起作用），并且没有显著的线粒体基因比例与UMI数量的相关性。

• 解读：这表明过滤后的数据没有包含高线粒体比例的异常细胞，细胞质量较高。

• 右图（nCount_RNA vs nFeature_RNA）：

• 现象：大部分细胞呈现出正相关趋势，UMI计数较高的细胞对应检测到更多的基因。

• 解读：说明过滤后数据的质量较好，检测到的基因数量和捕获到的RNA分子数量具有一致性。

3. 过滤后细胞数量：

• 过滤前细胞数：64030。

• 过滤后细胞数：59324(与原文一致)。

• 说明：过滤掉了线粒体基因比例过高或基因检测数异常的细胞，占比约为 8%。

4. 过滤后的意义：

1. 去除了可能是凋亡细胞、双细胞或污染的异常数据，减少了技术噪声的干扰。

2. 提高了后续分析（例如降维、分群）的可信度，使得数据更能反映真实的生物学特性。

3. 对于后续步骤（例如差异分析、通路富集分析），可以避免因异常细胞影响而导致的假阳性结果。

这里我们发现细胞数量与原文一致，说明我们的参数是可取的。继续后面的操作。

2.3 数据标准化：使用LogNormalize方法，对每个细胞的总转录量进行归一化（总转录量为10000）

2.4 识别高变异基因：使用默认的 "vst" 方法选择2000个高变异基因，用于后续分析

标准化这一步有点花时间，大家注意一下：

运行完出图：

2.5 PCA分析

PCA分析的意思就是告诉我们这个数据受多少个维度的影响。出图如下，我们可以看到这些细胞主要受到10个PC维度的影响（10以后坡度变缓了，原文在选择维度的时候也是用的10）：

2.6 使用 Harmony 进行批次效应校正

2.7 使用 dim=10 的 Harmony 结果运行UMAP

2.8 FindNeighbors 和 FindClusters 确定最佳分辨率

OK，这里需要注意，我们用了10个PC，这是因为原文就是这个参数（不过一般也都是10），所以我们用了这个，要是全新的数据，你就得不停的尝试，看哪个PC维度最好。维度越大，考虑的因素就越多。不同的PC，降维之后的细胞图形是不同的。

同时，我们设置了0.01, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1分辨率。也就是我们确定了PC之后调整的参数，这个时候调整这个分辨率，图形是不变的，只是会分辨率变了而已，细胞的群的数量变了而已。如下：不同的分辨率有不同的细胞群数量，但是整个图形是不变的：

OK，这个时候我们一般选择0.8进行试探并后续分析。以上就是细胞的降维，告一段落。我们需要后续继续看0.8（0到28，有29群细胞）这个参数对不对，就是要细胞命名，看合理吗？

#3细胞命名

# A. 选择分辨率

# B. 绘制DotPlot，检查特定基因在各个细胞簇的表达情况

# C. 绘制TSNE图

# D. 绘制UMAP图

下面那些细胞的标志分子是参考的原文的，那我们分析新数据的时候也是参考同领域的作为参考，然后一个一个分子调整的，自己需要有知识储备。那我们可以看到有29群（0到28）和是否有批次效应：

（每一块细胞都有样本分布说明没有样本批次）

可以看到没有批次效应。因此，到这里我们可以判断去批次正确。我们要继续判断10个PC和0.8的分辨率是否正确。应该怎么做呢？就是最后命名好细胞名字后做标志基因的热图，看每个marker是否准确的落在每一群细胞中。

# E.根据已有marker判断哪个群是什么细胞:

这个时候找每一群细胞的高表达基因会花点时间，也可以保存这个scRNA_innitial_markers输出为表格，然后一个基因一个基因的看，再确定哪一群细胞是叫什么细胞，大家注意：

那我这里已知知道0.8的分辨率是29群细胞。同时，原文是用CD79A和JCHAIN两个分子作为B细胞的marker，那我们用FeaturePlot(scRNA,features = c("CD79A", "JCHAIN"), raster=FALSE)看看这两个分子落在哪里。发现是14和21群，那么我就会备注好，其他细胞也是一样的这样做好备注，方便后面的命名。其他细胞群也是一样。那到这一步，我们已经知道哪一群细胞叫什么名字了，我们就需要命名了。如下：

F. 细胞命名

命名后细胞降维图如下，跟原文非常接近：

我们这里检测上皮细胞1.4万，基质细胞1.3万：

原文也是：

G. 细胞簇的标志基因表达热图

那最后，我们要可视化各种marker看看，这些群的命名是否真的正确。原文长这样，各个marker非常清晰：

运行代码后：

简直和原文一模一样：

这下我们放心了，细胞分群对了，意味着10PC和0.8的分辨率是可以的，这些参数发论文的时候都要提供的。后面就简单了。

H. 保存和看看基因的表达

普通的图长这样：FeaturePlot(scRNA,features = c("CD79A"), raster=FALSE,split.by = "group")：

很丑，用Nebulosa包很美观，表达差异非常明显：

OK，以上全套代码已经给大家整理好啦，免费发送哦。

🔽①长按下方二维码关注🔽

②输入关键词：7787