单细胞亚群合理命名是数据分析基石啊

2024-06-08 08:27:15 浏览数 (2)

前面我们diss了一个单细胞数据挖掘文章，详见：天啦，啥肿瘤有如此高比例的树突细胞，研究者们就是把单核细胞当做是了树突细胞，闹了个笑话。强烈建议新手如果是想利用好单细胞公共数据集，可以看我们的数据分析顾问服务哈， 否则很容易被diss，参考：单细胞，单细胞，我要diss你！

其实这样的笑话在单细胞数据挖掘文章里面层出不穷，比如另外一个数据挖掘文章；《Identification of Five Hub Genes Based on Single-Cell RNA Sequencing Data and Network Pharmacology in Patients With Acute Myocardial Infarction》，链接是：https://www.frontiersin.org/articles/10.3389/fpubh.2022.894129/full ，里面的bug就更可怕了，居然是把成纤维细胞错误的命名成为了单核细胞。如下所示：

把成纤维细胞错误的命名成为了单核细胞

而且上面的成纤维其实是平滑肌细胞或者周细胞。这样的话，这个文章后面就没法看了。一般来说，我们拿到了肿瘤相关的单细胞转录组的表达量矩阵后的第一层次降维聚类分群通常是：

  • immune (CD45 ,PTPRC),
  • epithelial/cancer (EpCAM ,EPCAM),
  • stromal (CD10 ,MME,fibro or CD31 ,PECAM1,endo)

参考我前面介绍过 CNS图表复现08—肿瘤单细胞数据第一次分群通用规则，这3大单细胞亚群构成了肿瘤免疫微环境的复杂。绝大部分文章都是抓住免疫细胞亚群进行细分，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类作为第二次细分亚群。但是也有不少文章是抓住stromal 里面的 fibro 和endo进行细分，并且编造生物学故事的。

而且就算是前面的特异性基因不认识，后面每个单细胞亚群可以做功能注释，也可以看到单核细胞不应该是有成纤维细胞的标志性的ECM和focal adhesion功能啊，这两个功能明明是成纤维细胞亚群的特异性的 ：

ECM和focal adhesion功能

上面的注释，从代码角度是很简单的：

代码语言：javascript复制
head(marker_cosg)
 symbols_list <-  as.list(as.data.frame(apply(marker_cosg$names,2,head,100)))
 symbols_list
 source('../com_go_kegg_ReactomePA_human.R') 
com_go_kegg_ReactomePA_human(symbols_list, pro='test' )

也就是说， 上面的symbols_list变量是每个单细胞亚群的top基因，可以是top100或者1000，需要自己多次测试哈。我这里给大家看看示范数据，可以看到在r编程语言里面是一个简单的列表结构 哦 ：

代码语言：javascript复制
>  symbols_list <-  as.list(as.data.frame(apply(marker_cosg$names,2,head,10)))
>  symbols_list
$`CD14 Mono`
 [1] "S100A9"   "CCL2"     "CTSB"     "PLA2G7"   "S100A8"   "LGALS3"  
 [7] "C15orf48" "FCN1"     "CTSL"     "CD63"    

$pDC
 [1] "TSPAN13"       "MAP1A"         "SCT"           "MYBL2"        
 [5] "SMPD3"         "CLEC4C"        "DNASE1L3"      "LILRA4"       
 [9] "RP11-117D22.2" "PTGDS"        
 
$Mk
 [1] "PF4"    "PPBP"   "SDPR"   "GNG11"  "ACRBP"  "NRGN"   "TUBB1" 
 [8] "TREML1" "GP9"    "SPARC" 
   
....

然后只需要有一个com_go_kegg_ReactomePA_human函数就可以批量的注释这些单细胞亚群的top基因列表啦，而且是在go和kegg以及Reactome 数据库哦。

细胞外基质（Extracellular Matrix, ECM）和焦点粘附（Focal Adhesion）通常与成纤维细胞（Fibroblasts）的功能密切相关，并且被认为是成纤维单细胞亚群特异性的生物学功能，原因如下：

  1. ECM的合成与维护：成纤维细胞是细胞外基质的主要生产者。ECM是由多种大分子组成的复杂网络，包括胶原蛋白、弹性蛋白、纤维连接蛋白、层粘连蛋白以及其他多种糖蛋白和蛋白多糖。这些分子共同为组织提供结构支持、营养输送、信号传递等功能。成纤维细胞负责合成、组装、维护和重塑ECM，以适应组织的需要。
  2. 焦点粘附的形成：焦点粘附是细胞与ECM之间相互作用的关键结构，它们是细胞表面蛋白质（如整合素）与ECM组分之间的连接点。焦点粘附不仅有助于细胞锚定到ECM上，还传递机械信号和化学信号，影响细胞的行为，如迁移、增殖和分化。成纤维细胞通过形成焦点粘附来响应外部环境变化，参与组织修复和再生。
  3. 组织结构与功能：成纤维细胞在维持组织结构和功能中发挥核心作用。它们通过调节ECM的组成和刚性来控制组织的力学特性。此外，成纤维细胞还通过分泌生长因子和细胞因子参与组织的修复和免疫反应。

前面我们分享了在单细胞转录组降维聚类分群的第一层次降维聚类分群后的，每个单细胞亚群细分的时候，是有 单细胞亚群的生物学命名的4个规则，如下所示 ：

  • 第一个规则：已知的生物学亚群
  • 第二个规则：顺序编号加上特异性高表达量基因
  • 第三个规则：生物学功能注释
  • 第四个规则：转录因子等基因集特异性亚群（更多的生物学功能数据库）

一般来说，第一层次降维聚类分群应该是使用第一个规则，也可以辅助第三个规则，高通量数据分析很多证据链是环环相扣的！

批量注释多个基因列表在多个数据库代码

上面的 source('../com_go_kegg_ReactomePA_human.R') 文件里面的代码如下所示：

代码语言：javascript复制
com_go_kegg_ReactomePA_human <- function(symbols_list ,pro){ 
  library(clusterProfiler)
  library(org.Hs.eg.db)
  library(ReactomePA)
  library(ggplot2)
  library(stringr)
  
  # 首先全部的symbol 需要转为 entrezID
  gcSample = lapply(symbols_list, function(y){ 
    y=as.character(na.omit(AnnotationDbi::select(org.Hs.eg.db,
                                  keys = y,
                                  columns = 'ENTREZID',
                                  keytype = 'SYMBOL')[,2])
    )
    y
  })
  gcSample
  
  # 第1个注释是 KEGG 
  xx <- compareCluster(gcSample, fun="enrichKEGG",
                       organism="hsa", pvalueCutoff=0.05)
  dotplot(xx)    theme(axis.text.x=element_text(angle=45,hjust = 1))   
    scale_y_discrete(labels=function(x) str_wrap(x, width=50)) 
  ggsave(paste0(pro,'_kegg.pdf'),width = 10,height = 8)
  
  # 第2个注释是 ReactomePA 
  xx <- compareCluster(gcSample, fun="enrichPathway",
                       organism = "human",
                  pvalueCutoff=0.05)
  dotplot(xx)    theme(axis.text.x=element_text(angle=45,hjust = 1))   
    scale_y_discrete(labels=function(x) str_wrap(x, width=50)) 
  ggsave(paste0(pro,'_ReactomePA.pdf'),width = 10,height = 8)
  
  # 然后是GO数据库的BP,CC,MF的独立注释
  # Run full GO enrichment test for BP 
  formula_res <- compareCluster(
    gcSample, 
    fun="enrichGO", 
    OrgDb="org.Hs.eg.db",
    ont     = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff  = 0.05,
    qvalueCutoff  = 0.05
  )
  
  # Run GO enrichment test and merge terms 
  # that are close to each other to remove result redundancy
  lineage1_ego <- simplify(
    formula_res, 
    cutoff=0.5, 
    by="p.adjust", 
    select_fun=min
  ) 
  pdf(paste0(pro,'_GO_BP_cluster_simplified.pdf') ,width = 15,height = 8)
  print(dotplot(lineage1_ego, showCategory=5)   theme(axis.text.x=element_text(angle=45,hjust = 1))   
          scale_y_discrete(labels=function(x) str_wrap(x, width=50)) )
  dev.off() 
  write.csv(lineage1_ego@compareClusterResult, 
            file=paste0(pro,'_GO_BP_cluster_simplified.csv'))
  # Run full GO enrichment test for CC 
  formula_res <- compareCluster(
    gcSample, 
    fun="enrichGO", 
    OrgDb="org.Hs.eg.db",
    ont     = "CC",
    pAdjustMethod = "BH",
    pvalueCutoff  = 0.05,
    qvalueCutoff  = 0.05
  )
  
  # Run GO enrichment test and merge terms 
  # that are close to each other to remove result redundancy
  lineage1_ego <- simplify(
    formula_res, 
    cutoff=0.5, 
    by="p.adjust", 
    select_fun=min
  ) 
  pdf(paste0(pro,'_GO_CC_cluster_simplified.pdf') ,width = 15,height = 8)
  print(dotplot(lineage1_ego, showCategory=5)   theme(axis.text.x=element_text(angle=45,hjust = 1))   
          scale_y_discrete(labels=function(x) str_wrap(x, width=50)) )
  dev.off() 
  write.csv(lineage1_ego@compareClusterResult, 
            file=paste0(pro,'_GO_CC_cluster_simplified.csv'))
  
  # Run full GO enrichment test for MF 
  formula_res <- compareCluster(
    gcSample, 
    fun="enrichGO", 
    OrgDb="org.Hs.eg.db",
    ont     = "MF",
    pAdjustMethod = "BH",
    pvalueCutoff  = 0.05,
    qvalueCutoff  = 0.05
  )
  
  # Run GO enrichment test and merge terms 
  # that are close to each other to remove result redundancy
  lineage1_ego <- simplify(
    formula_res, 
    cutoff=0.5, 
    by="p.adjust", 
    select_fun=min
  ) 
  pdf(paste0(pro,'_GO_MF_cluster_simplified.pdf') ,width = 15,height = 8)
  print(dotplot(lineage1_ego, showCategory=5)   theme(axis.text.x=element_text(angle=45,hjust = 1))   
          scale_y_discrete(labels=function(x) str_wrap(x, width=50)) )
  dev.off() 
  write.csv(lineage1_ego@compareClusterResult, 
            file=paste0(pro,'_GO_MF_cluster_simplified.csv'))
  
}

上面的代码所利用的 GO、KEGG和Reactome是三个广泛使用的生物信息数据库，它们提供了关于生物分子功能、生物途径和疾病相关的重要信息。

Gene Ontology (GO) 数据库

GO数据库是一个国际性的协作项目，旨在为所有的生物体提供详细的、结构化的基因和蛋白质功能信息的标准化词汇。GO数据库包含三个主要的本体（Ontology）：

  1. 生物过程（Biological Process, BP）：描述基因和蛋白质参与的生物学事件的序列，例如细胞分裂或信号传导。
  2. 分子功能（Molecular Function, MF）：描述基因和蛋白质的生物学活性，例如酶的催化作用或受体的配体结合。
  3. 细胞组分（Cellular Component, CC）：描述基因和蛋白质在细胞中的位置，例如细胞核、细胞质或细胞膜。

GO数据库通过提供标准化的术语和定义，使得生物学数据的注释和比较变得更加容易。

Kyoto Encyclopedia of Genes and Genomes (KEGG) 数据库

KEGG是一个综合数据库资源，提供基因组、生物途径、疾病和化学物质等方面的信息。KEGG数据库的主要特点包括：

  • 途径数据库（KEGG PATHWAY）：提供了细胞和生物体内的代谢途径和信号传导途径的图形化表示。
  • 基因组数据库（KEGG GENOME）：提供了基因组数据和基因组内基因的功能注释。
  • 疾病数据库（KEGG DISEASE）：提供了人类疾病的信息，包括遗传病和复杂疾病的基因和途径信息。
  • 化学物质数据库（KEGG COMPOUND 和 KEGG DRUG）：提供了化学物质和药物的信息，以及它们在生物体内的代谢和作用。

KEGG数据库通过整合多种类型的数据，为系统生物学研究提供了宝贵的资源。

Reactome 数据库

Reactome是一个免费的、开放访问的数据库，专注于生物途径和反应网络。Reactome数据库的特点包括：

  • 详细的途径信息：提供了细胞内部发生的详细生物化学过程的描述，包括代谢、信号传导和细胞过程。
  • 数据集成：整合了来自多个数据库的数据，包括蛋白质、蛋白质复合物、小分子和RNA的信息。
  • 用户友好的界面：提供了一个易于使用的界面，允许用户浏览和查询途径信息。
  • 数据分析工具：提供了工具来分析高通量数据，例如通过基因集富集分析（GSEA）来研究基因表达数据。

Reactome数据库是研究生物学途径和进行系统生物学分析的重要资源。

这三个数据库都是生物医学研究中不可或缺的资源，它们支持研究人员在基因功能、途径分析和疾病机制研究方面的工作，并有助于解释实验数据和发现新的生物学知识。

写在文末

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

代码语言：javascript复制
We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。

数据 数据挖掘 go 数据库 数据分析

0 人点赞