{"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":"code","source":"# This R environment comes with many helpful analytics packages installed\n# It is defined by the kaggle/rstats Docker image: https://github.com/kaggle/docker-rstats\n# For example, here's a helpful package to load\n\nlibrary(tidyverse) # metapackage of all tidyverse packages\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nlist.files(path = \"../input\")\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2023-05-28T04:50:37.799321Z","iopub.execute_input":"2023-05-28T04:50:37.801527Z","iopub.status.idle":"2023-05-28T04:50:39.146188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Source \n  - https://bioconductor.org/packages/release/bioc/html/topGO.html\n  - https://bioconductor.org/packages/devel/bioc/vignettes/topGO/inst/doc/topGO.pdf\n  - https://www.rdocumentation.org/packages/ontologyIndex/versions/2.10/topics/get_ontology\n  - https://rdrr.io/bioc/topGO/f/inst/doc/topGO.pdf\n  - https://bioconductor.org/packages/release/bioc/html/Biostrings.html\n  - https://stackoverflow.com/questions/21263636/read-fasta-into-a-dataframe-and-extract-subsequences-of-fasta-file\n  - https://mq-software-carpentry.github.io/r-ggplot-extension/02-categorical-data/index.html\n  - https://cran.r-project.org/web/packages/kableExtra/vignettes/awesome_table_in_html.html\n  - https://www.bioconductor.org/packages/release/data/experiment/vignettes/RforProteomics/inst/doc/RProtVis.html\n  - https://www.rdocumentation.org/packages/influential/versions/1.1.2/topics/graph_from_data_frame\n  - https://a-little-book-of-r-for-bioinformatics.readthedocs.io/en/latest/src/chapter11.html\n  - https://cran.r-project.org/web/packages/esquisse/vignettes/get-started.html\n  - https://www.bioconductor.org/packages/release/bioc/html/GOfuncR.html\n  - https://github.com/sgrote/GOfuncR","metadata":{}},{"cell_type":"markdown","source":"### Goal of the Competition\n- The goal of this competition is to predict the function of a set of proteins. You will develop a model trained on the amino-acid sequences of the proteins and on other data. Your work will help ​​researchers better understand the function of proteins, which is important for discovering how cells, tissues, and organs work. This may also aid in the development of new drugs and therapies for various diseases.\n\n### Context\n- Proteins are responsible for many activities in our tissues, organs, and bodies and they also play a central role in the structure and function of cells. Proteins are large molecules composed of 20 types of building-blocks known as amino acids. The human body makes tens of thousands of different proteins, and each protein is composed of dozens or hundreds of amino acids that are linked sequentially. This amino-acid sequence determines the 3D structure and conformational dynamics of the protein, and that, in turn, determines its biological function. \n\n### Background\n- The Gene Ontology (GO) is a concept hierarchy that describes the biological function of genes and gene products at different levels of abstraction (Ashburner et al., 2000). It is a good model to describe the multi-faceted nature of protein function.\n- GO is a directed acyclic graph. The nodes in this graph are functional descriptors (terms or classes) connected by relational ties between them (is_a, part_of, etc.). For example, terms 'protein binding activity' and 'binding activity' are related by an is_a relationship; however, the edge in the graph is often reversed to point from binding towards protein binding. This graph contains three subgraphs (subontologies): Molecular Function (MF), Biological Process (BP), and Cellular Component (CC), defined by their root nodes. Biologically, each subgraph represent a different aspect of the protein's function: what it does on a molecular level (MF), which biological processes it participates in (BP) and where in the cell it is located (CC). \n- We will use experimentally determined term-protein assignments as class labels for each protein. That is, if a protein is labeled with a term, it means that this protein has this function validated by experimental evidence. By processing these annotated terms, we can generate a dataset of proteins and their ground truth labels for each term. The absence of a term annotation does not necessarily mean a protein does not have this function, only that this annotation does not exist (yet) in the GO. A protein may be annotated by one or more terms from the same subontology, and by terms from more than one subontology.\n\n### Training Set\n- For the training set, we include all proteins with annotated terms that have been validated by experimental or high-throughput evidence, traceable author statement (evidence code TAS), or inferred by curator (IC). \n\n### Test Superset\n- The test superset is a set of protein sequences on which the participants are asked to predict GO terms.\n\n### Test Set\n- The test set is unknown at the beginning of the competition. It will contain protein sequences (and their functions) from the test superset that gained experimental annotations between the submission deadline and the time of evaluation.\n### File Descriptions\n- Gene Ontology: The ontology data is in the file go-basic.obo. This structure is the 2023-01-01 release of the GO graph. This file is in OBO format, for which there exist many parsing libraries.","metadata":{}},{"cell_type":"markdown","source":"### topGO\n- Enrichment Analysis for Gene Ontology\n- topGO package provides tools for testing GO terms while accounting for the topology of the GO graph. Different test statistics and different methods for eliminating local similarities and dependencies between GO terms can be implemented and applied.\n","metadata":{}},{"cell_type":"code","source":"if (!require(\"BiocManager\", quietly = TRUE))\n    install.packages(\"BiocManager\")\n\nBiocManager::install(\"topGO\")","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-05-28T04:50:39.150189Z","iopub.execute_input":"2023-05-28T04:50:39.184508Z","iopub.status.idle":"2023-05-28T04:55:03.893337Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"browseVignettes(\"topGO\")\nlibrary(topGO)","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-05-28T04:55:03.897334Z","iopub.execute_input":"2023-05-28T04:55:03.898929Z","iopub.status.idle":"2023-05-28T04:55:07.812342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"help(topGO)","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:55:07.815773Z","iopub.execute_input":"2023-05-28T04:55:07.818259Z","iopub.status.idle":"2023-05-28T04:55:08.04092Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"install.packages(\"ontologyIndex\")","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:55:08.044247Z","iopub.execute_input":"2023-05-28T04:55:08.046051Z","iopub.status.idle":"2023-05-28T04:55:17.539024Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"install.packages(\"influential\")","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:55:17.542155Z","iopub.execute_input":"2023-05-28T04:55:17.543854Z","iopub.status.idle":"2023-05-28T04:55:29.903504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"library(ontologyIndex)","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:55:29.906411Z","iopub.execute_input":"2023-05-28T04:55:29.908031Z","iopub.status.idle":"2023-05-28T04:55:29.945423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"file <- \"/kaggle/input/cafa-5-protein-function-prediction/Train/go-basic.obo\"\nget_OBO(\n  file,\n  propagate_relationships = \"is_a\",\n  extract_tags = \"minimal\",\n  merge_equivalent_terms = TRUE\n)","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:55:29.948293Z","iopub.execute_input":"2023-05-28T04:55:29.949861Z","iopub.status.idle":"2023-05-28T04:55:46.911043Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"library(influential)\nfile <- \"/kaggle/input/cafa-5-protein-function-prediction/Train/go-basic.obo\"\nobo_data <- get_OBO(\n  file,\n  propagate_relationships = \"is_a\",\n  extract_tags = \"minimal\",\n  merge_equivalent_terms = TRUE\n)\n\ngraph <- graph_from_data_frame(d=obo_data)\ngraph","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:55:46.915153Z","iopub.execute_input":"2023-05-28T04:55:46.916874Z","iopub.status.idle":"2023-05-28T04:56:02.92244Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"get_ontology(\n  file,\n  propagate_relationships = \"is_a\",\n  extract_tags = \"minimal\",\n  merge_equivalent_terms = TRUE\n)","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:56:02.926547Z","iopub.execute_input":"2023-05-28T04:56:02.92823Z","iopub.status.idle":"2023-05-28T04:56:17.090159Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if (!require(\"BiocManager\", quietly = TRUE))\n    install.packages(\"BiocManager\")\n\nBiocManager::install(\"Biostrings\")","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-05-28T04:56:17.094142Z","iopub.execute_input":"2023-05-28T04:56:17.095837Z","iopub.status.idle":"2023-05-28T04:57:54.364548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"library(\"Biostrings\")\n\nfastaFile <- readDNAStringSet(\"/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta\")\nseq_name = names(fastaFile)\nsequence = paste(fastaFile)\ndf <- data.frame(seq_name, sequence)","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:57:54.367289Z","iopub.execute_input":"2023-05-28T04:57:54.368778Z","iopub.status.idle":"2023-05-28T04:57:58.027377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(df,3)","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:57:58.030137Z","iopub.execute_input":"2023-05-28T04:57:58.031587Z","iopub.status.idle":"2023-05-28T04:57:58.058533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ggplot(data = df, aes(x=sequence)) +\n#     geom_bar(position = \"dodge\")","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:57:58.061206Z","iopub.execute_input":"2023-05-28T04:57:58.062652Z","iopub.status.idle":"2023-05-28T04:57:58.073994Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ggplot(df, aes(x=sequence)) + geom_boxplot(fill='steelblue')\n","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:57:58.07669Z","iopub.execute_input":"2023-05-28T04:57:58.078167Z","iopub.status.idle":"2023-05-28T04:57:58.089219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"require(data.table)\n\ntaxonomy_data<-as.data.frame(fread(\"/kaggle/input/cafa-5-protein-function-prediction/Train/train_taxonomy.tsv\"))\nhead(taxonomy_data,3)","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:57:58.091995Z","iopub.execute_input":"2023-05-28T04:57:58.093443Z","iopub.status.idle":"2023-05-28T04:57:58.465217Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"graph <- graph_from_data_frame(d=taxonomy_data)\ngraph\n","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:57:58.469328Z","iopub.execute_input":"2023-05-28T04:57:58.471033Z","iopub.status.idle":"2023-05-28T04:57:59.809794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"library(readr)\ntrain_terms <- read_tsv('/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv', col_names=FALSE)\nhead(train_terms,3)","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:57:59.813116Z","iopub.execute_input":"2023-05-28T04:57:59.814628Z","iopub.status.idle":"2023-05-28T04:58:04.217052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"graph <- graph_from_data_frame(d=train_terms)\ngraph","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:58:04.221295Z","iopub.execute_input":"2023-05-28T04:58:04.223387Z","iopub.status.idle":"2023-05-28T04:58:23.121407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# devtools::install_github(\"alesmascaro/BCDAG\")\n# library(BCDAG)\n# help(BCDAG)","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:58:23.124399Z","iopub.execute_input":"2023-05-28T04:58:23.12586Z","iopub.status.idle":"2023-05-28T04:58:23.13695Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# library(esquisse)\n# esquisse::esquisser(train_terms)","metadata":{"execution":{"iopub.status.busy":"2023-05-28T04:58:23.139758Z","iopub.execute_input":"2023-05-28T04:58:23.14118Z","iopub.status.idle":"2023-05-28T04:58:23.152022Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### GOfuncR \n - GOfuncR performs a gene ontology enrichment analysis based on the ontology enrichment software","metadata":{}},{"cell_type":"code","source":"if (!require(\"BiocManager\", quietly = TRUE))\n    install.packages(\"BiocManager\")\n\nBiocManager::install(\"GOfuncR\")","metadata":{"execution":{"iopub.status.busy":"2023-05-28T05:01:36.992139Z","iopub.execute_input":"2023-05-28T05:01:36.994095Z","iopub.status.idle":"2023-05-28T05:05:41.947945Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ?go_enrich","metadata":{"execution":{"iopub.status.busy":"2023-05-28T05:41:39.159905Z","iopub.execute_input":"2023-05-28T05:41:39.161796Z","iopub.status.idle":"2023-05-28T05:41:39.174379Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Directed acyclic graph \n\n# library(igraph)\n# library(GOfuncR)\n# # Graph object\n# graph <- graph.empty()\n\n# # Graph Nodes\n# for (i in 1:nrow(df)) {\n#   graph <- add.vertices(graph, 1, name = df$seq_name[i])\n# }\n\n# # Graph Edges\n# for (i in 1:nrow(obo_data)) {\n#   for (j in 1:nrow(obo_data)) {\n#     if (df$is_a[i] == df$seq_name[j]) {\n#       graph <- add.edges(graph, i, j)\n#     }\n#   }\n# }\n\n# plot(graph)\n","metadata":{"execution":{"iopub.status.busy":"2023-05-28T05:24:30.431117Z","iopub.execute_input":"2023-05-28T05:24:30.432835Z","iopub.status.idle":"2023-05-28T05:31:31.070161Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"###  biomaRt\n - biomaRt provides an interface to a growing collection of databases implementing the BioMart software suite ()","metadata":{}},{"cell_type":"code","source":"if (!require(\"BiocManager\", quietly = TRUE))\n    install.packages(\"BiocManager\")\n\nBiocManager::install(\"biomaRt\")","metadata":{"execution":{"iopub.status.busy":"2023-05-28T05:32:17.334973Z","iopub.execute_input":"2023-05-28T05:32:17.338482Z","iopub.status.idle":"2023-05-28T05:33:05.36929Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ?listDatasets()","metadata":{"execution":{"iopub.status.busy":"2023-05-28T06:21:47.654231Z","iopub.execute_input":"2023-05-28T06:21:47.656227Z","iopub.status.idle":"2023-05-28T06:21:47.670623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"library(biomaRt)\nlistEnsembl()\nensembl <- useEnsembl(biomart = \"ensembl\", dataset = 'hsapiens_gene_ensembl')\n## list the available datasets in this Mart\nlistDatasets(mart = ensembl)","metadata":{"execution":{"iopub.status.busy":"2023-05-28T06:40:00.336944Z","iopub.execute_input":"2023-05-28T06:40:00.339715Z","iopub.status.idle":"2023-05-28T06:40:13.261295Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"?listAttributes","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# data frame contaiNS the gene ID, gene symbol, and chromosome for all of the human genes.\nlibrary(biomaRt)\n\n# Connect to the Ensembl BioMart database\nensembl <- useEnsembl(biomart = \"genes\", dataset = \"hsapiens_gene_ensembl\")\n\n# Specify the attributes you want to retrieve\nattributes <- c(\"ensembl_gene_id\", \"hgnc_symbol\", \"chromosome_name\")\n\n# Retrieve the data\ndata <- getBM(attributes = attributes, mart = ensembl)\n\n# View the data\nhead(data)","metadata":{"execution":{"iopub.status.busy":"2023-05-28T05:43:55.867393Z","iopub.execute_input":"2023-05-28T05:43:55.872165Z","iopub.status.idle":"2023-05-28T05:44:14.095249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # Directed acyclic graph \n\n# library(igraph)\n# library(GOfuncR)\n# # Graph object\n# graph <- graph.empty()\n\n# # Graph Nodes\n# for (i in 1:nrow(data)) {\n#   graph <- add.vertices(graph, 1, name = data$ensembl_gene_id[i])\n# }\n\n# # Graph Edges\n# for (i in 1:nrow(data)) {\n#   for (j in 1:nrow(data)) {\n#     if (data$is_a[i] == data$ensembl_gene_id[j]) {\n#       graph <- add.edges(graph, i, j)\n#     }\n#   }\n# }\n\n# plot(graph)\n","metadata":{"execution":{"iopub.status.busy":"2023-05-28T05:46:47.837307Z","iopub.execute_input":"2023-05-28T05:46:47.843203Z","iopub.status.idle":"2023-05-28T05:47:36.803839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}