{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# An exploratory analysis of Gene Ontology (GO) terms predictions","metadata":{}},{"cell_type":"markdown","source":"<h3>What we do here:</h3>\n<div style=\"line-height:24px; font-size:16px\">\n    <ul style=\"list-style:circle\">\n<li>Examine the distribution of GO terms in the submisssion file\n<li>Visualise GO terms of a single protein in the form of a directed acyclic graph (DAG) - ancestor chart\n<li>Visualise separate  branches of predicted GO terms in the DAG\n    </ul>\n</div>","metadata":{}},{"cell_type":"markdown","source":"<hr>\n\n<h3>Data</h3> \n\n<div style=\"line-height:24px; font-size:14px\">\n\n\nTSV/qs file with CAFA5 submission and go-basic.obo provided by organizers.\n \n</div>\n<h3>Output</h3>  \n    \n<div style=\"line-height:24px; font-size:14px\">\n   HTML file with interactive ancestor chart of all terms of the selected protein (make your own by simply changing the name of the protein).\n</div>\n<hr>","metadata":{}},{"cell_type":"markdown","source":"**Load packages**","metadata":{}},{"cell_type":"code","source":"suppressPackageStartupMessages({\n    \n    library(stringr) #wrap text into paragraphs\n    library(data.table)\n    library(qs)\n    install.packages(\"ontologyIndex\")\n    library(ontologyIndex)    \n    library(visNetwork) #interactive plots\n})","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-08-23T21:39:13.225307Z","iopub.execute_input":"2023-08-23T21:39:13.248863Z","iopub.status.idle":"2023-08-23T21:39:20.44242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"compute_stats <- function(dt) {\n  data.table(\n    \"Unique proteins\" = length(unique(dt$EntryID)),\n    \"Unique terms\" = length(unique(dt$term)),\n    \"Min terms per prot\" = dt[, .(n_per_prot = .N), by = EntryID][, min(n_per_prot)],\n    \"Max terms per prot\" = dt[, .(n_per_prot = .N), by = EntryID][, max(n_per_prot)],\n    \"Number of predictions\" = nrow(dt),\n    \"Min.\" = round(min(dt[, pred]), 3),\n    \"Median\" = round(median(dt[, pred]), 3),\n    \"Mean\" = round(mean(dt[, pred]), 3),\n    \"Max.\" = round(max(dt[, pred]), 3)\n  )\n}","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-08-23T21:39:20.444536Z","iopub.execute_input":"2023-08-23T21:39:20.445584Z","iopub.status.idle":"2023-08-23T21:39:20.455283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load ontology**","metadata":{}},{"cell_type":"code","source":"ontology <- get_ontology(\n  \"/kaggle/input/cafa-5-protein-function-prediction/Train/go-basic.obo\",\n  propagate_relationships = c(\"is_a\", \"part_of\"),\n  extract_tags = \"everything\"\n)\nterms <- data.table(\n  term = ontology$id, \n  aspect = unlist(ontology$namespace)\n)","metadata":{"execution":{"iopub.status.busy":"2023-08-23T21:39:20.45725Z","iopub.execute_input":"2023-08-23T21:39:20.458247Z","iopub.status.idle":"2023-08-23T21:39:35.372409Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load CAFA5 submission file**","metadata":{}},{"cell_type":"code","source":"submission_file <- '/kaggle/input/cafa5-best-submits/mean_corr_filt_by_tax_v146_49_60_68_69_121_corr_go_mean450_corr.qs'\n\nif (grepl(\".qs\", submission_file)) {\n    \n    submit <- qread(file = submission_file)\n    \n} else if (grepl(\".tsv\", submission_file)) {\n    \n    submit <- fread(file = submission_file, sep = \"\\t\")\n}\n\n\ncat(\"The head of the data\")\nhead(submit)","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-08-23T21:39:35.37451Z","iopub.execute_input":"2023-08-23T21:39:35.375594Z","iopub.status.idle":"2023-08-23T21:39:55.889321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Set column names**","metadata":{}},{"cell_type":"code","source":"if (\"pred_cor\" %in% names(submit)) {\n    submit[, pred := pred_cor]\n    submit[, pred_cor := NULL]\n}\nif (\"aspect\" %in% names(submit)) {\n    submit[, aspect := NULL]\n}\n\nnames(submit) <- c(\"EntryID\", \"term\", \"pred\")","metadata":{"execution":{"iopub.status.busy":"2023-08-23T21:39:55.89321Z","iopub.execute_input":"2023-08-23T21:39:55.894669Z","iopub.status.idle":"2023-08-23T21:39:56.625354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Distributions of GO Terms and Proteins","metadata":{}},{"cell_type":"code","source":"compute_stats(submit)","metadata":{"execution":{"iopub.status.busy":"2023-08-23T21:39:56.628479Z","iopub.execute_input":"2023-08-23T21:39:56.629687Z","iopub.status.idle":"2023-08-23T21:40:08.039887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Add aspect**","metadata":{}},{"cell_type":"code","source":"submit[terms, aspect := i.aspect, on = \"term\"]\nsubmit <- submit[!is.na(aspect)]\nhead(submit)","metadata":{"execution":{"iopub.status.busy":"2023-08-23T21:40:08.041791Z","iopub.execute_input":"2023-08-23T21:40:08.042822Z","iopub.status.idle":"2023-08-23T21:40:18.322222Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Filter top N terms per ontology**","metadata":{}},{"cell_type":"code","source":"top_n = 30\n\nsubmit_filt <- submit[, .SD[order(-pred)][1:min(.N, top_n)], by = .(EntryID, aspect)]","metadata":{"execution":{"iopub.status.busy":"2023-08-23T21:40:18.324009Z","iopub.execute_input":"2023-08-23T21:40:18.324987Z","iopub.status.idle":"2023-08-23T21:42:55.216113Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"compute_stats(submit_filt)","metadata":{"execution":{"iopub.status.busy":"2023-08-23T21:42:55.218166Z","iopub.execute_input":"2023-08-23T21:42:55.219271Z","iopub.status.idle":"2023-08-23T21:42:57.261529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predicted terms of a single protein in the GO DAG","metadata":{}},{"cell_type":"code","source":"# root terms have level = 0\nget_go_level <- function(term, l = 1,\n                         roots = c(\"GO:0008150\",\n                                   \"GO:0005575\",\n                                   \"GO:0003674\")) {    \n  if (term %in% roots) {      \n      return(0)\n  }\n  \n  parents <- ontology$parents[[term]]\n  parents <- parents[!parents %in% roots]\n  \n  if (length(parents) == 0 ) {    \n    return(l)\n  }\n  \n  max_level <- l\n  \n  for (parent in parents) {\n    level <- get_go_level(parent, l = l + 1)\n    max_level <- max(max_level, level)\n  }  \n  return(max_level)\n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-08-23T21:42:57.263543Z","iopub.execute_input":"2023-08-23T21:42:57.264685Z","iopub.status.idle":"2023-08-23T21:42:57.275009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<div style=\"line-height:24px; font-size:16px;border-left: 5px solid silver; padding-left: 26px;\"> \n    📌 Task: \n    <div style=\"padding-left: 16px;\">\n     Visualise GO terms of a single protein of interest in the form of a directed acyclic graph (DAG) - aka ancestor chart (e.g. <a href = \"https://www.ebi.ac.uk/QuickGO/term/GO:0010977\">QuickGO - for terms </a>)\n    </div>\n         Please see the notebook describing GO structure <a href=\"https://www.kaggle.com/antoninadolgorukova/gene-ontology-explorer\">here</a>.\n</div>","metadata":{}},{"cell_type":"code","source":"prep_parents_graph <- function(dt, protein, terms = \"all\") {\n    \n    edges <- dt[EntryID == protein, .(term, pred)] \n    \n    if (terms == \"all\") {\n        \n        terms_of_interest <- unique(edges$term)\n    } else {\n        \n        terms_of_interest <- terms\n        \n    }\n    ancestors <- data.table(\n        term = unique(unlist(\n            lapply(terms_of_interest, function(x) ontology$ancestors[[x]])\n                   )),\n        pred = 0)    \n        \n    edges <- edges[term %in% c(terms_of_interest, ancestors$term)]    \n        \n    pred_terms <- unique(edges$term)\n    edges <- rbind(edges, ancestors[!term %in% edges$term])    \n\n\n    edges[, parents := lapply(term, function(x) ontology$parents[[x]])]\n    edges <- edges[\n      , .(from = term,\n          to = as.character(unlist(parents))),\n      by = .(term, pred)] # root have no parents and get NA\n    \n    edges <- edges[, .(pred, from, to)]\n\n    nodes <- data.table(\n        id = unique(edges$from),\n        label = sapply(unique(edges$from), function(x) ontology$name[[x]]),\n        level = sapply(unique(edges$from), function(x) get_go_level(x)),\n        shape = \"box\",\n        color.background = ifelse(unique(edges$from) %in% pred_terms, \"white\", \"beige\"), \n        color.border = \"black\",\n        font.color = \"black\"\n                      )\n    nodes[edges, pred := i.pred, on = c(\"id\" = \"from\")]                 \n    nodes[, label := str_wrap(paste0(pred, \"\\n\", id, \":\\n\", label),\n                              width = 10)]                \n\n    out <- list(nodes, edges)\n    names(out) <- c(\"nodes\", \"edges\")\n    return(out)\n}            ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-08-23T21:42:57.278384Z","iopub.execute_input":"2023-08-23T21:42:57.279624Z","iopub.status.idle":"2023-08-23T21:42:57.290286Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"protein = \"A0A023PXC2\" ","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-08-23T21:42:57.292382Z","iopub.execute_input":"2023-08-23T21:42:57.29344Z","iopub.status.idle":"2023-08-23T21:42:57.304293Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"graph <- prep_parents_graph(submit_filt, protein)\nplot <- visNetwork(graph$nodes, graph$edges, width=\"100%\", #height = \"350%\", width = \"2000px\", height = \"700px\"\n           main=paste0(\"The location of the protein of interest GO terms in the GO-DAG\")) %>%\n    visOptions(highlightNearest = list(enabled = TRUE, algorithm = \"hierarchical\"), selectedBy = \"label\") %>%\n    #visNodes(physics = FALSE) %>%\n    visHierarchicalLayout(direction = \"UD\")  #%>%, treeSpacing = 200\n    #visLegend(addEdges = dom_graph$leg_edges) #%>%\n    #visConfigure(enabled = TRUE)\nplot\nvisSave(plot, paste0(\"CAFA5_\", protein, \"_terms_ancestor_chart.html\"))","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-08-23T21:42:57.306503Z","iopub.execute_input":"2023-08-23T21:42:57.307575Z","iopub.status.idle":"2023-08-23T21:42:59.24042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<div style=\"line-height:24px; font-size:16px;border-left: 5px solid silver; padding-left: 26px;\"> \n    💡 Note: \n    <ul style=\"list-style:circle\">\n<li>Here you see top 30 predicted GO terms of the selected protein\n<li>The added parent terms (not present in the submission file) are displayed in colored boxes. They will get the maximum score of their child terms after the popagation implemented in the submission processing\n    </ul>\n</div>","metadata":{}},{"cell_type":"markdown","source":"# Track individual branches","metadata":{}},{"cell_type":"code","source":"graph <- prep_parents_graph(submit, protein, terms = \"GO:0031323\")\nplot <- visNetwork(graph$nodes, graph$edges, width=\"100%\", #height = \"350%\", width = \"2000px\", height = \"700px\"\n           main=paste0(\"Individual branches of GO terms in the GO-DAG\")) %>%\n    visOptions(highlightNearest = list(enabled = TRUE, algorithm = \"hierarchical\"), selectedBy = \"label\") %>%\n    #visNodes(physics = FALSE) %>%\n    visHierarchicalLayout(direction = \"UD\")  %>%#, treeSpacing = 200\n    #visLegend(addEdges = dom_graph$leg_edges) #%>%\n    #visConfigure(enabled = TRUE) %>%\n    visPhysics(solver = \"hierarchicalRepulsion\")\nplot","metadata":{"execution":{"iopub.status.busy":"2023-08-23T21:42:59.242491Z","iopub.execute_input":"2023-08-23T21:42:59.243546Z","iopub.status.idle":"2023-08-23T21:43:00.470014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]}]}