前些天我们众筹了这个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细胞圆润的有点让人吃惊。。。
这个单细胞转录组数据集的降维聚类分群结果显示了多个不同的细胞亚群,每个亚群由特定的标记基因定义。从文章里面可以很清晰的看到各个单细胞亚群对应的标记基因:
这些细胞负责产生细胞外基质,通常在组织结构和修复中起重要作用。
T细胞是适应性免疫系统的关键组成部分,参与细胞免疫反应。
标记基因:CD14、MS4A7(可能指CD20)、FCGR3A(可能指Fcγ受体)
髓系细胞包括多种类型的免疫细胞,如巨噬细胞、树突细胞和粒细胞。
上皮细胞构成组织的表面,如皮肤、消化道和呼吸道的内衬。
标记基因:MZB1(可能指XBP1)、IGKC、DERL3
B细胞是适应性免疫系统的另一个关键组成部分,负责产生抗体和提供抗原呈递。
内皮细胞(Endothelial Cells)
:
标记基因:TPSAB1(可能指Tryptase beta 1)、CPA3(可能指Carboxypeptidase A3)、MS4A2(可能指CD20)
肥大细胞在过敏反应和免疫防御中起作用,能够释放多种炎症介质。
这些细胞亚群的识别有助于理解组织内的细胞多样性和它们在健康和疾病中的作用。通过分析这些亚群的基因表达模式,研究人员可以探索不同细胞类型在生物学过程和疾病状态中的特定功能。
因为文章里面的上皮细胞和间质细胞是比较散乱的,所以就很难确定文章是否走了harmony类似的多个单细胞样品的数据整合过程,那我们自己读取一下cellranger出来的表达量矩阵然后降维聚类分群看看是否走harmony的不同选择可以对比看看:
seuratObj reduction = "pca" , dims = 1 :15 ) colnames([email protected] ) 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,则该细胞可能被认为是低质量的,可能会被移除。
这个阈值用于确保细胞具有足够的基因表达信息。如果一个细胞表达的基因数量少于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 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]
,但是这个策略的统计学理论很难站住脚。因为在数据处理中,过滤极端值(也称为离群值)的方法取决于数据的分布、研究目的和领域惯例。以下是几种常见的方法:
这种方法基于正态分布的假设,认为数据中约95%的值会落在均值的2个标准差之内。超出这个范围的值可能被视为离群值。