R语言KEGG分析

#install.packages('R2HTML')
library('R2HTML')
rm(list=ls())
setwd('D:/Bioinormatics Learning/Plantech 教程资料/Protein analysis/3.KEGG/Rice KEGG')

#读取背景基因（所有注释到KEGG的基因）和差异基因列表，以K号表示

allgene=read.delim('Rice_all_gene_KOs.txt',header=F,sep="\t",stringsAsFactors = F)#读取水稻所有pathways

head(allgene) #展示allgene数据

deggene=read.delim('Degs_KOs.txt',header=F,sep="\t",stringsAsFactors = F) #取差异基因的pathway信息

head(deggene) #展示deggene数据


pathinfo=read.delim('kegg.pathway.id',header=T,sep="\t",stringsAsFactors = F) #取所有的pathway信息

head(pathinfo) #展示pathinfo数据

#读取背景基因（所有注释到KEGG的基因）和差异基因列表，以K号表示

pathway=unique(pathinfo[,c('PID','PATHWAY')]) #unique()去除重复函数，获得单一的Pathway用于后面for循环计数每个pathway的Ko数目

head(pathway) #展示pathway数据

allgene=unique(allgene$V1)   #unique()去除重复函数,获得单一的Pathway用于后面for循环计数每个pathway的Ko数目

head(allgene) #展示allgene数据

deggene=unique(deggene$V1)   #unique()去除重复函数，获得单一的Pathway用于后面for循环计数每个pathway的Ko数目

head(deggene)  #展示deggene数据

#==============设定输出的文件名称============

outfile="kegg.enrichment"

#=======将基因与Pathway中有注释的基因取交集========

deginfo=pathinfo[pathinfo[,1]%in%deggene,] #相当于merge

head(deginfo)#展示deginfo数据

allinfo=pathinfo[pathinfo[,1]%in%allgene,]  #相当于merge

head(allinfo) #展示allinfo数据

k=length(unique(deginfo[,1]))  #差异基因总数，length函数求基因数目

k  #展示K数据

N=length(unique(allinfo[,1]))  #总基因数   #length函数求基因数目

N  #展示N数据

#=========对每个通路基于超几何分布计算富集p值=======

p<-data.frame() #新建data frame用于储存KO数目


urlColor<-NULL # used to color the genes in pathway
geneLink<-NULL # used to set hyperlink to DEGs in pathway

#通过循环计算每个KEGGpathways上的Ko号

for (i in seq(1,nrow(pathway))) {
 
  pid=pathway[i,'PID'] #获取Pathway的KO号码，例如：pathway[1,'PID']=》 "ko00010"；pathway[1,'PID']=》"ko00020"等
  
  allInPath=allinfo[allinfo[,'PID']%in%pid,] #获得每一个pathway所有的KO数目，具体见备注
  
  m=length(unique(allInPath[,1]))  #计算i Pathway基因总数
  
  degInPath=deginfo[deginfo[,'PID']%in%pid,] #获得每一个pathway所有的KO数目
  
  x=length(unique(degInPath[,1]))  #计算i Pathway中差异基因数
  
  p[i,1]=x #数据框p的第i行，第一列值
  p[i,2]=m #数据框p的第i行，第二列值
  p[i,3]=phyper(x-1,m,N-m,k,lower.tail=FALSE)  #lower.tail=FALSE,计算的是P[X>(x-1)]的概率,其中x-1：某一pathway中差异基因数目（或KO数），m：背景数（或KO数），N:所有基因数（或KO数），K：所有差异数（或KO数）
  
  urlColor[i]=apply(as.matrix(paste("/",t(degInPath[,1]),'%09red',sep="")),2,paste,collapse="")
  Link=paste('<a href=', shQuote(paste('https://www.kegg.jp/dbget-bin/www_bget?ko:',degInPath[,1],sep="")),'/','>', degInPath[,1],'</a>',sep="")
  if(x==0){
    geneLink[i]="--"
  }else{
    geneLink[i]=apply(as.matrix(Link),2,paste,collapse=", ")
  }
}

#============计算FDR，整理结果==============

fdr=p.adjust(p[,3],method="BH") #利用第三列的p值计算FDR值，！！！可以借鉴！！！！

output=cbind(pathway,p,fdr)#合并数据库框，产生新的数据框,因为for循环按照pathway内容从上到下依次运行的，所以顺序是相同的，可以进行列合并

head(output)#展示数据内容

colnames(output)=c('ID','Pathway','DEGsInPathway','GenesInPathway','Pvalue','FDR') #重新给数据框命名

ind=order(output[,"Pvalue"])#按照P值从小到大对结果排序
output=output[ind,]       #用于输出到表格文件
urlColor=urlColor[ind]
DEGs=geneLink[ind]

output2=cbind(output,DEGs) #用于输出到HTML文件

#==============导出到表格文件================

write.table(output,file=paste(outfile,".xls",sep=""),quote = F,row.names = F,col.names = T,sep="\t")

#==============导出结果到HTML文件============

#==============导出结果到HTML文件============
urlTitles <- output[,'ID']
url=paste('https://www.genome.jp/kegg-bin/show_pathway?',gsub("ko", "map", urlTitles),urlColor,sep="")

htmlOUT <- transform(output2, ID = paste('<a href = ', shQuote(url),'/','>', urlTitles, '</a>'))

target <- HTMLInitFile(Title = "KEGG Pathway Ernichment Results",outdir=getwd(),filename=outfile, BackGroundColor="#f8f8f8",useGrid = T,useLaTeX = F,HTMLframe = F)
write("<style>table{border-collapse:collapse;width:1280px;}a:link,a:visited {color: #3366FF;text-decoration: none; }a:hover,a:active {color: #ff3366;text-decoration: none; };</style>",target,append = T)
HTML(as.title("KEGG Pathway Ernichment Results"),file=target)
HTML(htmlOUT,file=target,innerBorder = 1,row.names = F,digits=3)
HTMLEndFile()
#============================================

<big>#------------------备注-------------------<big>

#------------------备注-------------------
pid=pathway[2,'PID']#获得第二个pathway

pid #pid内容

allInPath=allinfo[allinfo[,'PID']%in%pid,]#计算第二个pathway中所有背景Ko号

allInPath #展示allInPath内容

备注：超几何分布

image.png

参照Plantech 植物生物信息综合课程
详细信息可以观看https://ke.qq.com/course/1645755