ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

差异基因富集分析(R语言——GOKEGGGSEA)

差异基因富集分析(R语言——GOKEGGGSEA) 接着上次的内容上篇内容给大家分享了基因表达量怎么做分组差异分析从而获得差异基因集想了解的可以去看一下这篇主要给大家分享一下得到显著差异基因集后怎么做一下通路富集。#### 1.准备差异基因集我就直接把上次分享的拿到这边了。我们一般都把差异基因分为上调基因和下调基因分别做通路富集分析。下面上代码可能包含我的一些个人习惯勿怪。显著差异基因的筛选条件根据个人需求设置哈。 ##载入所需R包 library(readxl) library(DOSE) library(org.Hs.eg.db) library(topGO) library(pathview) library(ggplot2) library(GSEABase) library(limma) library(clusterProfiler) library(enrichplot) ##edger edger_diff - diff_gene_Group edger_diff_up - rownames(edger_diff[which(edger_difflogFC0.584962501),])edgerdiffdown−rownames(edgerdiff[which(edgerdifflogFC 0.584962501),]) edger_diff_down - rownames(edger_diff[which(edger_difflogFC0.584962501),])edgerd​iffd​own−rownames(edgerd​iff[which(edgerd​ifflogFC -0.584962501),]) ##deseq2 deseq2_diff - diff_gene_Group2 deseq2_diff_up - rownames(deseq2_diff[which(deseq2_difflog2FoldChange0.584962501),])deseq2diffdown−rownames(deseq2diff[which(deseq2difflog2FoldChange 0.584962501),]) deseq2_diff_down - rownames(deseq2_diff[which(deseq2_difflog2FoldChange0.584962501),])deseq2d​iffd​own−rownames(deseq2d​iff[which(deseq2d​ifflog2FoldChange -0.584962501),]) ##将差异基因集保存为一个list gene_diff_edger_deseq2 - list() gene_diff_edger_deseq2[[“edger_diff_up”]] - edger_diff_up gene_diff_edger_deseq2[[“edger_diff_down”]] - edger_diff_down gene_diff_edger_deseq2[[“deseq2_diff_up”]] - deseq2_diff_up gene_diff_edger_deseq2[[“deseq2_diff_down”]] - deseq2_diff_down #### 2.进行通路富集分析这里主要介绍普通的GOKEGGGSEA的简单富集。筛选显著富集通路的筛选条件也是根据自己的需求决定一般是矫正后P值小于0.05。我这里是省事写了各list循环。 for (i in 1:length(gene_diff_edger_deseq2)){ keytypes(org.Hs.eg.db) entrezid_all mapIds(x org.Hs.eg.db, keys gene_diff_edger_deseq2[[i]], keytype “SYMBOL”, #输入数据的类型 column “ENTREZID”)#输出数据的类型 entrezid_all na.omit(entrezid_all) #na省略entrezid_all中不是一一对应的数据情况 entrezid_all data.frame(entrezid_all) ##GO富集## GO_enrich enrichGO(gene entrezid_all[,1], OrgDb org.Hs.eg.db, keyType “ENTREZID”, #输入数据的类型 ont “ALL”, #可以输入CC/MF/BP/ALL #universe 背景数据集我没用到它。 pvalueCutoff 1,qvalueCutoff 1, #表示筛选的阈值阈值设置太严格可导致筛选不到基因。可指定 1 以输出全部 readable T) #是否将基因ID映射到基因名称。 GO_enrich_data data.frame(GO_enrich) write.csv(GO_enrich_data,paste(‘GO_enrich_’,names(gene_diff_edger_deseq2)[i], ‘.csv’, sep “”)) GO_enrich_data - GO_enrich_data[which(GO_enrich_dataKaTeX parse error: Expected group after _ at position 66: …aste(GO_enrich_̲,names(gene_di…p.adjust 0.05),] write.csv(KEGG_enrich_data, paste(‘KEGG_enrich_’,names(gene_diff_edger_deseq2)[i], filter.csv’, sep “”)) } #### 3.通路富集情况可视化这里只介绍一种简单的气泡图当然还有其他的自己去了解吧。 ##GOKEGG富集BPCCMFKEGG分面绘图需要分开处理一下富集结果里的ONTOLOGYL列修改 GO_enrich_data_BP - subset(GO_enrich_data, subset GO_enrich_dataONTOLOGYBP)GOenrichdataCC−subset(GOenrichdata,subsetGOenrichdataONTOLOGY BP) GO_enrich_data_CC - subset(GO_enrich_data, subset GO_enrich_dataONTOLOGYBP)GOe​nrichd​ataC​C−subset(GOe​nrichd​ata,subsetGOe​nrichd​ataONTOLOGY “CC”) GO_enrich_data_MF - subset(GO_enrich_data, subset GO_enrich_dataKaTeX parse error: Expected EOF, got # at position 26: … MF) #̲#提取GO富集BPCCMF的t…ONTOLOGY - “KEGG” KEGG_enrichresult - KEGG_enrichresult[,c(10,1:9)] ##整合GO KEGG富集结果 ego_GO_KEGG - GO_enrich ego_GO_KEGGresult - rbind(ego_GO_KEGGresult, KEGG_enrichresult[1:5,]) ego_GO_KEGGresultONTOLOGY−factor(egoGOKEGGresultONTOLOGY - factor(ego_GO_KEGGresultONTOLOGY−factor(egoG​OK​EGGresultONTOLOGY, levels c(“BP”, “CC”, “MF”,“KEGG”))##规定分组顺序 ##简单画图 pdf(“edger_diff_up_dotplot.pdf”, width 7, height 7) dotplot(ego_GO_KEGG, split “ONTOLOGY”, title“UP-GOKEGG”, label_format 60, color “pvalue”) facet_grid(ONTOLOGY~., scale “free_y”) theme(plot.title element_text(hjust 0.5, size 13, face “bold”), axis.text.x element_text(angle 90, hjust 1)) dev.off() #### 4.气泡图如图所示做了些处理真实图片左侧pathway是跟后面气泡一一对应的当然还有其他可视化方式那就需要各位自己去探索了谢谢5.GSEA富集分析这里也是做一下简单的GSEA ##GSEA官方网站下载背景gmt文件并读入 geneset - list() geneset[[“c2_cp”]] - read.gmt(“c2.cp.v2023.2.Hs.symbols.gmt”) geneset[[“c2_cp_kegg_legacy”]] - read.gmt(“c2.cp.kegg_legacy.v2023.2.Hs.symbols.gmt”) geneset[[“c2_cp_kegg_medicus”]] - read.gmt(“c2.cp.kegg_medicus.v2023.2.Hs.symbols.gmt”) geneset[[“c2_cp_reactome”]] - read.gmt(“c2.cp.reactome.v2023.2.Hs.symbols.gmt”) geneset[[“c3_tft”]] - read.gmt(“c3.tft.v2023.2.Hs.symbols.gmt”) geneset[[“c4_cm”]] - read.gmt(“c4.cm.v2023.2.Hs.symbols.gmt”) geneset[[“c5_go_bp”]] - read.gmt(“c5.go.bp.v2023.2.Hs.symbols.gmt”) geneset[[“c5_go_cc”]] - read.gmt(“c5.go.cc.v2023.2.Hs.symbols.gmt”) geneset[[“c5_go_mf”]] - read.gmt(“c5.go.mf.v2023.2.Hs.symbols.gmt”) geneset[[“c6”]] - read.gmt(“c6.all.v2023.2.Hs.symbols.gmt”) geneset[[“c7”]] - read.gmt(“c7.all.v2023.2.Hs.symbols.gmt”) ##进行GSEA富集分析这里也是写了个循环 gsea_results - list() for (i in names(gene_diff)){ geneList - gene_diff[[i]]KaTeX parse error: Expected EOF, got } at position 487: …s0) } }̲ ##批量绘图…ID[j], title gsea_results[[i]]resultID[j])pdf(paste(paste(names(gsearesults)[i],gsearesults[[i]]resultID[j]) pdf(paste(paste(names(gsea_results)[i], gsea_results[[i]]resultID[j])pdf(paste(paste(names(gsear​esults)[i],gsear​esults[[i]]resultID[j], sep “),”.pdf,sep “”)) print§ dev.off() } 6.GSEA富集最简单图形如下分享到此结束了希望对大家有所帮助。
返回列表