GO全称 Gene Ontology中文名叫基因本体论是生信中用于描绘基因功能的标准化词汇表其包括分子功能MF、生物过程BP、细胞组分CC在生信研究过程中我们常会做GO图用于体现细胞、组织富集的功能。而有时我们可能需要进一步深入研究这时我们会根据GO术语尽可能完整的拉出gene list本文为大家提供一个这样的脚本使用时仅需修改第一部分。# # 【模板版本】GO_fetch.R —— 物种 GO 基因清单批量抓取实验室参考脚本 # --------------------------------------------------------------------- # 用途给一份 GO 术语清单面板一次性抓取目标物种下这些术语注释到的基因 # 输出 CSV 供下游使用合并定量表、做桑基图 / 富集分析 / 韦恩图等。 # 怎么用只改「0. 参数区」和「1. GO 术语面板」两处其余代码不用动。 # # 重要说明 # * 面板里的 GO ID 请先在 GO / QuickGO 官网核对术语名、aspect、是否 obsolete # 再运行不要凭记忆写 ID写错不会报错只是抓不到基因。 # * 三种抓取方式自动选择method auto依次尝试 # 1) quickgo QuickGO REST API最全、最新需要 jsonlite 联网 # 2) bioconductor 物种注释包 GO.db离线装过 Bioconductor 就行 # 3) gaf 官方 GAF 注释文件只用 base R 联网兜底 # * use_descendants TRUE 时把 GO 子术语is_a / part_of的注释也算进来 # 口径与 clusterProfiler 的 GOALL / enrichGO 一致推荐 # 改成 FALSE 则只取该 GO 术语本身的直接注释。 # * 全脚本只用 base R不依赖 dplyr / tibble尽量减少环境缺包导致的失败。 # * 换物种要动的地方taxon_id第 0 节、GAF_URLS 里的物种文件第 5 节、 # 以及第 4 节的注释包名默认 org.Mm.eg.db 小鼠 Mus musculus, taxon 10090。 # # 输出文件均写入当前工作目录文件名 out_prefix _ 下面这些默认 out_prefix mouse # go_genes_annotated.csv gene x GO 长表含证据码、注释来源 # target_genes_unique.csv 唯一基因清单Gene 列大写去重 # target_genes_by_category.csv 宽表Gene x 类别(0/1) 命中 GO 证据码 # target_genes_manual.csv 只保留有人工证据非 IEA/IBA/ISS…的基因 # go_panel_counts.csv 面板里每条 GO 抓到的基因数自查用 # go_hits_in_参考表名.csv 仅当 ref_table 存在时输出它与 GO 清单的交集 # 默认参数下文件名与旧版完全一致例如 mouse_go_genes_annotated.csv # # --------------------------------------------------------------------- # 0. 参数区一般只需要改这里 # --------------------------------------------------------------------- method - auto # auto / quickgo / bioconductor / gaf use_descendants - TRUE # TRUE 包含 GO 子术语GOALL 口径 taxon_id - 10090 # 10090 Mus musculus # 输出文件名前缀换物种/换项目就改这里如 bovine、pig留 则不带前缀 out_prefix - mouse # 可选与同目录下的一张参考定量表求交集填文件名NULL 跳过这一步 ref_table - NULL ref_gene_col - 2 # 参考表里基因名在第几列 # 由 out_prefix 生成给下面拼文件名用不用改 fp - if (nzchar(out_prefix)) paste0(out_prefix, _) else # --------------------------------------------------------------------- # 1. GO 术语面板模板示例这里只留了 4 条泛素化 / 去泛素化相关术语 # 用前先把 GO ID 在 GO / QuickGO 官网核对一遍再照下面的格式加油/加行。 # category : 大类标签下游宽表会按它生成 0/1 列可写多类如 ubiquitination / autophagy # go_group : 小类标签只用于汇总打印如 ubiquitin_ligase、DUB、proteolysis # aspect : BP / MF / CC仅作备注不影响抓取 # include : 改成 FALSE 可临时剔除该行不必删 # --------------------------------------------------------------------- GO_terms - data.frame( category c( ubiquitination, # --- 泛素化 --- ubiquitination, ubiquitination, ubiquitination), go_group c( ubiquitination, ubiquitin_ligase, deubiquitination, ubiquitin_proteolysis), go_id c( GO:0016567, # protein ubiquitination GO:0061630, # ubiquitin protein ligase activity GO:0016579, # protein deubiquitination GO:0006511), # ubiquitin-dependent protein catabolic process go_name c( protein ubiquitination, ubiquitin protein ligase activity, protein deubiquitination, ubiquitin-dependent protein catabolic process), aspect c( BP, MF, BP, BP), include TRUE, stringsAsFactors FALSE ) # --------------------------------------------------------------------- # 2. 通用小工具 # --------------------------------------------------------------------- msg - function(...) cat(..., \n, sep ) have_pkg - function(pkg) requireNamespace(pkg, quietly TRUE) # 统一的注释表结构各抓取方式都返回这 6 列 empty_anno - function() { data.frame(panel_go_id character(0), matched_go_id character(0), symbol character(0), qualifier character(0), evidence character(0), assigned_by character(0), stringsAsFactors FALSE) } # 从 data.frame 里安全取一列字符向量缺列就返回 NA col_or_na - function(df, nm, n) { v - df[[nm]] if (is.null(v)) rep(NA_character_, n) else as.character(v) } # 证据码分类这些属于计算/自动注释其余算人工实验证据 computational_codes - c(IEA, IBA, ISS, ISO, ISA, ISM, IGC, RCA) is_manual_evidence - function(ev) { ev - toupper(trimws(as.character(ev))) vapply(strsplit(ifelse(is.na(ev), , ev), [;|,]), function(x) { x - trimws(x) x - x[nzchar(x)] length(x) 0 any(!(x %in% computational_codes)) }, logical(1)) } # --------------------------------------------------------------------- # 3. 方式一QuickGO REST API需要 jsonlite联网 # geneProductTypeprotein taxonId10090分页抓完整个术语的注释 # --------------------------------------------------------------------- quickgo_fetch - function(go_ids, use_desc TRUE, limit 100, verbose TRUE) { if (!have_pkg(jsonlite)) stop(未安装 jsonlite 包) base_url - https://www.ebi.ac.uk/QuickGO/services/annotation/search res_all - list() for (gid in go_ids) { page - 1L pages - NA_integer_ n_term - 0L repeat { q - paste0(base_url, ?goId, gid, taxonId, taxon_id, geneProductTypeprotein, limit, limit, page, page, if (use_desc) goUsagedescendantsgoUsageRelationshipsis_a,part_of else goUsageexact) con - url(q) js - tryCatch(jsonlite::fromJSON(con), error function(e) { try(close(con), silent TRUE); stop(e) }) try(close(con), silent TRUE) df - js$results if (is.null(df)) break if (!is.data.frame(df)) df - as.data.frame(df, stringsAsFactors FALSE) n - nrow(df) if (n 0) break res_all[[length(res_all) 1L]] - data.frame( panel_go_id gid, matched_go_id col_or_na(df, goId, n), symbol col_or_na(df, symbol, n), qualifier col_or_na(df, qualifier, n), evidence col_or_na(df, goEvidence, n), assigned_by col_or_na(df, assignedBy, n), stringsAsFactors FALSE) n_term - n_term n if (is.na(pages)) { pages - if (!is.null(js$pageInfo$total)) as.integer(js$pageInfo$total) else 1L } if (n limit || page pages) break page - page 1L Sys.sleep(0.15) # 对公共接口友好一点 } if (verbose) msg(sprintf( [QuickGO] %-12s %6d 条注释, gid, n_term)) } if (!length(res_all)) return(empty_anno()) do.call(rbind, res_all) } # --------------------------------------------------------------------- # 4. 方式二Bioconductor org.Mm.eg.db离线无需联网 # use_descendants TRUE - 用 GOALL 键 该术语及其所有子术语的基因 # use_descendants FALSE - 用 GO 键仅直接注释 # --------------------------------------------------------------------- bioc_fetch - function(go_ids, use_desc TRUE, verbose TRUE) { if (!have_pkg(org.Mm.eg.db) || !have_pkg(AnnotationDbi)) { stop(未安装 org.Mm.eg.db / AnnotationDbi) } keytype - if (use_desc) GOALL else GO cols - intersect(c(SYMBOL, ONTOLOGY, EVIDENCE), AnnotationDbi::columns(org.Mm.eg.db::org.Mm.eg.db)) res - suppressMessages(AnnotationDbi::select(org.Mm.eg.db::org.Mm.eg.db, keys go_ids, keytype keytype, columns cols)) if (is.null(res) || nrow(res) 0) return(empty_anno()) n - nrow(res) out - data.frame( panel_go_id as.character(res[[keytype]]), matched_go_id if (!is.null(res[[GO]])) as.character(res[[GO]]) else as.character(res[[keytype]]), symbol col_or_na(res, SYMBOL, n), qualifier involved_in, evidence col_or_na(res, EVIDENCE, n), assigned_by org.Mm.eg.db, stringsAsFactors FALSE) if (verbose) msg(sprintf( [org.Mm.eg.db] 共 %d 条 gene-GO 记录, nrow(out))) out } # --------------------------------------------------------------------- # 5. 方式三官方 GAF 注释文件只用 base R 联网兜底方案 # 注意GO 官网已把 GAF 按来源拆分并放进 gaf/ 子目录小鼠对应两个文件 # MOUSE-mod.gaf.gz MGI 人工注释 # MOUSE-uniprot.gaf.gz UniProt / Ensembl / InterPro / GO_Central 等 # 两个都下载后合并本体文件 go-basic.obo 用于展开子术语。 # --------------------------------------------------------------------- GAF_URLS - c( https://current.geneontology.org/annotations/gaf/MOUSE-mod.gaf.gz, https://current.geneontology.org/annotations/gaf/MOUSE-uniprot.gaf.gz ) OBO_URL - https://current.geneontology.org/ontology/go-basic.obo gaf_cols17 - c(DB, DB_Object_ID, DB_Object_Symbol, Qualifier, GO_ID, DB_Reference, Evidence_Code, With_From, Aspect, DB_Object_Name, DB_Object_Synonym, DB_Object_Type, Taxon, Date, Assigned_By, Annotation_Extension, Gene_Product_Form_ID) download_once - function(file, url, verbose TRUE) { if (file.exists(file) file.info(file)$size 0) return(invisible(file)) if (verbose) msg( 正在下载 , basename(file), ...) utils::download.file(url, file, mode wb, quiet !verbose) invisible(file) } # 从 go-basic.obo 里抽出 is_a / part_of 父子关系父 - 子 go_ontology_edges - function(obo_file) { lines - readLines(obo_file, warn FALSE) keep - grepl(^id: GO:, lines) | grepl(^is_a: GO:, lines) | grepl(^(relationship: part_of|part_of): GO:, lines) lines - lines[keep] ids - rep(NA_character_, length(lines)) is_id - grepl(^id: GO:, lines) ids[is_id] - sub(^id: (GO:[0-9]).*$, \\1, lines[is_id]) id_pos - which(!is.na(ids)) if (!length(id_pos)) stop(go-basic.obo 解析失败没有找到任何 term id) block - findInterval(seq_along(lines), id_pos) block_id - rep(NA_character_, length(lines)) ok - block 0 block_id[ok] - ids[id_pos][block[ok]] parent - rep(NA_character_, length(lines)) is_isa - grepl(^is_a: GO:, lines) is_par - grepl(^(relationship: part_of|part_of): GO:, lines) parent[is_isa] - sub(^is_a: (GO:[0-9]).*$, \\1, lines[is_isa]) parent[is_par] - sub(^(relationship: part_of|part_of): (GO:[0-9]).*$, \\2, lines[is_par]) e - !is.na(parent) !is.na(block_id) edges - data.frame(parent parent[e], child block_id[e], stringsAsFactors FALSE) edges - edges[edges$parent ! edges$child, , drop FALSE] unique(edges) } # 某术语的全部子术语is_a / part_of 向下遍历 go_descendants - function(term, edges) { seen - character(0) frontier - term repeat { kids - unique(edges$child[edges$parent %in% frontier]) kids - setdiff(kids, c(term, seen)) if (!length(kids)) break seen - c(seen, kids) frontier - kids } seen } gaf_read_one - function(url, verbose TRUE) { f - basename(url) download_once(f, url, verbose) gaf - utils::read.delim(gzfile(f), header FALSE, sep \t, quote , comment.char , fill TRUE, stringsAsFactors FALSE) if (ncol(gaf) 17) stop(GAF 列数异常, f, 只有 , ncol(gaf), 列GAF 2.x 应为 17 列) if (ncol(gaf) 17) gaf - gaf[, seq_len(17), drop FALSE] names(gaf) - gaf_cols17 gaf - gaf[!is.na(gaf$DB) !grepl(^!, gaf$DB), , drop FALSE] # 去掉文件头 gaf - gaf[!is.na(gaf$Taxon) grepl(10090, gaf$Taxon), , drop FALSE] # 只要小鼠 gaf - gaf[gaf$DB_Object_Type %in% c(protein, gene), , drop FALSE] gaf - gaf[!is.na(gaf$Qualifier) !grepl(NOT, gaf$Qualifier), , drop FALSE] # 去掉 NOT gaf[, c(GO_ID, DB_Object_Symbol, Qualifier, Evidence_Code, Assigned_By)] } # 两个小鼠 GAF 都读进来合并 gaf_read - function(verbose TRUE) { parts - lapply(GAF_URLS, function(u) gaf_read_one(u, verbose)) gaf - do.call(rbind, parts) if (is.null(gaf) || nrow(gaf) 0) stop(GAF 文件里没有小鼠注释) gaf } gaf_fetch - function(go_ids, use_desc TRUE, obo_file go-basic.obo, verbose TRUE) { gaf - gaf_read(verbose) if (nrow(gaf) 0) stop(GAF 文件里没有小鼠注释请删掉 MOUSE-*.gaf.gz 后重新运行) edges - NULL rows - vector(list, length(go_ids)) for (i in seq_along(go_ids)) { ids - go_ids[i] if (use_desc) { if (is.null(edges)) edges - go_ontology_edges(download_once(obo_file, OBO_URL, verbose)) ids - unique(c(ids, go_descendants(ids, edges))) } sub - gaf[gaf$GO_ID %in% ids, , drop FALSE] if (nrow(sub) 0) next rows[[i]] - data.frame(panel_go_id go_ids[i], matched_go_id sub$GO_ID, symbol sub$DB_Object_Symbol, qualifier sub$Qualifier, evidence sub$Evidence_Code, assigned_by sub$Assigned_By, stringsAsFactors FALSE) if (verbose) msg(sprintf( [GAF] %-12s %6d 条注释, go_ids[i], nrow(sub))) } rows - rows[!vapply(rows, is.null, logical(1))] if (!length(rows)) return(empty_anno()) do.call(rbind, rows) } # --------------------------------------------------------------------- # 6. 抓取调度按 method 选择方式auto 时依次尝试三种方式 # --------------------------------------------------------------------- fetch_annotations - function(go_ids, method auto, use_descendants TRUE) { if (identical(method, auto)) { cand - character(0) if (have_pkg(jsonlite)) cand - c(cand, quickgo) if (have_pkg(org.Mm.eg.db) have_pkg(AnnotationDbi)) cand - c(cand, bioconductor) cand - c(cand, gaf) # 兜底只需 base R 联网 } else { cand - method } msg([信息] 依次尝试抓取方式, paste(cand, collapse - )) errs - character(0) for (m in cand) { res - tryCatch({ if (identical(m, quickgo)) quickgo_fetch(go_ids, use_descendants) else if (identical(m, bioconductor)) bioc_fetch(go_ids, use_descendants) else if (identical(m, gaf)) gaf_fetch(go_ids, use_descendants) else stop(未知的 method: , m) }, error function(e) { errs - c(errs, paste0( - , m, 失败: , conditionMessage(e))) NULL }) if (!is.null(res) nrow(res) 0) { res$source - m msg([OK] 抓取成功方式 , m, 共 , nrow(res), 条 gene-GO 注释) return(res) } } stop(所有抓取方式都失败\n, paste(errs, collapse \n), \n请检查网络或先安装 jsonlite / org.Mm.eg.db 后重试。) } # --------------------------------------------------------------------- # 7. 主流程抓取 整理 输出 # --------------------------------------------------------------------- panel - GO_terms[GO_terms$include %in% TRUE, , drop FALSE] if (nrow(panel) 0) stop(GO_terms 里没有任何 include TRUE 的条目) msg() msg(开始抓取 , nrow(panel), 条 GO 术语taxon , taxon_id, , if (use_descendants) 包含子术语 else 仅术语本身, ) msg() anno - fetch_annotations(panel$go_id, method method, use_descendants use_descendants) # 7.1 清洗去空 symbol、去 NOT 注释生成大写基因名 anno$symbol - trimws(as.character(anno$symbol)) anno - anno[!is.na(anno$symbol) nzchar(anno$symbol), , drop FALSE] bad_not - !is.na(anno$qualifier) grepl(NOT, anno$qualifier, ignore.case TRUE) anno - anno[!bad_not, , drop FALSE] anno$gene - toupper(anno$symbol) # 7.2 贴上面板里的分类信息 idx - match(anno$panel_go_id, panel$go_id) anno$category - panel$category[idx] anno$go_group - panel$go_group[idx] anno$go_name - panel$go_name[idx] anno$aspect - panel$aspect[idx] # 7.3 证据码区分人工证据 / 自动计算证据 anno$evidence - toupper(trimws(as.character(anno$evidence))) anno$evidence[is.na(anno$evidence)] - anno$evidence_class - ifelse(is_manual_evidence(anno$evidence), manual, computational) # 7.4 去重、排序、写长表 anno - anno[, c(category, go_group, panel_go_id, go_name, aspect, matched_go_id, symbol, gene, evidence, evidence_class, assigned_by, source, qualifier)] names(anno)[names(anno) panel_go_id] - go_id anno - unique(anno) anno - anno[order(anno$category, anno$go_id, anno$gene, anno$evidence), , drop FALSE] rownames(anno) - NULL anno_file - paste0(fp, go_genes_annotated.csv) write.csv(anno, anno_file, row.names FALSE) msg([写出] , anno_file, (, nrow(anno), 行)) # 7.5 唯一基因清单Gene 列全大写方便与其他表格直接 merge genes - sort(unique(anno$gene)) genes_file - paste0(fp, target_genes_unique.csv) write.csv(data.frame(Gene genes, stringsAsFactors FALSE), genes_file, row.names FALSE) msg([写出] , genes_file, (, length(genes), 个基因)) # 7.6 宽表Gene x 类别(0/1) 命中的 GO 证据码 cats - unique(panel$category) wide - data.frame(Gene genes, stringsAsFactors FALSE) for (cc in cats) wide[[cc]] - as.integer(genes %in% anno$gene[anno$category cc]) cat_go - lapply(cats, function(cc) { m - tapply(anno$go_id[anno$category cc], anno$gene[anno$category cc], function(x) paste(sort(unique(x)), collapse ;)) unname(m[genes]) }) names(cat_go) - paste0(cats, _GO_ids) wide - cbind(wide, as.data.frame(cat_go, stringsAsFactors FALSE)) wide$n_categories - rowSums(wide[, cats, drop FALSE]) m_group - tapply(anno$go_group, anno$gene, function(x) paste(sort(unique(x)), collapse ;)) m_go - tapply(anno$go_id, anno$gene, function(x) paste(sort(unique(x)), collapse ;)) m_ev - tapply(anno$evidence, anno$gene, function(x) paste(sort(unique(x[nzchar(x)])), collapse ;)) manual_genes - sort(unique(anno$gene[anno$evidence_class manual])) wide$evidence_class - ifelse(genes %in% manual_genes, manual, computational) wide$go_groups - unname(m_group[genes]) wide$GO_ids - unname(m_go[genes]) wide$evidence_codes - unname(m_ev[genes]) wide - wide[, c(Gene, cats, n_categories, evidence_class, go_groups, GO_ids, evidence_codes)] wide_file - paste0(fp, target_genes_by_category.csv) write.csv(wide, wide_file, row.names FALSE) msg([写出] , wide_file, (, nrow(wide), 行)) # 7.7 只保留有人工证据的基因剔除纯 IEA/IBA/ISS 等自动注释 manual_tbl - wide[wide$Gene %in% manual_genes, , drop FALSE] manual_file - paste0(fp, target_genes_manual.csv) write.csv(manual_tbl, manual_file, row.names FALSE) msg([写出] , manual_file, (, nrow(manual_tbl), 个基因)) # 7.8 面板统计每条 GO 抓到多少基因方便判断哪条该留、哪条该删 cnt_all - tapply(anno$gene, anno$go_id, function(x) length(unique(x))) cnt_man - tapply(anno$gene[anno$evidence_class manual], anno$go_id[anno$evidence_class manual], function(x) length(unique(x))) panel_out - panel[, c(category, go_group, go_id, go_name, aspect)] n_all - as.integer(cnt_all[panel_out$go_id]); n_all[is.na(n_all)] - 0L n_man - as.integer(cnt_man[panel_out$go_id]); n_man[is.na(n_man)] - 0L panel_out$n_genes - n_all panel_out$n_genes_manual - n_man panel_file - paste0(fp, go_panel_counts.csv) write.csv(panel_out, panel_file, row.names FALSE) msg([写出] , panel_file) # 7.9 可选与一张参考定量表求交集ref_table 留 NULL 或文件不存在就跳过 if (!is.null(ref_table) nzchar(ref_table) file.exists(ref_table)) { mt - tryCatch(utils::read.csv(ref_table, stringsAsFactors FALSE), error function(e) NULL) if (!is.null(mt) ncol(mt) ref_gene_col) { gcol - toupper(trimws(as.character(mt[[ref_gene_col]]))) hit - wide[wide$Gene %in% unique(gcol), , drop FALSE] hits_file - paste0(fp, go_hits_in_, sub(\\.csv$, , basename(ref_table)), .csv) write.csv(hit, hits_file, row.names FALSE) msg([写出] , hits_file, (与 , ref_table, 交集 , nrow(hit), 个基因)) } } else { msg([跳过] ref_table 未填或文件不存在未做交集统计) } # --------------------------------------------------------------------- # 8. 汇总输出 # --------------------------------------------------------------------- msg() msg( 汇总 ) msg(抓取方式 : , paste(unique(anno$source), collapse , )) msg(基因-GO 注释 : , nrow(anno), 行) msg(唯一基因数 : , length(genes)) for (cc in cats) msg(sprintf( %-22s %5d 个基因, cc, sum(wide[[cc]]))) msg(人工证据基因 : , length(manual_genes)) msg() msg(按小类 (go_group) 的唯一基因数) print(tapply(anno$gene, anno$go_group, function(x) length(unique(x)))) msg() msg(说明) msg( * 唯一基因清单, genes_file, 的 Gene 列全大写可直接 merge) msg( * 更严格的清单, manual_file, 剔除纯 IEA/IBA/ISS 等自动注释) msg( * 想缩小范围把 GO_terms 里某行 include 改成 FALSE或 use_descendants 改成 FALSE) msg( * 每条 GO 抓到多少基因见 , panel_file)