r语言：基于单基因批量相关性分析的gsea的使用场景
1.已经确定研究的基因，但是想探索他潜在的功能，可以通过跟这个基因表达最相关的基因来反推他的功能，这种方法在英语中称为guilt of association，协同犯罪。
2.我们的注释方法依赖于tcga大样本，既然他可以注释基因，那么任何跟肿瘤相关的基因都可以被注释，包括长链非编码rna。
这个方法以前阐述过：
单基因批量相关性分析的妙用
但是这个方法有个小缺陷，并不知道最后富集的通路是正向影响还是反向影响，也就是无法判断方向。判断方向的工具也不是没有，gsea就是一个。所以，我想能不能把批量相关性分析和gsea结合起来。
gsea需要的gene set是现成的没有问题，但是genelist没有，这里我们可以把所有基因跟单个基因的相关性系数当做logfc，有正有负，就解决了genelist的问题。这个想法不是我的，是我的一个学员的，不过他要解决的是microrna把基因的问题。
下面来实战一下：
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
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分析
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跟阿拉法干扰素正相关，这个事情没什么好说的吧。
