基于单基因批量相关性分析的GSEA—科研工具箱

有这样的使用场景么？

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

下面来实战一下：

1.首先加载数据

这个数据是我下载了TPM数据，然后提取出乳腺癌的数据得来的。

load(file = \'BRCA_mRNA_exprSet.Rdata\')

exprSet <- mRNA_exprSet

test <- exprSet[1:10,1:10]

2.写一个函数批量计算相关性

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

batch_cor <- function(gene){
  y <- as.numeric(exprSet[gene,])
  rownames <- rownames(exprSet)
  do.call(rbind,future_lapply(rownames, function(x){
    dd  <- cor.test(as.numeric(exprSet[x,]),y,type=\'spearman\')
    data.frame(gene=gene,mRNAs=x,cor=dd$estimate,p.value=dd$p.value )
  }))
}

3.并行化运行函数

以PCDC1这个基因为例

library(future.apply)
plan(multiprocess)
system.time(dd <- batch_cor(\'PDCD1\'))

这是返回的结果

4.制作genelist

gene <- dd$mRNAs## 转换
library(clusterProfiler)
gene = bitr(gene, fromType=\'SYMBOL\', toType=\'ENTREZID\', OrgDb=\'org.Hs.eg.db\')## 去重
gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)
gene_df <- data.frame(logFC=dd$cor,
                      SYMBOL = dd$mRNAs)
gene_df <- merge(gene_df,gene,by=\'SYMBOL\')
## geneList 三部曲
## 1.获取基因logFC
geneList <- gene_df$logFC
## 2.命名
names(geneList) = gene_df$ENTREZID
## 3.排序很重要
geneList = sort(geneList, decreasing = TRUE)

5.运行GSEA分析

从GESA(https://www.gsea-msigdb.org/gsea/downloads.jsp)的官网上，下载一个gmt文件

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

### 看整体分布library(ggplot2)
dotplot(y,showCategory=12,split=\'.sign\')+facet_grid(~.sign)

本次结果中全是激活的

6.特定通路作图

yd <- data.frame(y)

library(enrichplot)
gseaplot2(y,\'HALLMARK_INTERFERON_ALPHA_RESPONSE\',color = \'red\',pvalue_table = T)

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

如需要批量进行几个基因的分析，可以包进一个函数会更方便。

好了，我们又掌握了一个特别强悍，实用的技能。我是进哥哥，有问题随时留言讨论。