基因的功能推断之大队列表达量相关性排序后gsea分析

前面我们强调了，基因功能推断的数据分析的重要性 ，而且我们已经演示了第一个基因的功能推断方法就是：大队列高低分组后差异分析然后功能富集，但实际情况下，我们的表达量并不是离散变量，它高低分组是我们认为划分的，未必不可以高中低的3分组，或者更多可能性，比如根据它的表达量的密度分布曲线来看最佳的阈值进行分组。实际上既然表达量是连续值，我们就可以直接使用连续变量的统计学方法，连续型变量就可以看相关性。

接下来我们就演示大队列表达量相关性排序后gsea分析，这个时候以2023的这个文章为例：《Comprehensive analysis reveals CCDC60 as a potential biomarker correlated with prognosis and immune infiltration of head and neck squamous cell carcinoma》，研究者们首先说明了Coiled-coil domain containing 60 (CCDC60) 这个基因在头颈癌是一个典型的抑癌基因：

• 首先呢，它在癌症样品里面的表达量所显著的低于正常的癌旁组织
• 然后呢，它表达量越高的病人的预后反而越好

很明显，Coiled-coil domain containing 60 (CCDC60) 这个基因也不是什么明星基因，那么就需要进行基因的功能推断。让我们一起来看看作者所如何做的吧，方法学是：

可能是研究者们不怎么懂代码，所以使用了LinkedOmics这个数据库，实际上自己很容易下载tcga数据库里面的头颈癌队列里面的转录组表达量矩阵，然后计算目标基因跟所有的其它基因的相关性，然后根据相关性进行排序，就能定位到跟目标基因的表达量最相关的那些基因：

有了相关性排序列表，同样的可以进行gsea分析或超几何分布检验，针对go或者kegg等生物学功能数据库。下面是一个示例代码：

load(  file = 'step1-output.Rdata')
dat[1:4,1:4] 
head(phe)  
this_gene = 'CCDC60'
cc = apply(dat, 1, function(x) cor( dat[this_gene,],x))
df = as.data.frame(cc) 
head(df)
df$SYMBOL = rownames(df)
df2 <- bitr(unique(df$SYMBOL), fromType = "SYMBOL",
            toType = c( "ENTREZID"),
            OrgDb = org.Hs.eg.db) 
DEG=merge(df2,df,by='SYMBOL')
colnames(DEG)
DEG =na.omit(DEG)
data_all_sort <- DEG %>%  #排序
  arrange(desc(cc))
geneList = data_all_sort$cc #把foldchange按照从大到小提取出来
names(geneList) <- data_all_sort$ENTREZID #给上面提取的foldchange加上对应上ENTREZID
head(geneList) #排序好的基因序列，而且是entrezeID的形式
R.utils::setOption( "clusterProfiler.download.method",'auto' )
kkgsea <- gseKEGG(geneList     = geneList,
                  organism     = 'hsa', 
                  minGSSize    = 10,
                  maxGSSize    = 500,
                  pvalueCutoff = 1,
                  pAdjustMethod = "none" ) #进行gseKEGG富集分析