Usage

  1. Usage:
  3.      enricher(
  4.        gene,
  5.        pvalueCutoff = 0.05,
  6.        pAdjustMethod = "BH",
  7.        universe,
  8.        minGSSize = 10,
  9.        maxGSSize = 500,
  10.        qvalueCutoff = 0.2,
  11.        TERM2GENE,
  12.        TERM2NAME = NA
  13.      )
  15. Arguments:
  17.     gene: a vector of gene id
  19. pvalueCutoff: pvalue cutoff on enrichment tests to report
  21. pAdjustMethod: one of "holm", "hochberg", "hommel", "bonferroni", "BH",
  22.           "BY", "fdr", "none"
  24. universe: background genes. If missing, the all genes listed in the
  25.           database (eg TERM2GENE table) will be used as background.
  27. minGSSize: minimal size of genes annotated for testing
  29. maxGSSize: maximal size of genes annotated for testing
  31. qvalueCutoff: qvalue cutoff on enrichment tests to report as
  32.           significant.  Tests must pass i) ‘pvalueCutoff’ on unadjusted
  33.           pvalues, ii) ‘pvalueCutoff’ on adjusted pvalues and iii)
  34.           ‘qvalueCutoff’ on qvalues to be reported.
  36. TERM2GENE: user input annotation of TERM TO GENE mapping, a data.frame
  37.           of 2 column with term and gene
  39. TERM2NAME: user input of TERM TO NAME mapping, a data.frame of 2 column
  40.           with term and name

Example

  1. #!/usr/bin/env Rscript
  2. library(ggplot2)
  3. library(optparse)
  4. suppressMessages(library(showtext))
  5. suppressMessages(library(clusterProfiler))
  7. description <- "
  8.     \033[1;32mDifferential Enrichment Analysis and Plot for KEGG/GO\033[0m
  9. "
  11. # =========================================== 参数解析 =========================================================
  12. option_list = list(
  13.   make_option(c('-b', '--background'), help = 'The background file contains geneid and function information'),
  14.   make_option(c('-g', '--diff_genes'), help = 'The different genes'),
  15.   make_option(c('-t', '--type'), help = 'The type of your function, choices: kegg, go'),
  16.   make_option(c('-x', '--ontology'), help = 'The ontology for GO enrichment, choices: BP, CC, MF'),
  18.   make_option(c('-o', '--outdir'), help = 'The output directory [default: %default]', default = '.'),
  19.   make_option(c('-n', '--number'), help = 'The max number used to plot [default: %default]', default = 10, type='integer'),
  20.   make_option(c('-p', '--pcutoff'), help = 'The pvalue(p.adjust) cutoff [default: %default]', default = 0.05, type='double'),
  21.   make_option(c('-m', '--method'), help = 'The adjustment method [default: %default]', default = 'BH'),
  23.   make_option(c('-T', '--title'), help = 'The title of the picture'),
  24.   make_option(c('-X', '--xlab'), help = 'The xlab of the picture'),
  25.   make_option(c('-Y', '--ylab'), help = 'The ylab of the picture'),
  26.   make_option(c('-W', '--width'), type = 'integer', help = 'The width of the picture'),
  27.   make_option(c('-H', '--height'), type = 'integer', help = 'The height of the picture')
  29. )
  31. parser <- OptionParser(description = description,
  32.                        option_list = option_list,
  33.                        epilogue = 'contact: suqingdong@novogene.com')
  34. args <- parse_args(parser)
  36. if (is.null(args$background) || is.null(args$diff_genes)) {
  37.   print_help(parser)
  38.   quit('no')
  39. }
  40. # ============================================================================================================
  42. # # just for testing
  43. # args = data.frame(diff_genes = 'diff_gene.list',
  44. #                   background = 'geneid_kegg.tsv',
  45. #                   pcutoff = 0.05,
  46. #                   method = 'BH')
  47. # print(args)
  48. type = args$type
  49. ontology = args$ontology
  50. pcutoff = args$pcutoff
  51. method = args$method
  52. number = args$number
  55. # showtext
  56. # showtext_auto()
  57. # font_add('Times New Roman', './times.ttf')
  60. # 差异基因：每行一个基因ID
  61. diff_genes <- read.table(args$diff_genes)[, 1]
  64. # 背景基因：geneID <TAB> Term
  65. bg_data <- read.table(args$background, sep = '\t')
  67. colnames(bg_data) = c('geneid', 'term')
  69. # term去重，避免重复获取
  70. term_list <- unique(bg_data$term)
  72. # type未指定是自动判断
  73. if (is.null(args$type)) {
  74.   if (startsWith(term_list[1], 'GO:')) {
  75.     type = 'go'
  76.   } else {
  77.     type = 'kegg'
  78.   }
  79. }
  82. # 添加注释
  83. if (type == 'kegg') {
  84.   ko_names <- ko2name(term_list)
  85.   merged_data <- merge(bg_data, ko_names, by.x = 'term', by.y = 'ko')
  86.   term2gene <- data.frame(merged_data$term, merged_data$geneid)
  87.   term2name <- data.frame(merged_data$term, merged_data$name)
  88. } else {
  89.   go_terms <- go2term(term_list)
  90.   go_onts <- go2ont(term_list)
  91.   md <- merge(bg_data, go_terms, by.x = 'term', by.y = 'go_id')
  92.   merged_data <- merge(md, go_onts, by.x = 'term', by.y = 'go_id')
  94.   if( ! is.null(ontology) ) {
  95.     merged_data <- merged_data[merged_data$Ontology==ontology, ]
  96.     type <- paste(type, ontology, sep='.')
  97.   }
  99.   term2gene <- data.frame(merged_data$term, merged_data$geneid)
  100.   term2name <- data.frame(merged_data$term, merged_data$Term)
  101. }
  104. # 富集分析
  105. result = enricher(
  106.   diff_genes,
  107.   TERM2GENE = term2gene,
  108.   TERM2NAME = term2name,
  109.   pvalueCutoff = pcutoff,
  110.   qvalueCutoff = 1,
  111.   pAdjustMethod = method  # BH == FDR
  112. )
  114. # 结果转数据框
  115. result_df = data.frame(result)
  117. # 按p.adjust正序排序, 分值小的在上
  118. result_df = result_df[order(result_df$p.adjust, decreasing = F), ]
  120. # 保存文件
  121. outfile = sprintf('%s.result.xls', type)
  122. write.table(result_df, outfile, sep='\t', quote = F, row.names = F)
  125. # 截取前N行
  126. result_df = head(result_df, number)
  128. # Ratio列转换成数值
  129. calc_ratio = function(x) {
  130.   res = as.numeric(unlist(strsplit(x, "/")))
  131.   res[1] / res[2]
  132. }
  133. result_df$GeneRatio = unlist(lapply(result_df$GeneRatio, calc_ratio))
  135. # Description列设置因子，用于控制显示顺序
  136. # 按p.adjust逆序排序, 分值小的图片中在上
  137. result_df = result_df[order(result_df$p.adjust, decreasing = T), ]
  138. result_df$Description = factor(result_df$Description, levels=result_df$Description)
  140. # 样式
  141. family = 'Times New Roman'
  142. size = 28
  143. angle = 0
  145. theme.style <- theme(
  146.   text = element_text(family = family),
  147.   axis.text = element_text(size = size),
  148.   axis.text.y = element_text(angle = angle),
  149. )
  152. # 绘图
  153. # barplot
  154. outbar = sprintf('%s.barplot.png', type)
  156. barplot <-
  157.   ggplot(result_df, aes(x = Count, y = Description, fill = p.adjust)) +
  158.   geom_bar(stat = 'identity') +
  159.   theme_bw() +
  160.   theme.style +
  161.   scale_fill_gradient(low = "red", high = "blue")
  163. ggsave(outbar, barplot, type = 'cairo')
  164. sprintf('save file: %s', outbar)
  167. # dotplot
  168. outdot = sprintf('%s.dotplot.png', type)
  170. dotplot <- 
  171.   ggplot(result_df, aes(x=GeneRatio, y=Description, color=p.adjust, size=Count)) +
  172.   geom_point() +
  173.   scale_color_gradient(low="red", high="blue") +
  174.   theme_bw() +
  175.   theme.style
  177. ggsave(outdot, dotplot, type='cairo')
  178. sprintf('save file: %s', outdot)