查看原文
其他

什么时候P值大于0.05也无所谓呢

生信技能树 生信技能树 2022-06-07
号外:绝大部分生信技能树粉丝都没有机会加我微信,已经多次满了5000好友,所以我开通了一个微信小号,前100名添加我,仅需150元即可,3折优惠期机会不容错过哈。我的微信小号二维码在:0元,10小时教学视频直播《跟着百度李彦宏学习肿瘤基因组测序数据分析》

我在生信技能树分享了一个教程:不要怀疑,你的基因就是没办法富集到统计学显著的通路,然后在生信菜鸟团给了一个解决方案:差异基因没办法富集到通路你就放弃了吗,但是他们都是基于转录组表达量来的,这并不意味着GO/KEGG这样的生物学功能数据库的注释仅仅是只能针对转录表达的数据。

实际上,GO/KEGG数据库针对的是基因,不需要是表达量上下调的基因,也可以是突变与否的基因。比如下面的表格:

样本-基因突变与否二值矩阵

这是一个突变与否的二值矩阵,每一个病人的每一个基因,是否有突变,都列出来了。那么这个时候有一个分析就顺理成章了,就是我们通常是想关心具体的某些通路是否突变,而不是某一个基因,因为目前的药物领域还有很多是靶向通路的。上面的样本-基因突变与否二值矩阵,就需要进行GO/KEGG数据库注释,作为样本-通路突变与否二值矩阵。

这个时候就完全没有超几何分布检验的事情了,因为哪怕一个通路里面仅仅是只有一个基因发生了突变,我们都会认为这个通路是突变的,无需要计算P值。因为这个时候每个病人突变的基因数量本来就很少,通常是10个左右,这个数量级是无需超几何分布检验的。但是我们仍然是可以使用超几何分布检验,仅仅是需要把P值的阈值去掉即可:

genes=c("MLL2","BCL2","MYD88","HIST1H1E","PIM1",
        "CREBBP","CARD11","SPEN","TP53","ARID1A",
        "TNFRSF14","SOCS1","CDKN2A","NOTCH2","ARID1B",
        "SETD1B","GNA13","SMARCA4","SGK1","MGA","DUSP2","CD70",
        "ATM","BTG2","ZNF608","STAT6","KLHL6","MTOR","BIRC6","IRF8",
        "IRF4","PIK3CD","SETD2","TET2","B2M","TNFAIP3","ARID5B","EZH2",
        "ETS1","ZNF292","FAM5C","MCL1","EP300","MLL3","MEF2B","MYC","TBL1XR1",
        "DCAF6","DDX10","PAX5","KLHL14","GNAS","RB1","INO80","CHD8","BTG1","CD79B",
        "DNMT3A","FOXO1","PTPN6","BCL6","BCL7A","TMEM30A","ZEB2","JAK3","MECOM","EBF1",
        "ANKRD17","ETV6","ZFAT","NCOR1","IKZF3","MSH2","ATR","CHD1","NFKBIE","PRDM1",
        "CD83","STAT3","TAF1","ZBTB7A","CCND3","MAP4K4","NF1","UBE2A","DDX3X","IL16",
        "SETD5","JUNB","JAK1","PTPRK","CHST2","GOLGA5","NFKBIA","RHOA","ACTB","CDC73",
        "PIM2","BCL10","BTK","CD58","CBLB","MARK1","HNRNPU","SF3B1","UBR5","MYB",
        "TIPARP","NFKB2","MET","RRAGC","IKBKB","TOX","FAS","DICER1","TMSB4X","GNAI2",
        "MAGT1","PHF6","POU2F2","XPO1","WAC","YY1","CD22","HRAS","BRAF","FOXP1",
        "HIST1H2BC","RARA","S1PR2","FUBP1","STAT5B","KRAS","RUNX1","TCL1A","CXCR4",
        "PTEN","FBXW7","MAP2K1","ZFX","BTBD3","CIITA",
        "HNRNPD","LIN54","MSH6","PIK3R1","KCMF1","SYK","TGFBR2","CASP8")
patient=sample(genes,10)
library(clusterProfiler)
library(org.Hs.eg.db)
gene.df<-bitr(patient,fromType="SYMBOL",
              toType=c("ENSEMBL","ENTREZID"),
              OrgDb=org.Hs.eg.db)
head(gene.df)
genelist=gene.df$ENTREZID
kegg <- enrichKEGG(genelist, organism = 'hsa',
                   pvalueCutoff = 1,
                   qvalueCutoff =1)
kegg=setReadable(kegg,org.Hs.eg.db,keyType = 'ENTREZID')
head(kegg)
dt=kegg@result
dotplot(kegg,showCategory=30)

可以看到,随便去10个基因,都是可以在kegg里面找到113条通路的,虽然是,大量的通路仅仅是被“富集”到一个基因,如果从统计学角度来看,当然是没有“富集”的统计学意义啦。

“富集”的统计学结果

但是这113个通路里面,是有我们感兴趣的肿瘤免疫相关的,就可以制作  样本-通路突变与否的二值矩阵啦。

步骤如下:

  • 对每个病人的突变基因,都走上面的流程拿到全部的kegg通路
  • 然后把全部的病人的kegg通路结果数据塑性,长变宽
  • 如果要求更多的通路,就试试看msigdb啦

这个时候的超几何分布检验对我们来说,仅仅是起通路注释的作用,所以P值显著与否没有关系!

文末友情宣传

强烈建议你推荐我们生信技能树给身边的博士后以及年轻生物学PI,帮助他们多一点数据认知,让科研更上一个台阶:
推荐阅读







您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存