前些天我们众筹了这个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这样的多个单细胞样品的数据整合过程，就是右图，每个单细胞亚群比如t细胞和b细胞都是可以很明显的清晰的独立成群，但是具体到每个群里面又确实是可以看到样品之间的差异。而且呢，上皮细胞和成纤维又确实是很明显的有个体差异，就好奇怪，因为上皮细胞通常情况下是可以解释为恶性上皮细胞的肿瘤病人之间的异质性所以有个体差异，成纤维我就暂时解释不清楚。
然后呢我们的左边的图，很明显是harmony之后的效果，无论是什么样的单细胞亚群都很清晰的独自成群了，而且每个群内部的每个单细胞样品个体也很好的很均匀的混合起来了。说明harmony确实是抹去了单细胞样品本身差异，但是同时也抹去了肿瘤病人的癌细胞的异质性，这个就需要单独去讨论了。
但是，这个并不是重点，因为多个单细胞样品的数据整合与否，都有可取之处，真正的麻烦的地方是细胞数量很难使用作者在文章里面的公布的质量控制标准来达到，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%。这种方法适用于非正态分布的数据，或者当你认为数据的极端值是由特定的生物学过程或技术因素引起的。
选择哪种方法取决于以下几个因素：
数据分布：如果数据近似正态分布，使用2SD或3SD可能更合适。如果数据分布偏斜或有多个峰值，可能需要考虑其他方法。
研究目的：如果你的目标是去除由于技术或生物学原因产生的异常值，可能需要更具体的标准。
领域惯例：不同的研究领域可能有不同的标准和惯例。查看相关文献或与领域专家交流可以帮助确定最合适的方法。
数据敏感性：在某些情况下，过度过滤可能会移除重要的生物学信号。因此，需要权衡过滤的严格性与数据完整性。
数据量：样本量的大小也会影响过滤策略。在样本量较小的情况下，过滤可能会对结果产生较大影响。
在实际操作中，你可能需要尝试不同的方法，并结合数据的可视化（如箱线图、散点图）来评估哪种方法最适合你的数据。此外，过滤前后的数据都应该进行仔细检查，以确保过滤过程没有引入新的偏差。
因为作用最大的就是 nCount_RNA < fivenum(sce.all.filt$nCount_RNA)[4] ，虽然我不知道它的合理性，但确实是达到了文章类似的细胞数量，从36.6万多的细胞数量过滤到227,836 个单细胞，如下所示：
如果你也想做单细胞转录组数据分析，最好是有自己的计算机资源哦，比如我们的2024的共享服务器交个朋友福利价仍然是800，而且还需要有基本的生物信息学基础，也可以看看我们的生物信息学马拉松授课（买一得五） ，你的生物信息学入门课。而且这个周六日我们在长沙线下授课哦：千呼万唤，让我们长沙线下约起
