Posted on

有这样的使用场景么?

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

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 <- 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跟阿拉法干扰素正相关,这个事情没什么好说的吧。

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

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

##对于多个通路绘制在一起:::
pathway=c("KEGG_ENDOCYTOSIS","KEGG_MAPK_SIGNALING_PATHWAY")
gseaplot2(y,pathway,color = c('red',"blue"),pvalue_table = F)

11 Replies to “基于单基因批量相关性分析的GSEA”

  1. error in if(abs(max.ES)>abs(min.ES)){:missing value where TRUE/FALSE needed,请问老师,这个是什么意思呢,最后的出图出不来。望老师帮帮我。

    1. 这样看不出来诶。。。确定你所要画的通路在富集结果里面吗 并且拼写正确?可以加我微信具体讨论

    1. 您好 发给你数据太大,这个不一定需要TPM数据,只要是任意标准化的表达矩阵都可以,你可以从XENA上下载任意肿瘤的表达矩阵进行分析,你先试试 有问题加我微信讨论

  2. 您好!非常感谢您提供便捷的学习资料。想问下您:我在运行“plan(multiprocess);system.time(dd <- batch_cor('PDCD1'))”时,显示“ Error: ‘node$session_info$process$pid == pid’ is not TRUE ” ,这个问题应该怎么解决。谢谢

          1. 你好,王老师,多个通路的颜色没办法自己绘制,是不是要加个颜色的代码。

发表评论

邮箱地址不会被公开。 必填项已用*标注