重要的并不是整合与否，而应该是质量控制

生信技能树

发布于 2024-11-21 10:18:13

1.7K0

前些天我们众筹了这个2024-HRA004702-胃癌-卵巢转移-单细胞数据集：45 tumor samples from 18 GC patients with OM.

是2024年9月30日，浙江省肿瘤医院程向东教授团队联合北京大学生物医学前沿创新中心白凡教授团队合作在Nature Communications上发表了题为The estrogen response in fibroblasts promotes ovarian metastases of gastric cancer的研究成果，该研究发现雌激素调节卵巢成纤维细胞会促进胃癌的卵巢转移，可能是胃癌卵巢转移好发于年轻女性的重要原因。对应的数据集链接是：https://ngdc.cncb.ac.cn/gsa-human/browse/HRA004702

标题：Application of single cell sequencing in metastatic gastric cancer
发布日期：2024-07-14
描述信息：Single cell sequencing was performed on samples of ovarian metastasis of gastric cancer to construct a single cell panorama of ovarian metastasis of gastric cancer
访问方式：   Open access
BioProject编号:PRJCA017230

因为是公开了fastq文件格式的测序数据，大概耗费了7个T的磁盘空间，一周的下载数据和计算时间，所以很容易走cellranger的定量流程拿到表达量矩阵。但是，拿到了表达量矩阵先小伙伴们反馈说压根就拿不到文章那样的降维聚类分群结果，如下所示：

文章那样的降维聚类分群结果

从上面的UMAP图可以看到，作者得到了227,836 个单细胞，然后第一层次降维聚类分群后的上皮细胞和间质细胞是比较散乱的，其它单细胞亚群就可以比较好的把来源于不同病人不同部位的单细胞样品聚集在一起。而且里面的t细胞圆润的有点让人吃惊。。。

这个单细胞转录组数据集的降维聚类分群结果显示了多个不同的细胞亚群，每个亚群由特定的标记基因定义。从文章里面可以很清晰的看到各个单细胞亚群对应的标记基因：

基质细胞（Stromal Cells）：
- 标记基因：COL1A1、COL1A2、DCN
- 这些细胞负责产生细胞外基质，通常在组织结构和修复中起重要作用。
T细胞（T Cells）：
- 标记基因：CD3E、CD3D、CD3G
- T细胞是适应性免疫系统的关键组成部分，参与细胞免疫反应。
髓系细胞（Myeloid Cells）：
- 标记基因：CD14、MS4A7（可能指CD20）、FCGR3A（可能指Fcγ受体）
- 髓系细胞包括多种类型的免疫细胞，如巨噬细胞、树突细胞和粒细胞。
上皮细胞（Epithelial Cells）：
- 标记基因：EPCAM、KRT8、KRT19
- 上皮细胞构成组织的表面，如皮肤、消化道和呼吸道的内衬。
浆细胞（Plasma Cells）：
- 标记基因：MZB1（可能指XBP1）、IGKC、DERL3
- 浆细胞是B细胞的终末分化形式，负责产生抗体。
B细胞（B Cells）：
- 标记基因：CD79A、MS4A1、CD19
- B细胞是适应性免疫系统的另一个关键组成部分，负责产生抗体和提供抗原呈递。
内皮细胞（Endothelial Cells）：
- 标记基因：VWF、PECAM1、ENG
- 内皮细胞构成血管的内衬，参与血液和淋巴循环。
肥大细胞（Mast Cells）：
- 标记基因：TPSAB1（可能指Tryptase beta 1）、CPA3（可能指Carboxypeptidase A3）、MS4A2（可能指CD20）
- 肥大细胞在过敏反应和免疫防御中起作用，能够释放多种炎症介质。

这些细胞亚群的识别有助于理解组织内的细胞多样性和它们在健康和疾病中的作用。通过分析这些亚群的基因表达模式，研究人员可以探索不同细胞类型在生物学过程和疾病状态中的特定功能。

因为文章里面的上皮细胞和间质细胞是比较散乱的，所以就很难确定文章是否走了harmony类似的多个单细胞样品的数据整合过程，那我们自己读取一下cellranger出来的表达量矩阵然后降维聚类分群看看是否走harmony的不同选择可以对比看看：

seuratObj <- RunUMAP(sce.all.int,  
                     reduction = "pca",
                     dims = 1:15 )
colnames(sce.all.int@meta.data)
p3=DimPlot(sce.all.int, reduction = "umap", 
        group.by = "celltype" ) #+ NoLegend()
p4=DimPlot(seuratObj, reduction = "umap" , 
          group.by = "celltype")+ NoLegend()

p1=DimPlot(sce.all.int, reduction = "umap", 
           group.by = "orig.ident" )+ NoLegend()
p2=DimPlot(seuratObj, reduction = "umap" , 
           group.by = "orig.ident")+ NoLegend()

(p1+p2)/(p3+p4)
ggsave('harmony-or-not.pdf',width = 12,height = 12)

详细的代码在百度云网盘链接: https://pan.baidu.com/s/1QRFWje5tI6Nodw3I3EX5Tg?pwd=7xp4 提取码: 7xp4 ，这里我们仅仅是展示harmony与否的UMAP的差异：

harmony与否的UMAP的差异

很明显的是，如果不走harmony这样的多个单细胞样品的数据整合过程，就是右图，每个单细胞亚群比如t细胞和b细胞都是可以很明显的清晰的独立成群，但是具体到每个群里面又确实是可以看到样品之间的差异。而且呢，上皮细胞和成纤维又确实是很明显的有个体差异，就好奇怪，因为上皮细胞通常情况下是可以解释为恶性上皮细胞的肿瘤病人之间的异质性所以有个体差异，成纤维我就暂时解释不清楚。

然后呢我们的左边的图，很明显是harmony之后的效果，无论是什么样的单细胞亚群都很清晰的独自成群了，而且每个群内部的每个单细胞样品个体也很好的很均匀的混合起来了。说明harmony确实是抹去了单细胞样品本身差异，但是同时也抹去了肿瘤病人的癌细胞的异质性，这个就需要单独去讨论了。

但是，这个并不是重点，因为多个单细胞样品的数据整合与否，都有可取之处，真正的麻烦的地方是细胞数量很难使用作者在文章里面的公布的质量控制标准来达到，36.6万多的细胞数量过滤到227,836 个单细胞：

36.6万多的细胞数量过滤到227,836 个单细胞

从36.6万多的细胞数量过滤到227,836 个单细胞，确实是很严格的质量控制了，但是上面写的 maximum percentage mito = 10%, minimum number of UMIs = 1000, and minimum nGene = 500. ：

最大线粒体基因表达百分比（Maximum percentage mito）：
- 这个阈值用于识别可能受到线粒体DNA污染或具有高线粒体基因表达的细胞，这可能是由于细胞应激或损伤。设置为10%意味着如果一个细胞的线粒体基因表达量占其总表达量的10%以上，则该细胞可能会被移除或标记为异常。
最小UMI数量（Minimum number of UMIs）：
- UMI（Unique Molecular Identifiers）是用于标记单个mRNA分子的分子条形码，有助于准确估计基因表达水平。最小UMI数量为1000意味着如果一个细胞的UMI计数少于1000，则该细胞可能被认为是低质量的，可能会被移除。
最小基因数量（Minimum nGene）：
- 这个阈值用于确保细胞具有足够的基因表达信息。如果一个细胞表达的基因数量少于500，则该细胞可能被认为是低质量的，可能会被移除。

这些过滤阈值有助于从数据集中移除可能影响分析结果的低质量细胞，例如那些受到污染、死亡或处于非典型状态的细胞。通过应用这些阈值，研究人员可以提高后续分析的可靠性和生物学解释的有效性。需要注意的是，这些阈值可能需要根据具体的实验设计、数据质量和研究目的进行调整。在设置这些阈值时，应考虑数据的分布和生物学背景。

初始矩阵读入后，可以看到其实每个样品细胞数量仍然是在5到10k之间，问题不大，基本上没什么双细胞的可能性。

> tmp
        [,1]  [,2]
P01_OM 26591  9653
P01_PM 23848  6969
P01_PT 23078  2695
P02_AS 24022 10104
P02_OM 26501  7984
P02_PM 25401  8293
P02_PT 22542  3728
P03_OM 25249  7277
P03_PM 25308  9511
P03_PT 23062  7037
P04_OM 23706  5830
P04_PM 23691  5598
P04_PT 24104  6626
P05_AS 24756  7674
P05_OM 24072  5462
P05_PT 24004  4532
P06_AS 17451  3841
P06_PT 17385  1685
P07_AS 22084  5121
P07_OM 18518  2797
P09_AD 19461  4118
P09_AS 19201 10919
P09_PT 20720  3269
P10_AS 19973  7971
P10_PT 20715  4473
P11_OM 27215  7725
P11_PM 27581  9678
P11_PT 26387  7011
P12_AD 21277  4832
P13_AD 20859  4064
P14_AD 21774  5066
P15_AD 22864  7583
P17_AS 25084  9346
P17_OM 27155  8546
P17_PT 24893 15448
P18_AS 24528 11413
P18_OM 26289 11252
P18_PT 26054 13154
P19_OM 26288 11826
P19_PT 25159 10617
P20_OM 25568  6676
P20_PT 22662  3607
P21_OM 25165 11636
P21_PT 22638  3912
P22_OM 26131  8315
P22_PT 25450  7028
P23_OM 27676 10125
P23_PT 25012  6443
P24_OM 26163  9735
P24_PT 26606  8118

如果大家使用同样的阈值，其实压根就没办法过滤到文献一样的结果，唯一剩下的就是在每个样品里面的使用 R package scDblFinder (v.1.4.0) ，我测试了一下，如果想达到文献里面的过滤程度是需要卡上限而不是下限，比如 nFeature_RNA 和nCount_RNA的上限的指标：

 
sce.all.filt
fivenum(sce.all.filt$nCount_RNA)
plot(sort(sce.all.filt$nCount_RNA))
tmp=subset(sce.all.filt,subset = percent_mito < 10 &
             nCount_RNA < fivenum(sce.all.filt$nCount_RNA)[4]  &
             nFeature_RNA <7500&
             nCount_RNA > 1000 & 
             nFeature_RNA>500)
tmp

sce.all=tmp
sort(table(tmp$orig.ident))
sce.all.filt =tmp

其实，上面的质量控制指标里面，作用最大的就是 nCount_RNA < fivenum(sce.all.filt$nCount_RNA)[4] ，但是这个策略的统计学理论很难站住脚。因为在数据处理中，过滤极端值（也称为离群值）的方法取决于数据的分布、研究目的和领域惯例。以下是几种常见的方法：

Mean + 2SD（均值加2个标准差）：
- 这种方法基于正态分布的假设，认为数据中约95%的值会落在均值的2个标准差之内。超出这个范围的值可能被视为离群值。
Mean + 3SD（均值加3个标准差）：
- 类似于2SD方法，但更为宽松，允许更多的数据点被保留。在正态分布中，约99.7%的值会落在均值的3个标准差之内。
去除Top 25% ：
- 这种方法不依赖于标准差，而是直接去除数据中最高的25%。这种方法适用于非正态分布的数据，或者当你认为数据的极端值是由特定的生物学过程或技术因素引起的。

选择哪种方法取决于以下几个因素：