100行R语言代码手搓一个KEGG富集分析程序

library(magrittr)
library(dplyr)
library(readr)
library(qvalue)
library(enrichplot)

rm(list= ls())

# 需要富集的基因：207个
data(geneList, package="DOSE")
geneList
gene <-names(geneList)[abs(geneList)>2]
gene
length(gene)
gene_sets <- data.frame(gene_sets=gene)
gene_sets
head(gene_sets)

utils::data(list="kegg_category", package="clusterProfiler")
get("kegg_category", envir = .GlobalEnv)
kegg_category
kegg_category <- select(kegg_category,1:3)

# 下载表格数据
kegg_rest <-function(rest_url){
  message('Reading KEGG annotation online: "', rest_url,'"...')

  reader = readLines
  params =list()
  content <- do.call(reader, args =c(rest_url, params))
  content %<>% strsplit(.,"\t")%>% do.call('rbind', .)
  res <- data.frame(from=content[,1],
                    to=content[,2])
return(res)
}

# 下载 path - gene 列表
url ='https://rest.kegg.jp/link/hsa/pathway'
keggpathid2extid.df <- kegg_rest(url)
keggpathid2extid.df
keggpathid2extid.df[,1]%<>% gsub("[^:]+:","", .)
keggpathid2extid.df[,2]%<>% gsub("[^:]+:","", .)

# 下载 path - name 列表
url ='https://rest.kegg.jp/list/pathway/hsa'
keggpathid2name.df <- kegg_rest(url)
keggpathid2name.df
keggpathid2name.df[,2]<- sub("\\s-\\s[a-zA-Z ]+\\(\\w+\\)$","", keggpathid2name.df[,2])

# 过滤：无名称的通路不用于后续分析
kegg <- inner_join(keggpathid2extid.df, keggpathid2name.df, by ='from')|>
  rename(path=1, gene=2, name=3)|>
  select(1,3,2)|>
  left_join(gene_sets, by =c("gene"="gene_sets"), keep =TRUE)
kegg_enrich <- filter(kegg,!is.na(gene_sets))
kegg
kegg_enrich

# 统计所有基因数, 各通路中基因的个数, 所有富集上的基因数, 富集上基因的通路
N <-length(unique(kegg$gene))
M <- count(kegg, path, name ='M')
n <-length(unique(kegg_enrich$gene_sets))
k <- count(kegg_enrich, path, name ='k')
N
M
n
k
# 整理信息
enrich_count <- left_join(k, M, by ="path")|>
  mutate(O = N-M, n = n)
enrich_count

enrich_stats <- enrich_count |>
  mutate(
    GeneRatio = sprintf("%s/%s", k, n),
    BgRatio = sprintf("%s/%s", M, N),
    RichFactor = k / M,
    FoldEnrichment = RichFactor * N / k,
    mu = M * n / N,
    sigma = mu *(N - n)*(N - M)/ N /(N-1),
    zScore =(k - mu)/sqrt(sigma),
    pvalue = phyper(k-1, M, O, n, lower.tail =FALSE),
    p.adjust = p.adjust(pvalue, method="BH"),
    Count = k
)|>
  arrange(pvalue)|>
  select(-k,-M,-O,-n,-mu,-sigma)
enrich_stats

qobj <- tryCatch(qvalue(p=enrich_stats$pvalue, lambda=0.05, pi0.method="bootstrap"), error=function(e)NULL)
if(inherits(qobj,"qvalue")){
  qvalues <- qobj$qvalues
}else{
  qvalues <-NA
}

enrich_stats$qvalue <- qvalues
enrich_stats

gene_id <- kegg_enrich |>
  group_by(path)|>
  summarize(geneID = paste(gene_sets, collapse ='/'))
gene_id

rval <- enrich_stats |>
  mutate(id=sub("^\\D+","", path))|>
  left_join(gene_id, by ="path")|>
  left_join(keggpathid2name.df, by =c("path"="from"))|>
  left_join(kegg_category, by ="id")|>
  select(-id)|>
  rename(
    ID=path,
    Description=to
)|>
  select(category, subcategory, ID, Description, everything())|>
  relocate(Count, .after = geneID)
rval

write_tsv(rval,'kegg_result.tsv')