1.已经确定研究的基因，但是想探索他潜在的功能，可以通过跟这个基因表达最相关的基因来反推他的功能，这种方法在英语中称为guilt of association，协同犯罪。
    2.我们的注释方法依赖于TCGA大样本，既然他可以注释基因，那么任何跟肿瘤相关的基因都可以被注释，包括长链非编码RNA。

    这个方法以前阐述过：
    单基因批量相关性分析的妙用
    但是这个方法有个小缺陷，并不知道最后富集的通路是正向影响还是反向影响，也就是无法判断方向。判断方向的工具也不是没有，GSEA就是一个。所以，我想能不能把批量相关性分析和GSEA结合起来。
    GSEA需要的gene set是现成的没有问题，但是genelist没有，这里我们可以把所有基因跟单个基因的相关性系数当做LogFC，有正有负，就解决了geneList的问题。这个想法不是我的，是我的一个学员的，不过他要解决的是microRNA把基因的问题。

    下面来实战一下：
    1.首先加载数据
    这个数据是我下载了TPM数据，然后提取出乳腺癌的数据得来的。

    1. load(file = "BRCA_mRNA_exprSet.Rdata")
    3. exprSet <- mRNA_exprSet
    4. test <- exprSet[1:10,1:10]

    image.png
    2.写一个函数批量计算相关性
    这个函数只要输入一个基因，他就会批量计算这个基因跟其他编码基因的相关性，返回相关性系数和p值。

    1. batch_cor <- function(gene){
    2. y <- as.numeric(exprSet[gene,])
    3. rownames <- rownames(exprSet)
    4. do.call(rbind,future_lapply(rownames, function(x){
    5. dd  <- cor.test(as.numeric(exprSet[x,]),y,type="spearman")
    6. data.frame(gene=gene,mRNAs=x,cor=dd$estimate,p.value=dd$p.value )
    7. }))
    8. }

    3.并行化运行函数
    以PCDC1这个基因为例

    1. library(future.apply)
    2. plan(multiprocess)
    3. system.time(dd <- batch_cor("PDCD1"))

    这是返回的结果
    image.png
    4.制作genelist gene

    1. gene <- dd$mRNAs
    3. ## 转换
    4. library(clusterProfiler)
    5. gene = bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
    6. ## 去重
    7. gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)
    9. gene_df <- data.frame(logFC=dd$cor,
    10. SYMBOL = dd$mRNAs)
    11. gene_df <- merge(gene_df,gene,by="SYMBOL")
    13. ## geneList 三部曲
    14. ## 1.获取基因logFC
    15. geneList <- gene_df$logFC
    16. ## 2.命名
    17. names(geneList) = gene_df$ENTREZID
    18. ## 3.排序很重要
    19. geneList = sort(geneList, decreasing = TRUE)

    5.运行GSEA分析

    1. library(clusterProfiler)
    3. ## 读入hallmarks gene set，从哪来？
    4. hallmarks <- read.gmt("h.all.v6.2.entrez.gmt")
    5. # 需要网络
    6. y <- GSEA(geneList,TERM2GENE =hallmarks)

    作图看整体分布

    1. ### 看整体分布
    2. library(ggplot2)
    3. dotplot(y,showCategory=12,split=".sign") facet_grid(~.sign)

    本次结果中全是激活的
    image.png
    6.特定通路作图

    1. yd <- data.frame(y)
    2. library(enrichplot)
    3. gseaplot2(y,"HALLMARK_INTERFERON_ALPHA_RESPONSE",color = "red",pvalue_table = T)

    PCDC1跟阿拉法干扰素正相关，这个事情没什么好说的吧。