GSEA基因基富集分析

r 同时被 2 个专栏收录
14 篇文章 0 订阅
6 篇文章 0 订阅

###GSEA富集分析中,不需要提取差异基因,只需要将所有基因的表达情况按照一定顺序排列(一般按log2FD)之后根据对照组和实验组中所有基因在红色(蓝色)富集,从而得出对照组或者实验组所富集到的通路。因此,GSEA分析可以很有效的避免基因过滤从而导致一些基因被筛除掉。

library(DO.db)
require(DOSE)
library(clusterProfiler)
library(AnnotationHub)
library(readr)
library(pheatmap)
library(tidyverse)
library(DESeq2)
library(ggplot2)
library(export)
library(enrichplot)
library(Rgraphviz)
library(org.Hs.eg.db)
#all_entrez.csv的第一列是entrezid,第二列是FoldChange的值。

DEG <- read.csv(file.choose(),sep = "\t",header = TRUE)
fix(DEG)
DEG1 <- DEG$X  #第一列向量化
DEG.gene_symbol = as.character(DEG1)#选择基于列表并且将其向量化

#进行基因名称转换匹配

DEG.entrize_id <- mapIds(x= org.Hs.eg.db,
                         keys = DEG.gene_symbol,
                         keytype = "SYMBOL",
                         column = "ENTREZID")
luffy <- data.frame(DEG.entrize_id)
fix(luffy)
luffy$lgfc = " "
luffy$lgfc = DEG$log2FoldChange
write.csv(DEG.entrize_id,file = "GESA-GO_CC.csv")

#去除未匹配到的NA值

DEG.enter_id <- na.omit(DEG.entrize_id)

#gse需要单独做数据格式

geneList <- luffy[,2]
names(geneList)=as.factor(luffy[,1])
geneList <- sort(geneList,decreasing = TRUE)
fix(geneList)

#gseGO进行GSEA分析

###gseBP <- gseGO(geneList=geneList,ont="BP",OrgDb=maize,keyType = 'ENTREZID',nPerm = 50000,minGSSize = 100,maxGSSize = 6000,pvalueCutoff = 0.05,verbose = FALSE)

############# GSEA CC 模式 start

ego3 <- gseGO(geneList = geneList,OrgDb = org.Hs.eg.db,ont = "BP",nPerm = 1000,minGSSize = 100,maxGSSize = 1000,pvalueCutoff = 0.05,verbose = FALSE)
write.csv(ego3,file = "GESA-GO_CC.csv")

#显示前4组信息

gseaplot2(ego3,geneSetID = 1:4, ES_geom = "dot",pvalue_table = TRUE)

############# GSEA BP 模式 start

ego2 <- gseGO(geneList = geneList,OrgDb = org.Hs.eg.db,ont = "CC",pvalueCutoff = 0.05,verbose = FALSE)
write.csv(ego2,file = "GESA-GO_BP.csv")
#ridgeline plot for expression distribution of GSEA result
ridgeplot(ego2)
out_img(filename = "ridgeplot_BP",pic_width = 12,pic_height = 12)

#只显示值最高的一组的信息

gseaplot2(ego2,geneSetID = 1:4, ES_geom = "dot",pvalue_table = TRUE)

############# GSEA MF 模式 start

ego4 <- gseGO(geneList = geneList,OrgDb = org.Hs.eg.db,ont = "MF",pvalueCutoff = 0.05,verbose = FALSE)
write.csv(ego4,file = "GESA-GO_MF.csv")
#ridgeline plot for expression distribution of GSEA result
ridgeplot(ego4)
out_img(filename = "ridgeplot_MF",pic_width = 12,pic_height = 12)

#显示前4组信息

gseaplot2(ego4,geneSetID = 1:4, ES_geom = "dot",pvalue_table = TRUE)
  • 0
    点赞
  • 1
    评论
  • 2
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

GSEA富集分析,1、准备三个文件第一行:#1.2,表示版本号,自己准备文件时照抄就行; 第二行:两个数分别表示gene NAME的数量和样本数量(矩阵列数-2); 矩阵:第一列是NAME;第二列Description,没有的话可以全用na或任意字符串填充;后面的就是基因在不同样本中标准化后的表达数据了 (部分统计量metrics for ranking genes计算需要log转换后的数据,后面会有提及。其它情况是否为log转换的数据都可用,GSEA关注的是差异,只要可比即可)。 #其次是样品分组信息(通常用.gmt作为后缀) 第一行:三个数分别表示:34个样品,2个分组,最后一个数字1是固定的; 第二行:以#开始,tab键分割,分组信息(有几个分组便写几个,多个分组在比较分析时,后面需要选择待比较的任意2组);(样品分组中NGT表示正常耐糖者,DMT表示糖尿病患者,自己使用时替换为自己的分组名字) 第三行:样本对应的组名。样本分组信息的第三行,同一组内的不同重复一定要命名为相同的名字,可以是分组的名字。例如相同处理的不同重复在自己试验记录里一般是Treat6h_1、Treat6h_2、Treat6h_3,但是在这里一定都要写成一样的值Treat6h。与表达矩阵的样品列按位置一一对应,名字相同的代表样品属于同一组。如果是样本分组信息,上图中的0和1也可以对应的写成NGT和DMT,更直观。但是,如果想把分组信息作为连续表型值对待,这里就只能提供数字。 3. 预定义基因集(gmx or gmt)——非必需文件(需要注意第一列的基因集名称必须是唯一的) 通常用.gmt作为后缀。若采用GSEA预定义的MSigDB数据库中的功能基因分析,则无需自己定义该文件。每一行为一个功能基因集,第一列为基因集的名称,第二列为简单描述,第三列及以后列为该功能基因集所包含的基因symbol。基因集包含多少个基因,就列出多少个基因。文件以tab作为分隔符。
相关推荐
©️2020 CSDN 皮肤主题: 数字20 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值