GO和KEGG富集分析

写在前面

我们《复现SCI文章系列教程》专栏现在是免费开放，推出这个专栏差不多半年的时间，但是由于个人的精力和时间有限，只更新了一部分。后续的更新太慢了。因此，最终考虑后还是免费开放吧，反正不是什么那么神秘的东西。原本就是一个套路的文章，此外，这篇文章也相对比较简单。在此章节以前，还有一个WGCNA的分析，你若需要可以看**WGCNA分析 | 全流程分析代码**

目前全部开放链接：

SCI文章复现 | GEO文章套路，数据下载和批次效应处理
差异分析和PPI网路图绘制教程

原付费：复现SCI文章系列教程文章

订阅《复现SCI文章系列教程》须知
复现SCI文章系列 | 第一篇文章复现：1. 文章讨论与文章分析套路讲解
2.1 材料与方法 (IF 7.3)
2.2 数据集下载 (IF 7.3)
2.3 数据去重和标准化(附送去批次效应)
2.4 差异分析
2.5 加权基因共表达分析(WGCNA)
2.6 PPI网络分析

本期推文内容

2.7.1 章节总结

在前的教程中，我们已经获得差异基因(2.4 差异分析)和获得与纤维化相关的模块基因。此教程，我们做功能富集分析。但是，此数据问题依旧是很大的影响因素，严重影响后续的分析。

2.7.2 文章结果内容

GO和KEGG富集分析结果
分析结果图

2.7.3 取交集

根据文章分析流程，将DEGs和WGNCA分析获得的结果去交集，获得的交集基因进行后续分析。

在差异分析中，我们获得600多个DEGs，在WGCNA分析中，与纤维化相关的模块为“yellow”。

共有200多个基因。

WGCNA分析全部结果。

2.7.4 获得交集并绘制韦恩图

汇总数据

绘制vennplot，可以是哟TBtools绘制韦恩图，很是方便。或是直接使用我们前提教程进行绘制，基于R语言绘制VennPlot图 | 可以绘制大于等于7个类别的码.

根据筛选出来的关键模块，采用GS>0.2，MM>0.8筛选高度关联基因并与差异基因取交集。我们在自己做分析时，或在写论文时，其实这些参数可以写进论文中，对读者是比较友好的。

但是这里的结果真的是出乎我的意料。我们继续往下看。

绘制维恩图代码

setwd("E:\\小杜的生信筆記\\2023-复现期刊文章系列教程\\复现文章一分析\\05.功能富集分析")
library(VennDiagram) 
library(ggplot2) 
library(venn) 
library(RColorBrewer) 

readFlie=function(input,type,row=T,header=T){ 
  # input 为读入文件的路径，type为读入文件的类型，格式为‘.txt’或‘.csv’,row=T,将文件的第一列设置为列名 
  library(data.table,quietly = TRUE) 
  if(type=='txt'){ 
    dat = fread(input,header = header,sep='\t',stringsAsFactors = F,check.names = F) 
    if(row){ 
      dat = as.data.frame(dat,stringsAsFactors = F) 
      rownames(dat) = dat[,1] 
      dat = dat[,-1] 
    }else{ 
      dat = as.data.frame(dat,stringsAsFactors = F) 
    } 
  }else{ 
    dat = fread(input,header = header,sep=',',stringsAsFactors = F,check.names = F) 
    if(row){ 
      dat = as.data.frame(dat,stringsAsFactors = F) 
      rownames(dat) = dat[,1] 
      dat = dat[,-1] 
    }else{ 
      dat = as.data.frame(dat,stringsAsFactors = F) 
    } 
  } 
  return(dat) 
} 

## 绘制venn图 
wn_venn=function(list,col='black'){ 
  # 定义颜色体系 
  library(RColorBrewer,quietly = TRUE) 
  #corlor = brewer.pal(8,'Dark2') 
  corlor = brewer.pal(8,"Accent") 
  # 绘制Venn图 
  library(VennDiagram, quietly=TRUE) 
  library(venn,quietly = TRUE) 
  if(length(list)<=11){ 
    if(length(list)<=4){ 
      graphics=venn.diagram(list,filename=NULL,fill = corlor[1:length(list)], 
                            col = col,alpha = 0.5, cat.cex = 1.5,rotation.degree = 0) 
      grid.draw(graphics) 
    }else if(length(list)==5){ 
      graphics=venn(list, zcolor = corlor[1:length(list)],box=F,ellipse =TRUE,cexil = 1, cexsn = 1) 
    }else{ 
      graphics=venn(list, zcolor = corlor[1:length(list)],box=F,cexil = 1, cexsn = 1) 
    } 
    return(graphics) 
  }else{ 
    print('The function only supports data of dimension 7 and below.') 
  } 
} 
## 保存图片,只支持ggplot对象 
savePlots=function(path,plot,type=c('pdf','png','tiff')[1],width=10,height=8,dpi=300){ 
  # path表示保存图片路径,需要加上相应的文件扩展名称 
  library(ggplot2) 
  if(type=='pdf'){ 
    ggsave(filename = path,plot = plot,width = width,height = height,device = 'pdf') 
  }else if(type=='png'){ 
    ggsave(filename = path,plot = plot,width = width,height = height,device = 'png',dpi = dpi) 
  }else{ 
    ggsave(filename = path,plot = plot,width = width,height = height,device = 'tiff',dpi = dpi) 
  } 
} 


df <- readFlie("03_DEG_WGCNA交集数据.txt", type = "txt", row = F)
head(df)

df_list = list('DEGs'=sample(df$DEGs),'WGCNA'=sample(df$WGCNA)) 

# 绘制venn图 
## 4维veen图 
venn_2 = wn_venn(df_list[1:2]) 

# 保存图片
savePlots(path = 'DEGs_WGCNA_VennPlot.jpg',plot = venn_2,type = 'jpg',width = 8,height = 8,dpi = 300)

结果

所以基因无任何交集。这个是什么原因呢？？？

现在我们又重新回去做分析。

WGNCA重新分析代码

WGCNA分析 | 全流程分析代码 | 代码一
WGCNA分析 | 全流程分析代码 | 代码二
WGCNA分析 | 全流程代码分享 | 代码三
WGCNA分析 | 全流程分析代码 | 代码四
WGCNA分析 | 全流程分析代码 | 代码五(最新版本)

寻找方案一

我们开始做WGCNA分析时，有过滤数值。

使不过滤的数据进行WGCNA分析。

寻找方案二

WGCNA分析，不去离散样本，所有样本进行分析。

我们这里不在坐重复演练，大体的流程就是这样。我们只需要学会这个思路就可以，不需要和文章一模一样的结果。

正常结果应该是这样的，一定会有交集。我们的只需要调整相关的参数即可。但是更具前面的分析就是离谱。

猜测

数据问题
我有仔细回去看了原文作者使用的数据，这是count值吗？但是这个数据怎么能这么小？以及与我们获得表达量直接是不对等，若是原文作者将FPKM值转换成count值，但依旧是不对等。
数据标准化或归一化的原因
自己认为，我们原始获得数据，已经是进行标准化后的数据了，应该是作者上传时就已经处理过。最大值在13…,但是原文作者依旧是进行标准化处理，那么我们也进行处理吧。
不知道，能力有限的原因？？？？？？

GO富集分析

setwd("E:\\小杜的生信筆記\\2023-复现期刊文章系列教程\\复现文章一分析\\05.功能富集分析")
library(stringr)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
diffSig <- read.table("关键基因ID.txt",header=F,sep="\t",row.names=1)
DEG.gene_symbol = as.character(rownames(diffSig))
DEG.entrez_id = mapIds(x = org.Hs.eg.db,keys = DEG.gene_symbol,
                       keytype = "SYMBOL",column = "ENTREZID")
gene = bitr(DEG.gene_symbol, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
DEG.entrez_id = na.omit(DEG.entrez_id)
DEG.entrez_id
write.csv(gene, "01_78.gene.SYMBOL_WNTREZID.csv")
#write.table(DEG.entrez_id,"DEG.entrez_id.txt",sep="\t")

##  GO 富集分析
ego <- enrichGO(
  gene  = gene$ENTREZID,
  keyType = "ENTREZID",
  OrgDb   = org.Hs.eg.db,
  ont     = "all",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.05,
  readable      = TRUE)

pdf(file = "GO.pdf", width = 10, height = 9)
barplot(ego, drop = TRUE, showCategory = 10,split="ONTOLOGY") +
  scale_y_discrete(labels=function(x) str_wrap(x, width=80))+
  facet_grid(ONTOLOGY~., scale='free')
dev.off()
#
dotplot(ego, showCategory = 10)
# 导出GO数据
write.csv(ego, "05.GO.result.csv")

## KEGG 富集分析
ekegg <- enrichKEGG(
  gene = gene$ENTREZID,
  keyType = "kegg",
  organism  = 'hsa',
  pvalueCutoff  = 0.05,
  pAdjustMethod  = "BH")

pdf(file = "KEGG_bar.pdf", width = 9, height = 8)
barplot(ekegg, showCategory = 10)
dev.off()
pdf(file = "KEGG_dotplot.pdf", width = 9, height = 8)
dotplot(ekegg, showCategory = 10)
dev.off()
## 导出KEGG富集分析数据
write.csv(ekegg, "06.KEGG.result.csv")