多种场合我都是高举生物信息学数据分析的靠谱的大旗，详见前面的推文：离大谱了，生信转湿实验？。仍然每次跟网友们的讨论的结局都是一地鸡毛，大家总会是把生物信息学数据分析跟常见的临床数据挖掘预后模型那些水文混为一谈。

最简单的一个生物信息学数据分析案例，可以是针对一个表达量矩阵以及里面的样品分组进行差异分析，然后拿到统计学显著的上下调基因列表，因为基因数量很多，所以我们默认会做一些生物学功能数据库注释。这里面的生物信息学算法就差异分析和注释算法， 它本身是非常可靠的！

数据分析有错是基本上不可能的，基本上都是新手才会翻车，比如表达量矩阵的格式，是否取多次log会影响变化倍数这个指标，以及上下调弄反的，而且这些错误很容易检验，不存在实验室和地理位置特异性。抛开这个不谈，但确实是有很多人对这个结果的解读是有误的， 比如几十个到几百个基因进行生物学功能数据库注释后有十几个通路，其中选择那个作为后续分析的指导确实是有主观层面的误差。

真正不靠谱的地方在于表达量矩阵产出这个实验阶段，你完全是没办法想象不同的实验室不同操作人员会引入多少误差。另外一个值得商榷的地方就是人类知识库的残缺，这个是泛生命科学领域的大bug，不能怪我们生物信息学数据分析。

比如前两天我们的马拉松授课交流群就有学员提到了他在复现一个中山大学的昼夜节律相关数据挖掘文章的时候出现了源头的背景知识冲突。文章是：《Circadian rhythm-related genes index: A predictor for HNSCC prognosis, immunotherapy efficacy, and chemosensitivity》，首先是一个简单的TCGA数据库的头颈癌的差异分析（502 HNSCC samples and 44 normal samples from the TCGA- HNSCC dataset），然后拿这个差异分析的上下调基因集合去跟circadian rhythm-related genes (CRRGs) 基因列表进行取交集操作。在文章的附件给出来了这个 SUPPLEMENTARY TABLE 1，是 1424 CRRGs were obtained through Genecards database.

昼夜节律相关基因列表到底是什么

文献里面的昼夜节律相关基因列表获取，方法学是：We obtained Circadian rhythm-associated 1424 genes by Genecards (https://www.genecards.org/) (Supplementary Table 1) ，我们下载该列表后简单的看看；

rm(list = ls()) 
CRGs_gene = rio::import('./Table_1_Circadian rhythm-related genes index A predictor for HNSCC prognosis, immunotherapy efficacy, and chemosensitivity.xlsx')
head(CRGs_gene)
CRGs_gene = unique(CRGs_gene[,1])
save(CRGs_gene,file = 'gene_from_genecards.Rdata')

从Genecards数据库获取昼夜节律相关基因列表，就是一个简单的关键词搜索即可，然后基因会按照搜索关键词的相关性排序，卡一个阈值即可：

> head(CRGs_gene) 
 Gene Symbol Description
1 CLOCK Clock Circadian Regulator
2 PER2 Period Circadian Regulator 2
3 PER3 Period Circadian Regulator 3
4 PER1 Period Circadian Regulator 1
5 CRY1 Cryptochrome Circadian Regulator 1
6 ARNTL Aryl Hydrocarbon Receptor Nuclear Translocator Like
> dim(CRGs_gene)
[1] 1424 7

但是如果大家简单的搜索，就可以看到华中科技大学的一个在线数据库“《Circadian Gene DataBase》”，发表在Nucleic Acids Res. 2017:，CGDB: a database of circadian genes in eukaryotes.，可以很简单的从这个课题组的整理好的数据库里面拿到昼夜节律相关基因列表 ，因为网页里面的人类的昼夜节律相关基因列表是以蛋白质序列方式给出来的，所以需要简单的转换 ：

rm(list = ls()) 
# # download.file(url = 'http://cgdb.biocuckoo.org/download/all/Homo%20sapiens%20(Human).txt',
# # destfile = 'CRGs_gene.txt')
# 
CRGs_gene = readLines('CRGs_gene.txt')
CRGs_gene = CRGs_gene[grepl('>',CRGs_gene)]
CRGs_gene = stringr::str_split(CRGs_gene,'\\|',simplify = T)[,2]
head(CRGs_gene)
tail(CRGs_gene)

ens_crgs = CRGs_gene[grepl('^ENSP',CRGs_gene)]
uniport_crge = CRGs_gene[!CRGs_gene %in% ens_crgs]
ens_crgs = stringr::str_split(ens_crgs,'\\.',simplify = T)[,1]

#gene转换为symbol
library(clusterProfiler)
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db)
ens_crgs2 = bitr(ens_crgs,
 fromType = 'ENSEMBLPROT',
 toType = 'SYMBOL',
 OrgDb = 'org.Hs.eg.db')
# 76.63% of input gene IDs are fail to map...
uniport_crge2 = bitr(uniport_crge,
 fromType = 'UNIPROT',
 toType = 'SYMBOL',
 OrgDb = 'org.Hs.eg.db')
# 14.03% of input gene IDs are fail to map...

CRGs_gene = unique(c(ens_crgs2$SYMBOL,uniport_crge2$SYMBOL))
head(CRGs_gene)
save(CRGs_gene,file = 'gene_from_cgdb.Rdata')

接下来就可以看看这两个来源的昼夜节律相关基因列表的交集如何了：

rm(list = ls()) 
load(file = 'gene_from_genecards.Rdata')
genecards=CRGs_gene
load(file = 'gene_from_cgdb.Rdata')
cgdb=CRGs_gene 
genes<- list(genecards = genecards,
 cgdb = cgdb ) #输入数据需要是
#设定颜色：长度要与list内容长度一致（这里提供备选了几种）
colors<- c("#2874C5", "#f87669", "#e6b707", "#868686", "#92C5DE", "#F4A582", "#66C2A5", "#FC8D62", "#8DA0CB", 
 "#E78AC3", "#A6D854", "#FFD92F", "#E5C494", "#B3B3B3") [1:length(lengths(genes))] 
p = VennDiagram::venn.diagram(x = genes,
 filename = NULL, #不生成文件（还是必写参数，烦）
 fill = colors, #设定颜色(长度要与分组数相等)
 cat.col = colors, cat.cex = 1.5, #标签颜色，字体大小
 col = colors, #边框颜色
 alpha = 0.1 #调整透明度
) 
grid::grid.draw(p)
library(ggplot2)
#转换为ggplot对象
p<- ggplotify::as.ggplot(cowplot::as_grob(p))
p
ggsave(p,file = './Venn.jpg')

int=intersect(genecards,cgdb)
symbols_list = list(
 int=int,
 only_genecards=setdiff(genecards,int),
 only_cgdb=setdiff(cgdb,int)
)
lapply(symbols_list, length)

可以看到，两者之间的交集居然是就各自的十分之一不到：

$int
[1] 166

$only_genecards
[1] 1258

$only_cgdb
[1] 1105

> lapply(symbols_list, head)
$int
[1] "PER2" "PER3" "PER1" "CRY1" "CRY2" "NPAS2"

$only_genecards
[1] "CLOCK" "ARNTL" "TIMELESS" "ARNTL2" "CIPC" "CSNK1D"

$only_cgdb
[1] "ZNF251" "C1orf54" "MGAM" "KIR2DL5B" "KIR2DS1" "H2BC6"

在各个癌症里面的做昼夜节律相关基因相关预后模型的数据挖掘文章非常多，做乳酸化，线粒体，衰老相关基因列表就更多了。如果大家根本就不是从同一个基因列表开始，那凭什么保证这些数据挖掘文章的可靠性和合理性呢？

这些基因简单的kegg数据库 注释：

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
pro='test'
# 第1个注释是 KEGG 
xx <- compareCluster(gcSample, fun="enrichKEGG",
 organism="hsa", pvalueCutoff=0.3)
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)

也很难下结论，感觉就是不同来源的数据库确实是需要有benchmark了，这个也是大量的生物信息学课题了：

比如大家如果是想在各个癌症里面的做转录因子相关基因相关预后模型的数据挖掘文章，就需要明确人类的两万多个蛋白质编码基因里面到底是哪些基因是转录因子呢？大概有1500到2000个转录因子基因，但是具体是什么，又是有不同来源的数据库了。比如在生信菜鸟团公众号看到了有一个介绍：TCGA数据挖掘常见基因集合，首先是Cancer Manag Res. 2020的文章《Prognostic and Predictive Value of a 15 Transcription Factors (TFs) Panel for Hepatocellular Carcinoma》就提到了：

• The TF list was downloaded from the Human Transcription Factors website (http://humantfs.ccbr.utoronto.ca/.

然后是2021的文章《A Transcription Factor-Based Risk Model for Predicting the Prognosis of Prostate Cancer and Potential Therapeutic Drugs》提到一个出处：

• Atotal of 1665 transcription factors were obtained from the Animal TFDB database。（AnimalTFDB 3.0 ），链接：http://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/

这两个数据库关于转录因子的收录，都是接近于2000个基因。然后是刘小乐的Cistrome数据库，详见：http://cistrome.org/db/#/stat

我下载里面的 human_factor_full_QC.txt 文件，然后统计了一下，在人类这个研究领域，有chip-seq数据的转录因子是1359个，略低于上面的两个网页数据库里面的1600~2000的数量。

还有 RcisTarget包里面的motifAnnotations_hgnc数据，也有自己的转录因子基因合辑。

如果做各个癌症里面的做转录因子相关基因相关预后模型的数据挖掘文章之前压根就没有了解过什么是转录因子，没有读下面的这两个综述文章：

• 首先是2018的CELL文章：《The Human Transcription Factors》
• 然后是2020的NBT文章：《A comprehensive library of human transcription factors for cell fate engineering》

做出来的数据挖掘文章又价值几何呢？

不同算法的一致性呢

数据源头的质量问题，以及人类知识数据库的不完善，都是会影响生物信息学数据分析结果的，但是这些是目前很难完全责怪生物信息学本身，反而是算法的多样性往往是生物信息学被人诟病的一个点了。比如之前我们介绍过irGSEA：基于秩次的单细胞基因集富集分析整合框架，针对17种常见的Functional Class Scoring (FCS)方法进行了benchmark。

单细胞转录组数据分析最基础的就是给每个细胞一个身份，通常是降维聚类分群后然后对每个亚群进行描述，首先可以描述每个亚群的高表达量的特异性基因，然后可以对基因进行生物学功能数据库注释。也可以直接对亚群进行打分，比如gsea和gsva算法，针对免疫或者代谢等基因集进行打分，也是拿到每个亚群的特异性生物学功能。如果选择了不同打分算法，往往是拿到了不同结果，有时候确实是很难解释。

未解之谜很正常

很多人对生物信息学数据处理的另外一个误解是，它应该是有百分百的确定性的的分析结果。其实也不是，比如不同单细胞转录组数据集的降维聚类分群其实都会有 热激蛋白的亚群，细胞增殖亚群，干扰素亚群，金属离子酶亚群，线粒体或者核糖体亚群，或者低质量亚群。。。但是每个数据分析人员在处理这些细胞状态的变化的时候都有自己的人为干预。

当然了，这个其实又是归咎到了人类知识数据库的不完善，本来就暂时还没有定论，关于细胞类型和细胞状态的讨论一直在持续：

• 1.通俗地说，细胞类型是体内执行特定功能的一类细胞。通常细胞类型被认为是相对稳定的，就像皮肤中的细胞不会自发地变成神经元。
• 2.我们现在知道细胞实际上具有很强的可塑性，虽然皮肤细胞自然不会变为神经元，但可以通过特殊处理诱导它们变形为神经元。
• 3.此外细胞确实相对频繁地切换其功能，淋巴结中的T细胞可以“激活”以对抗感染，我们可以称这种瞬态细胞状态为新的“细胞类型”吗？
• 4.最后，还有如何处理患病细胞的问题，不再能够执行神经元功能的神经元仍然是神经元吗？“患病”神经元有自己的细胞类型定义吗？我们用什么标准来确定细胞是否“患病”？

感兴趣的可以仔细品读：https://mbernste.github.io/posts/cell_types_cell_states/