---
title: 'CAFA 6 Protein Data Host species and GO Terms'
author: 'Sam Lycett'
date: '`r Sys.Date()`'
output:
  html_document:
    number_sections: true
    toc: true
---

```{r, echo=FALSE}
library(knitr)
knitr::opts_chunk$set(echo = FALSE)
```

## Summary

Competition data from: 
<https://www.kaggle.com/competitions/cafa-6-protein-function-prediction>

The purpose of this R markdown is to explore the data from a biological perspective and calculate some statistics.

Please see my R notebook for details of how additional files were generated, and the additional files are now uploaded as a separate dataset.

<ul>
<li>R notebook:        <https://www.kaggle.com/code/redberry/cafa6-protein-data-shape>
<li>Additional files:  <https://www.kaggle.com/datasets/redberry/additional-data-for-cafa6-protein-data-shape>
</ul>

The result from <https://www.kaggle.com/code/redberry/cafa6-protein-data-shape> is file 
train_terms_with_info.csv, and this can now be used in CAFA 6 Protein Function Prediction models.

Note - 14 Dec 2025: I have now re-run this notebook to take in the new slightly updated CAFA 6 files (hopefully), 
and the subsequent re-run and re-uploaded additional data files.


```{r, echo=FALSE}
# Set up libraries
outPath <- "/kaggle/working/"
lib     <- paste0(outPath,"R-site-library/")

# do this once
if (!dir.exists(lib)) {
    dir.create(lib)
    install.packages("seqinr",lib=lib)
    install.packages("ontologyIndex",lib=lib)
    install.packages("ontologySimilarity",lib=lib)    
}
```

Specialist R packages are required to read the protein sequence fasta files (seqinr), 
and to read the GO terms ontology files (ontologyIndex, ontologySimilarity).
These are downloaded into `r lib` associated with this R markdown.

```{r, echo=TRUE}
# load specialist packages
library(seqinr,lib=lib)
library(ontologyIndex,lib=lib)
library(ontologySimilarity,lib=lib)
```


## Input Data files

```{r, echo=FALSE}
rootPath  <- "../input/"
dataPath  <- paste0(rootPath,"cafa-6-protein-function-prediction/")
trainPath <- paste0(dataPath,"Train/")
testPath  <- paste0(dataPath,"Test/")
train_fnames <- dir(trainPath)
test_fnames  <- dir(testPath)
other_fnames <- dir(dataPath)
other_fnames <- other_fnames[grep("tsv$",other_fnames)]

additionalPath <- paste0(rootPath,"additional-data-for-cafa6-protein-data-shape/")
```

The data files in the competition dataset are:

In top directory <b>cafa-6-protein-function-prediction/</b>: <b>
`r paste(other_fnames, collapse=" , ")`</b>

In <b>Train</b> subdirectory: <b>
`r paste(train_fnames, collapse=" , ")`</b>

In <b>Test</b> subdirectory: <b>
`r paste(test_fnames, collapse=" , ")`</b>

The additional files in `r additionalPath` are: <b>
`r paste(list.files(path = additionalPath), collapse=" , ")`</b>

## Protein Train Terms

```{r}
train_terms <- read.csv(paste0(trainPath,"train_terms.tsv"),sep="\t")
uids_tbl    <- table(train_terms$EntryID)
uids        <- names(uids_tbl)
uterms_tbl  <- table(train_terms$term)
uterms      <- names(uterms_tbl)
```

The train_terms.tsv data file contains `r length(train_terms[,1])` entries, of which there are `r length(uids)` unique protein ids, `r length(uterms)` unique GO terms, and the first 5 entries are:

```{r, echo=FALSE}
kable(train_terms[1:5,],caption="First 5 entries in train_terms.tsv")
```

There are 3 main ontologies (aspects) for the GO terms:  
C = cellular compartment (CCO), 
F = molecular function (MFO), 
P = biological process (BPO).  
These are distributed across the terms as follows:

```{r, echo=FALSE}
kable(t(table(train_terms$aspect)),caption="Distribution of aspects in train_terms.tsv")
```

For each unique id protein, there are a maximum of `r max(uids_tbl)` GO terms and a minimum of `r min(uids_tbl)`.  

`r format(100*length(which(uids_tbl<=3))/length(uids_tbl),digits=4)`% of unique proteins have <=3 GO terms.

Only `r length(which(uids_tbl>=100))` unique proteins (number not percent) have >=100 GO terms.

Only `r length(which(uids_tbl>=200))` unique proteins (number not percent) have >=200 GO terms, these are `r names(which(uids_tbl>=200))`.

The distribution of number of GO terms per protein is:

```{r, echo=FALSE}
h <- hist(uids_tbl, breaks= (1:(max(uids_tbl)+1))-0.5, plot=FALSE )
plot(h$mids,log10(h$count), pch=1, cex=0.5,
     xlab="Number of GO terms per protein", 
     ylab="Log10(Number of proteins)",
     main="Distribution of GO terms per unique protein")
```


## Host species of proteins

The file train_taxonomy.tsv contains the protein ids and the host species id (taxon) of these proteins.

```{r}
train_taxonomy <- read.csv(paste0(trainPath,"train_taxonomy.tsv"),sep="\t",header=FALSE)
```

```{r, echo=FALSE}
tinds <- match(uids,train_taxonomy[,1])
if (any(!is.finite(tinds))) {
  print("Warning cannot find all of the unique protein ids")
} else {
  if (all(uids==train_taxonomy[tinds,1])) {
    uids_taxons <- train_taxonomy[tinds,2]
  } else {
    stop("Problem matching protein uids with train_taxonomy")
  }
  taxon_tbl <- table(train_taxonomy[,2])
}
```


There are `r length(train_taxonomy[,1])` entries for the unique proteins in train_taxonomy.tsv, and these match the list of unique proteins in train_terms.tsv.

The host species taxon code with the most entries has `r max(taxon_tbl)` proteins, and this is taxon id=`r names(which.max(taxon_tbl))`

Uniprot has a species list of taxon id and host organism names, this can be found at
<https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/docs/speclist.txt>

```{r}
# now process the speclist (see R markdown, but not shown in this html)
# this creates the following smaller semi-colon (;) separated table files
species_list_file <- "speclist.txt"
real_file    <- "real_organism_codes.csv"
virtual_file <- "virtual_organism_codes.csv"
codes_file   <- "type_codes.csv"
```

The species list file (`r species_list_file`) has already been downloaded and is now in the additional data set: `r additionalPath`.

The species list file has already been processed to smaller semi-colon (;) separated table files: 
`r real_file`, `r virtual_file`, and `r codes_file`.

These files can be found in 
<https://www.kaggle.com/datasets/redberry/additional-data-for-cafa6-protein-data-shape>

```{r, echo=TRUE}
# read in processed additional data from speclist
real_taxon    <- read.csv(paste0(additionalPath,real_file),sep=";")
virtual_taxon <- read.csv(paste0(additionalPath,virtual_file),sep=";")
type_codes    <- read.csv(paste0(additionalPath,codes_file),sep=";")
```

```{r, echo=FALSE}
kable(type_codes,caption="Kingdom/Organism Type Codes")
```

```{r, echo=FALSE}
kable(real_taxon[grep(9606,real_taxon[,3])[1:4],],caption="Example 4 lines of real_taxon, these are specific species")
```

```{r, echo=FALSE}
kable(virtual_taxon[c(24,56,64,230),], caption="Example 4 lines of virtual_taxon, these are species groups")
```

```{r, echo=FALSE}
# real species full list
tbl <- table(factor(real_taxon$Type, levels=type_codes$Type))
tbl <- cbind(type_codes$Kingdom,tbl)

# by protein
colnames(train_taxonomy) <- c("Protein","Organism")
tinds <- match(train_taxonomy$Organism, real_taxon$Taxon)
if (any(is.na(tinds)) | !all(train_taxonomy$Organism==real_taxon$Taxon[tinds])) {
  stop("Cannot match all train_taxonomy organisms to real_taxon$Taxon")
}
train_org_types      <- factor(real_taxon$Type[tinds],levels=type_codes$Type)
tbl_train_orgs_types <- table(train_org_types)
train_host_names     <- real_taxon$Scientfic_Name[tinds]
train_taxonomy_with_host <- cbind(train_taxonomy,train_host_names)
colnames(train_taxonomy_with_host)[3] <- "Scientfic_Name"
top_train_host_names <- sort(table(train_host_names),decreasing = TRUE)

# by hosts in train
uhosts_in_train <- unique(train_taxonomy$Organism)
tinds           <- match(uhosts_in_train, real_taxon$Taxon)
uhosts_org_types<- factor(real_taxon$Type[tinds],levels=type_codes$Type)
tbl_uhosts_types <- table(uhosts_org_types)

tbl <- cbind(tbl,tbl_uhosts_types,tbl_train_orgs_types)
colnames(tbl) <- c("Kingdom","Number of real species","Number of species in train","Proteins per Kingdom in train")
```

In the train_taxomony, the top 10 species and the number of proteins in the training set are:

```{r, echo=FALSE}
kable(top_train_host_names[1:10],caption="Top 10 host species by number of proteins in train_taxonomy")
```

The distribution of the species (taxon) in the real full species list, as compared to the train_taxonomy is the following table.

From this we can see that there are no Viruses, only one Archea, and a few Bacteria in the training set.  This may be problematic for the future performance (real world) against unknown proteins, since the unknown proteins could come from a metagenomic study of environmental samples which may contain archea, bacteria and possibly viruses (although depends on how the proteins/DNA/RNA are extracted), and there is alot of metagenomic 'dark matter'.

```{r, echo=FALSE}
kable(tbl, caption="Number of species taxon in the real species list compared to the training set")
```

## Host species in test superset

```{r}
test_taxonomy <- read.csv(paste0(testPath,"testsuperset-taxon-list.tsv"),sep="\t")
```

The test superset taxon list contains `r length(test_taxonomy$ID)` entries, and `r length(unique(test_taxonomy$Species))` unique Species; the first 5 are:

```{r, echo=FALSE}
kable(test_taxonomy[1:5,], caption="First 5 entries of test taxonomy")
```

```{r, echo=FALSE}
match_train_organism_to_test_superset <- match(unique(train_taxonomy$Organism),unique(test_taxonomy$ID))
num_train_organism_in_test_superset   <- length(which(is.finite(match_train_organism_to_test_superset)))  # 1381
num_train_organism_not_in_test        <- length(which(!is.finite(match_train_organism_to_test_superset))) # 0
match_test_superset_to_train_organism <- match(unique(test_taxonomy$ID), unique(train_taxonomy$Organism))
num_test_superset_in_train_organism   <- length(which(is.finite(match_test_superset_to_train_organism)))  # 1381
num_test_not_in_train_organism        <- length(which(!is.finite(match_test_superset_to_train_organism))) # 7072
```

The number of unique Organisms in the training data (proteins) which match those in the test superset are `r num_train_organism_in_test_superset`, and the number of unique Organisms (from train) which are not in test superset are `r num_train_organism_not_in_test`.  There are `r num_test_not_in_train_organism` unique listed Species (Organisms) in the test superset which are not in the training data.

## Host species in test proteins

```{r, echo=FALSE}
# read in details of test set (dont actually need the proteins at this point)
test_taxonomy_with_host <- read.csv(paste0(additionalPath,"test_taxonomy_with_host.csv"),
                                   sep=";")
num_test_seqs <- length(test_taxonomy_with_host[,1])

# the unique hosts in the test proteins are listed in the test superset
# all(unique(test_taxonomy_with_host$Organism)==unique(test_taxonomy$ID))
```

There are `r num_test_seqs` protein sequences in the test data set; these have accession numbers and also the id of their host organisms in the sequence names.  The unique host organisms are listed in the test superset (see above).

```{r, echo=FALSE}
# real species full list (repeat from above in Train)
tbl <- table(factor(real_taxon$Type, levels=type_codes$Type))
tbl <- cbind(type_codes$Kingdom,tbl)

# by protein of test
tinds <- match(test_taxonomy_with_host$Organism, real_taxon$Taxon)
if (any(is.na(tinds)) | !all(test_taxonomy_with_host$Organism==real_taxon$Taxon[tinds])) {
  stop("Cannot match all test_taxonomy_with_host organisms to real_taxon$Taxon")
}
test_org_types       <- factor(real_taxon$Type[tinds],levels=type_codes$Type)
tbl_test_orgs_types  <- table(test_org_types)
test_host_names      <- real_taxon$Scientfic_Name[tinds]

# check these match
xinds <- which(test_host_names!=test_taxonomy_with_host$Scientfic_Name)
# there are some full stops differences, but otherwise OK

top_test_host_names <- sort(table(test_host_names),decreasing = TRUE)

# by hosts in test
uhosts_in_test   <- unique(test_taxonomy_with_host$Organism)
tinds            <- match(uhosts_in_test, real_taxon$Taxon)
uhosts_org_types <- factor(real_taxon$Type[tinds],levels=type_codes$Type)
tbl_uhosts_types <- table(uhosts_org_types)

tbl <- cbind(tbl,tbl_uhosts_types,tbl_test_orgs_types)
colnames(tbl) <- c("Kingdom","Number of real species","Number of species in test","Proteins per Kingdom in test")

```

In the test proteins data set, the top 10 species and the number of proteins in the test set are:

```{r, echo=FALSE}
kable(top_test_host_names[1:10],caption="Top 10 host species by number of proteins in test proteins dataset")
```

The distribution of the species (taxon) in the real full species list, as compared to the test data set is the following table.

From this we can see that there are no Viruses, only one Archea, and a few Bacteria in the test set.

```{r, echo=FALSE}
kable(tbl, caption="Number of species taxon in the real species list compared to the test set")
```

```{r, echo=FALSE}
# repeat of the test superset analysis as above
uspecies_in_test     <- sort(unique(test_taxonomy_with_host$Organism))
uspecies_in_train    <- sort(unique(train_taxonomy_with_host$Organism))
common_species       <- intersect(uspecies_in_test,uspecies_in_train)   # 1381
species_in_test_only <- setdiff(uspecies_in_test,uspecies_in_train)     # 0
species_in_train_only<- setdiff(uspecies_in_train,uspecies_in_test)     # 7072

# what are the organism types
tinds_common        <- match(common_species, real_taxon$Taxon)
common_types        <- factor(real_taxon$Type[tinds_common],levels=type_codes$Type)
tbl_common_types    <- table(common_types)
tinds_test          <- match(species_in_test_only, real_taxon$Taxon)
test_only_types     <- factor(real_taxon$Type[tinds_test],levels=type_codes$Type)
tbl_test_only_types <- table(test_only_types)

# real species full list (repeat from above in Train)
tbl <- table(factor(real_taxon$Type, levels=type_codes$Type))
tbl <- cbind(type_codes$Kingdom,tbl,tbl_common_types,tbl_test_only_types)

colnames(tbl) <- c("Kingdom","Number of real species","Number of species common to train and test","Number of species in test only")

kable(tbl, caption="Number of species in both the train and test sets compared to only the test set")
```

## GO Terms

The Gene Ontology terms (GO Terms) are within file go-basic.obo, this is a text based format, and can be read by R package ontologyIndex <https://cran.r-project.org/web/packages/ontologyIndex/> from paper Greene et al 2017 <https://doi.org/10.1093/bioinformatics/btw763>.

There are 3 root terms in the ontology:
<ul>
<li>GO:0008150	biological_process
<li>GO:0005575	cellular_component
<li>GO:0003674	molecular_function
</ul>


From the ontologyIndex R package description ("Introduction to ontologyX"), an ontology_index object (within R) is vector of lists, indexed by the IDs of the terms, and it has the following elements:
<ul>
<li>id
<li>name
<li>parents   (list of GO terms)
<li>children  (list of GO terms)
<li>ancestors (list of GO terms)
<li>obsolete  (TRUE/FALSE)
</ul>

The numbers of direct parents, children and ancestors of a term can be large.  The ontologySimilarity R package has a useful function, descendants_IC, which "calculates information content of terms based on frequency with which it is an ancestor of other terms", and this is a useful measure of the "depth" (detailed-ness) of a term from the roots.

```{r, echo=TRUE}
# read ontology file
ontology <- get_ontology(paste0(trainPath,"go-basic.obo"))

# Get information content based on number of descendants each term has
ontology_info <- descendants_IC(ontology)

# define root terms for ontology
bpo_term <- "GO:0008150"
cco_term <- "GO:0005575"
mfo_term <- "GO:0003674"
```

```{r, echo=FALSE}
# define custom function for ease of printing go term information
go_term_info_for_print <- function(gt, ontology, ontology_info, limitTerms=2) {
  if (limitTerms>0) {
    parTxt <- paste(ontology$parents[[gt]],  collapse=" ; ")
    numPar <- length(ontology$parents[[gt]])
    if (numPar > limitTerms) {
      parTxt <- paste(numPar,ontology$parents[[gt]][1],sep="+ ")
    }
    
    childTxt <- paste(ontology$children[[gt]],  collapse=" ; ")
    numChild <- length(ontology$children[[gt]])
    if (numChild > limitTerms) {
      childTxt <- paste(numChild,ontology$children[[gt]][1],sep="+ ")
    }
    
    ancTxt <- paste(ontology$ancestors[[gt]],collapse=" ; ")
    numAnc <- length(ontology$ancestors[[gt]])
    if (numAnc > limitTerms) {
      ancTxt <- paste(numAnc,ontology$ancestors[[gt]][1],sep="+ ")
    }
  }
  
  info_val_txt <- format(ontology_info[which(names(ontology_info)==gt)],digits=3) 
  
  gt_details <- c(ontology$id[[gt]],ontology$name[[gt]],
                  parTxt,
                  childTxt,
                  ancTxt,
                  ontology$obsolete[[gt]],
                  info_val_txt)
  names(gt_details) <- c("id","name","parents","children","ancestors","obsolete","descendants_IC")
  return(gt_details)
}

# print the details of the selected terms, and include the descends info
sel_terms  <- c(bpo_term,cco_term,mfo_term,
                train_terms$term[1],train_terms$term[2])
sel_info   <- matrix(0,length(sel_terms),7)
for (j in 1:length(sel_terms)) {
  res <- go_term_info_for_print(sel_terms[j],ontology,ontology_info)
  if (j==1) {
    colnames(sel_info) <- names(res)
    rownames(sel_info) <- sel_terms
  }
  sel_info[j,] <- res
}
```

Details of the root terms, and the first two GO terms in the list:

```{r, echo=FALSE}
kable(sel_info, caption="GO term information for 5 selected terms", row.names=FALSE)
```

```{r, echo=FALSE}
######################################################
# make protein train terms with attached information
######################################################

# combine information
# do this once
if (!file.exists(paste0(additionalPath,"train_terms_with_info.csv"))) {
  # matching terms into ontology_info for the descendants_IC
  tinds1                 <- match(train_terms$term,names(ontology_info))
  
  # matching proteins in train terms to proteins in hosts
  tinds2                 <- match(train_terms$EntryID,train_taxonomy_with_host$Protein)
  
  # combining with train_terms
  train_terms_with_info <- cbind(train_terms,
                                 ontology_info[tinds1],
                                 ontology$name[train_terms$term],
                                 train_taxonomy_with_host$Organism[tinds2],
                                 train_taxonomy_with_host$Scientfic_Name[tinds2])
  colnames(train_terms_with_info)[4] <- "term_info"
  colnames(train_terms_with_info)[5] <- "term_name"
  colnames(train_terms_with_info)[6] <- "Organism"
  colnames(train_terms_with_info)[7] <- "Scientific_Name"
  train_terms_with_info <- as.data.frame(train_terms_with_info)
  
  tinds3 <- match(train_terms_with_info$Organism,real_taxon$Taxon)
  train_terms_with_info <- cbind(train_terms_with_info,real_taxon$Type[tinds3])
  colnames(train_terms_with_info)[8] <- "Type"
  write.table(train_terms_with_info, file=paste0(additionalPath,"train_terms_with_info.csv"),
              row.names=FALSE, col.names=TRUE, sep=";")
}
```

The GO terms from the training data have been linked with 
the term descendents information and Organisation scientific name, 
and these have been written as a semi colon (;) separated file to train_terms_with_info.csv in the additional data files.

```{r, echo=TRUE}
train_terms_with_info <- read.csv(paste0(additionalPath,"train_terms_with_info.csv"),sep=";")
```

The first 5 terms are:

```{r}
kable(train_terms_with_info[1:5,], caption="Protein training terms with additional information")
```

The file train_terms_with_info.csv can now be used in CAFA 6 Protein Function Prediction models.