{"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":"# Exploratory analysis of train and test data","metadata":{}},{"cell_type":"markdown","source":"

What we do here:

\n
\n \n Useful notebook that inspired this one\n \n
","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle"}},{"cell_type":"markdown","source":"
\n\n

Data

\n\n
\n
    \n
  1. TSV file with CAFA5 training set. This set contains 142246 proteins with annotated GO terms (31466 in total) that have been validated by experimental or high-throughput evidence, traceable author statement (evidence code TAS), or inferred by curator (IC).\n
  2. fasta file with test dataset (test superset), containing 141864 proteins and their sequences\n
  3. TSV file with 90 taxons for proteins of the test superset\n
\n
\n

Output

\n \n
\n RDS/qs/csv file with test superset in a data.table format \n
\n
","metadata":{}},{"cell_type":"code","source":"if (!require(\"BiocManager\", quietly = TRUE))\n install.packages(\"BiocManager\")\nBiocManager::install(version = \"3.12\")\n\nsuppressPackageStartupMessages({ #takes around 5 min\n \n #library(stringr) #wrap text into paragraphs\n library(dplyr) #basic\n library(tidyr) #basic\n library(data.table) #work with data.table\n library(scales) #percent()\n library(tictoc) #measure processing time\n library(qs)\n \n BiocManager::install(\"Biostrings\")\n library(Biostrings) # read fasta\n \n # plots\n library(ggplot2)\n library(ggvenn) #venn diagrams with labels\n library(patchwork) #arrange plots\n})","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Useful functions**","metadata":{}},{"cell_type":"code","source":"#function for figure size adjusment\nfig <- function(width, heigth) {\n options(repr.plot.width = width, repr.plot.height = heigth)\n}\n\n# Output the number of unique values in each column\ndata.summary <- function(dt) {\n \n #print()\n cat(\"\\nClass(es): \", paste(class(dt), collapse = \", \"),\n \"\\nDimentions: \", ncol(dt), \" columns x \",\n nrow(dt), \" rows\", sep = \"\")\n # Calculate the number of unique values \n unique_counts <- data.table(\n column = names(dt),\n unique_count = lapply(dt, function(x) length(unique(x))))\n cat(\"\\nNumber of unique values in columns:\\n\\n\")\n print(unique_counts)\n}\n\n# Output descriptive statistics for one numeric column of a df\nstat_calcs <- function(df, colname) {\n \n out <- df %>%\n summarise_at(colname, .funs = list(\n \"Min\" = min, \n \"10%\" = function(x)quantile(x, .1), \n \"25%\" = function(x)quantile(x, .25), \n \"Median\" = median, \n \"Mean\" = mean, \n \"75%\" = function(x)quantile(x, .75), \n \"90%\" = function(x)quantile(x, .90), \n \"Max\" = max) ) %>% \n pivot_longer(cols = everything(), \n names_to = \"stat\", \n values_to = colname ) %>% \n as.data.frame\n out[, 2] <- round(out[, 2], 0)\n out$stat <- factor(out$stat, levels = out$stat)\n out\n}","metadata":{"execution":{"iopub.status.busy":"2023-05-25T18:48:31.594585Z","iopub.execute_input":"2023-05-25T18:48:31.597191Z","iopub.status.idle":"2023-05-25T18:48:31.619894Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data loading","metadata":{}},{"cell_type":"markdown","source":"**Loading the training dataset**","metadata":{}},{"cell_type":"code","source":"train_set_all <- readRDS(\"/kaggle/input/cafa5-supp-pre-calcs-for-ml/datasets/train_set_supp_GOdb3.15.rds\")\n\ncat(\"The head of the combined train data with additional info (parents, children, levels) added from the GO.db 3.15 database:\")\nhead(train_set_all) %>% mutate(sequence = paste0(substring(sequence, 1, 10), \"...\") )\ndata.summary(train_set_all)","metadata":{"execution":{"iopub.status.busy":"2023-05-25T18:48:34.106081Z","iopub.execute_input":"2023-05-25T18:48:34.107885Z","iopub.status.idle":"2023-05-25T18:49:06.007602Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Loading a fasta file with test data and taxon information**","metadata":{}},{"cell_type":"code","source":"test_fasta <- readAAStringSet(\"/kaggle/input/cafa-5-protein-function-prediction/Test (Targets)/testsuperset.fasta\")\n\n#cat(\"Test superset fasta file structure:\\n\")\n#str(test_fasta)\n\n# Convert to a data frame\nproteins <- gsub(\"\\t.*\", \"\", names(test_fasta))\ntaxons <- gsub(\".*\\t\", \"\", names(test_fasta))\nsequences <- paste(test_fasta)\n\ntest_set <- data.table(EntryID = proteins,\n taxonomyID = taxons, \n sequence = sequences)\n# Add taxon info\ntaxon_list <- read.table(file = \"/kaggle/input/cafa-5-protein-function-prediction/Test (Targets)/testsuperset-taxon-list.tsv\", sep = '\\t', header = TRUE)\n#data.summary(taxon_list)\n\n# Merge data\ntaxon_list$ID <- as.character(taxon_list$ID)\nnames(taxon_list)[names(taxon_list) == \"ID\"] <- \"taxonomyID\"\ntest_set <- data.table(test_set, key = \"taxonomyID\")\ntaxon_list <- data.table(taxon_list, key = \"taxonomyID\")\n\ntest_set[taxon_list, species := Species]\n\n# Output\ncat(\"The head of the test set extracted from the fasta file:\")\nhead(test_set) %>% mutate(sequence = paste0(substring(sequence, 1, 10), \"...\") )\n\ndata.summary(test_set)","metadata":{"execution":{"iopub.status.busy":"2023-05-25T18:51:26.278756Z","iopub.execute_input":"2023-05-25T18:51:26.280583Z","iopub.status.idle":"2023-05-25T18:51:27.732698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Sanity check**","metadata":{}},{"cell_type":"code","source":"# Remove duplicates\ndupl <- test_set[duplicated(test_set$EntryID), ]$EntryID\n\ncat(\"Duplicated protein in the test superset\")\ntest_set[test_set$EntryID %in% dupl, ]\n\ntest_set <- test_set[!duplicated(test_set$EntryID), ]","metadata":{"execution":{"iopub.status.busy":"2023-05-25T18:51:49.45194Z","iopub.execute_input":"2023-05-25T18:51:49.453696Z","iopub.status.idle":"2023-05-25T18:51:49.533831Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Additional calculations**","metadata":{}},{"cell_type":"code","source":"# Test set\ntest_set[, seq_length := nchar(sequence) ]\ntest_set[, prot_per_tax := .N, by = taxonomyID]\ntest_set[, prot_per_seq := .N, by = sequence]\n\n# Train set\n\ntrain_set_all[, prot_per_tax := uniqueN(EntryID), by = taxonomyID]\ntrain_set_all[, prot_per_seq := uniqueN(EntryID), by = sequence]\ntrain_set_all[, ont_per_prot := uniqueN(aspect), by = EntryID]\ntrain_set_all[, term_per_prot := uniqueN(term), by = EntryID]\n\n# Remove columns analysed in previously: https://www.kaggle.com/code/antoninadolgorukova/cafa5-train-set-go-terms-eda\ntrain_set <- copy(train_set_all)\ntrain_set <- train_set[, c('term', 'aspect', 'level', 'missing_parents', 'has_children'):= NULL] %>% unique\n\n# Output\ncat(\"The head of the resulting test set (n = \", nrow(test_set), \" proteins in rows)\")\nhead(test_set) %>% mutate(sequence = paste0(substring(sequence, 1, 10), \"...\") )\n\ncat(\"\\nThe head of the resulting training set (n = \", nrow(train_set), \" proteins in rows)\")\nhead(train_set) %>% mutate(sequence = paste0(substring(sequence, 1, 10), \"...\") )","metadata":{"execution":{"iopub.status.busy":"2023-05-25T18:51:52.867632Z","iopub.execute_input":"2023-05-25T18:51:52.869371Z","iopub.status.idle":"2023-05-25T18:52:09.769762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Store test data as RDS for future use**","metadata":{}},{"cell_type":"code","source":"#save RDS to for easy use later\nqsave(test_set, \"test_set.qs\")\nsaveRDS(test_set, \"test_set.rds\")\nfwrite(test_set, \"test_set.csv\", sep = \";\")","metadata":{"execution":{"iopub.status.busy":"2023-05-25T18:52:09.772783Z","iopub.execute_input":"2023-05-25T18:52:09.774305Z","iopub.status.idle":"2023-05-25T18:52:17.291517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # read csv\n# check_test <- fread(\"/kaggle/working/test_set.csv\")\n# check_test %>% head %>% mutate(sequence = paste0(substring(sequence, 1, 10), \"...\") )\n\n# # read qs\n# check_test <- qread(\"/kaggle/working/test_set.qs\")\n# check_test %>% head %>% mutate(sequence = paste0(substring(sequence, 1, 10), \"...\") )\n\n# # read RDS\n# check_test <- readRDS(\"/kaggle/working/test_set.rds\")\n# check_test %>% head %>% mutate(sequence = paste0(substring(sequence, 1, 10), \"...\") )","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Proteins and their sequences","metadata":{}},{"cell_type":"code","source":"cat(\"Number of protein IDs in the train set:\", n_distinct(train_set$EntryID))\ncat(\"\\nNumber of proteins IDs in the test set:\", n_distinct(test_set$EntryID))\n\ncat(\"\\n\\nNumber of unique sequences in the train set:\", n_distinct(train_set$sequence))\ncat(\"\\nNumber of unique sequences in the test set:\", n_distinct(test_set$sequence))","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:07:06.776863Z","iopub.execute_input":"2023-05-23T12:07:06.778253Z","iopub.status.idle":"2023-05-23T12:07:06.828637Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 10)\nvenn_diag <- list(\n \"Train set\" = unique(train_set$EntryID), \n \"Test set\" = unique(test_set$EntryID))\n\np1 <- ggvenn(venn_diag, show_elements = FALSE, label_sep = \"\\n\",\n text_size = 6, set_name_size = 8,\n fill_color = c(\"#cc7722\", \"#01796f\")) +\n ggtitle(paste(\"Intersection of proteins in the datasets\")) +\n theme(plot.title = element_text(size = 20, face = \"bold\", hjust = 0.5))\n\nvenn_diag <- list(\n \"Train set\" = unique(train_set$sequence), \n \"Test set\" = unique(test_set$sequence))\n\np2 <- ggvenn(venn_diag, show_elements = FALSE, label_sep = \"\\n\",\n text_size = 6, set_name_size = 8,\n fill_color = c(\"#cc7722\", \"#01796f\")) +\n ggtitle(paste(\"Intersection of protein sequences in the datasets\")) +\n theme(plot.title = element_text(size = 20, face = \"bold\", hjust = 0.5))\np1 + p2","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:07:09.976667Z","iopub.execute_input":"2023-05-23T12:07:09.978222Z","iopub.status.idle":"2023-05-23T12:07:13.450088Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Proteins with the same sequence and different name**","metadata":{}},{"cell_type":"code","source":"cat(\"The top of the subset with the same sequences for different protein names (n = \", n_distinct(train_set[prot_per_seq > 1]$sequence), \")\")\nsetorder(train_set[prot_per_seq > 1], sequence) %>% mutate(sequence = paste0(substring(sequence, 1, 50), \"...\") ) %>% head(7)","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:07:23.265334Z","iopub.execute_input":"2023-05-23T12:07:23.267251Z","iopub.status.idle":"2023-05-23T12:07:23.324618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dat <- rbind(train_set %>% select(\"sequence\", \"prot_per_seq\") %>% unique %>% mutate(dataset = \"Train\"),\n test_set %>% select(\"sequence\", \"prot_per_seq\") %>% unique %>% mutate(dataset = \"Test\"))\n\ncat(\"Distribution of the number of proteins per sequence\")\ndstats <- merge(stat_calcs(dat[dat$dataset == \"Train\", ], colname = \"prot_per_seq\") %>%\n dplyr::rename_with(.fn = ~paste0(\"Train (n = \",\n n_distinct(train_set$sequence),\n \")\"),\n .cols = prot_per_seq),\n stat_calcs(dat[dat$dataset == \"Test\", ], colname = \"prot_per_seq\") %>%\n dplyr::rename_with(.fn = ~paste0(\"Test (n = \",\n n_distinct(test_set$sequence),\n \")\"), .cols = prot_per_seq) ) %>%\n arrange(stat)\ndstats","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:07:27.610982Z","iopub.execute_input":"2023-05-23T12:07:27.612927Z","iopub.status.idle":"2023-05-23T12:07:27.904044Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#-------------train set--------------------#\ncat(\"\\033[1mTrain dataset\\033[0m\\n\")\n\ncat(\"\\nSequences with only 1 protein name:\", nrow(dat[dataset == \"Train\" & prot_per_seq == 1, ]), \"(\", \n percent(nrow(dat[dataset == \"Train\" & prot_per_seq == 1, ])/n_distinct(train_set$sequence), accuracy = 0.01),\")\" )\n\ncat(\"\\nSequences with 2 protein names or more:\",\n nrow(dat[dataset == \"Train\" & prot_per_seq >= 2, ]), \"(\", \n percent(nrow(dat[dataset == \"Train\" & prot_per_seq >= 2, ])/n_distinct(train_set$sequence), accuracy = 0.01),\")\")\n\nstat <- dstats[dstats$stat == \"Max\", \"Train (n = 138924)\"]\ncat(\"\\nSequences with \", stat,\n \"(max) protein names:\", nrow(dat[dataset == \"Train\" & prot_per_seq == stat, ]), \"(\", \n percent(nrow(dat[dataset == \"Train\" & prot_per_seq == stat, ])/n_distinct(train_set$sequence), accuracy = 0.01),\")\" )\n\n#-------------test set--------------------#\ncat(\"\\n\\n\\033[1mTest dataset\\033[0m\\n\")\n\ncat(\"\\nSequences with only 1 protein name:\", nrow(dat[dataset == \"Test\" & prot_per_seq == 1, ]), \"(\", \n percent(nrow(dat[dataset == \"Test\" & prot_per_seq == 1, ])/n_distinct(test_set$sequence), accuracy = 0.01),\")\" )\n\ncat(\"\\nSequences with 2 protein names or more:\",\n nrow(dat[dataset == \"Test\" & prot_per_seq >= 2, ]), \"(\", \n percent(nrow(dat[dataset == \"Test\" & prot_per_seq >= 2, ])/n_distinct(test_set$sequence), accuracy = 0.01),\")\")\n\nstat <- dstats[dstats$stat == \"Max\", \"Test (n = 139155)\"]\ncat(\"\\nSequences with\", stat,\n \"(max) protein names:\", nrow(dat[dataset == \"Test\" & prot_per_seq == stat, ]), \"(\", \n percent(nrow(dat[dataset == \"Test\" & prot_per_seq == stat, ])/n_distinct(test_set$sequence), accuracy = 0.01),\")\" )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-05-23T12:07:31.446127Z","iopub.execute_input":"2023-05-23T12:07:31.447759Z","iopub.status.idle":"2023-05-23T12:07:31.615217Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(\"Number of protein sequenses in datasets by different amount of protein names\")\nsummary <- dat %>% group_by(dataset, prot_per_seq) %>%\n summarise(n = n(), .groups = \"drop_last\") %>%\n mutate(freq = n / sum(n))\n\nsummary %>% mutate(freq = scales::percent(freq, accuracy = 0.01)) %>%\n pivot_wider(names_from = \"dataset\", values_from = c(\"n\", \"freq\"),\n values_fill = \"0\",\n values_fn = paste, names_vary = \"slowest\") %>%\n arrange(prot_per_seq)","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:07:38.261257Z","iopub.execute_input":"2023-05-23T12:07:38.262943Z","iopub.status.idle":"2023-05-23T12:07:38.425231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**There are a maximum of 17 proteins with the same sequence in both datasets**","metadata":{}},{"cell_type":"code","source":"test_set[prot_per_seq == 17, ] %>% mutate(sequence = paste0(substring(sequence, 1, 10), \"...\") )","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:07:52.589045Z","iopub.execute_input":"2023-05-23T12:07:52.590697Z","iopub.status.idle":"2023-05-23T12:07:52.636932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_set[prot_per_seq ==17, ] %>% mutate(sequence = paste0(substring(sequence, 1, 10), \"...\") )","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:07:59.389565Z","iopub.execute_input":"2023-05-23T12:07:59.391209Z","iopub.status.idle":"2023-05-23T12:07:59.442192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Checking the number of unique sequences with more than 1 protein name**","metadata":{}},{"cell_type":"code","source":"# Subset sequences with more than 1 protein name\nmulti_prot <- train_set_all[prot_per_seq > 1, ]\nmulti_prot <- multi_prot[, c('level', 'missing_parents', 'has_children', 'seq_length', 'prot_per_tax'):= NULL] %>% unique\ndata.summary(multi_prot)\n\n#multi_prot %>% mutate(sequence = paste0(substring(sequence, 1, 50), \"...\") ) %>% arrange(sequence)\nmulti_prot[, term_per_seq := uniqueN(term), by = sequence]","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:08:10.521041Z","iopub.execute_input":"2023-05-23T12:08:10.522725Z","iopub.status.idle":"2023-05-23T12:08:10.905103Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"-> There are 2731 sequences with more than 1 protein name (6053 proteins per 2731 sequences in total)","metadata":{}},{"cell_type":"markdown","source":"**Does annotaion of such proteins differ?**","metadata":{}},{"cell_type":"code","source":"# combine all GO terms of each protein in one row\n#sequence_terms <- train_set[, .(terms = paste(unique(term), collapse = \", \")), by = sequence]\nmulti_prot_stats <- multi_prot %>% group_by(sequence, EntryID, prot_per_seq) %>%\n summarise(annotation = list(term), .groups = \"drop\")\n\ncat(\"\\nTwo sequences with different protein names and their annotations: \")\nmulti_prot_stats %>% head(7) %>%\n mutate(sequence = paste0(substring(sequence, 1, 50), \"...\") )","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:08:16.210028Z","iopub.execute_input":"2023-05-23T12:08:16.213073Z","iopub.status.idle":"2023-05-23T12:08:16.302369Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"multi_ann_seq <- multi_prot_stats %>% \n group_by(sequence) %>%\n summarise(annotations_per_seq = n_distinct(annotation)) %>%\n filter(annotations_per_seq > 1) %>%\n #mutate(sequence = paste0(substring(sequence, 1, 50), \"...\") ) %>%\n nrow\n\ncat(\"For\", multi_ann_seq, \"of\", n_distinct(multi_prot_stats$sequence), \"(\",\n scales::percent(multi_ann_seq/n_distinct(multi_prot_stats$sequence)) ,\")\", \"sequences with more than 1 protein ID,\nthere are more than 1 annotation\")","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:08:28.70879Z","iopub.execute_input":"2023-05-23T12:08:28.710459Z","iopub.status.idle":"2023-05-23T12:08:28.880842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Example**","metadata":{}},{"cell_type":"code","source":"multi_prot_stats[multi_prot_stats$sequence == \"GILLDKLKNFAKTAGKGVLQSLLNTASCKLSGQC\", ]","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:08:35.686233Z","iopub.execute_input":"2023-05-23T12:08:35.688124Z","iopub.status.idle":"2023-05-23T12:08:35.712227Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Sequence length","metadata":{}},{"cell_type":"code","source":"cat(\"Number of amino acids in protein sequences varies:\")\ncat(\"\\nTrain set:\", min(train_set$seq_length), \"to\", max(train_set$seq_length))\ncat(\"\\nTest set:\", min(test_set$seq_length), \"to\", max(test_set$seq_length))","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:08:53.632409Z","iopub.execute_input":"2023-05-23T12:08:53.634416Z","iopub.status.idle":"2023-05-23T12:08:53.663839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dat <- rbind(train_set %>% select(\"sequence\", \"seq_length\") %>% unique %>% mutate(dataset = \"Train\"),\n test_set %>% select(\"sequence\", \"seq_length\") %>% unique %>% mutate(dataset = \"Test\"))\n\ncat(\"Distribution of the number of amino acids per sequence\")\ndstats <- merge(stat_calcs(dat[dat$dataset == \"Train\", ], colname = \"seq_length\") %>%\n dplyr::rename_with(.fn = ~paste0(\"Train (n = \",\n n_distinct(train_set$sequence),\n \")\"),\n .cols = seq_length),\n stat_calcs(dat[dat$dataset == \"Test\", ], colname = \"seq_length\") %>%\n dplyr::rename_with(.fn = ~paste0(\"Test (n = \",\n n_distinct(test_set$sequence),\n \")\"), .cols = seq_length) ) %>%\n arrange(stat)\ndstats","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:08:57.864652Z","iopub.execute_input":"2023-05-23T12:08:57.866271Z","iopub.status.idle":"2023-05-23T12:08:58.12905Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 10)\nggplot(dat, aes(x = seq_length, fill = dataset)) +\n geom_histogram(binwidth = 100, color = \"black\", alpha = 0.4, position = \"identity\") +\n scale_fill_manual(values = c(\"#cc7722\", \"#01796f\")) +\n labs(x = \"Sequence length (bin = 100)\", y = \"Count\",\n title = paste0(\"Distribution of the number of amino acids in protein sequences (n = \",\n n_distinct(test_set$EntryID), \" proteins)\")) +\n scale_x_continuous(breaks = c(seq(0, max(test_set$seq_length), 10000),\n 1000,\n max(test_set$seq_length))) +\n theme_bw(base_size = 20) +\n theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(),\n legend.position = c(0.95, 0.75))\n\nggplot(dat, aes(x = seq_length, fill = dataset)) +\n geom_histogram(binwidth = 0.05, color = \"black\", alpha = 0.4, position = \"identity\") +\n scale_x_log10() +\n scale_fill_manual(values = c(\"#cc7722\", \"#01796f\")) +\n labs(x = \"Sequence length (log10 scale)\", y = \"Count\") +\n theme_bw(base_size = 20) +\n theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(),\n legend.position = c(0.95, 0.75))","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:09:02.489449Z","iopub.execute_input":"2023-05-23T12:09:02.491132Z","iopub.status.idle":"2023-05-23T12:09:04.207804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#-------------train set--------------------#\ncat(\"\\033[1mTrain dataset\\033[0m\\n\")\n\nstat <- dstats[dstats$stat == \"Min\", \"Train (n = 138924)\"]\ncat(\"\\nSequences with only\", stat,\n \"(min) amino acids:\", nrow(dat[dataset == \"Train\" & seq_length <= stat, ]), \"(\", \n percent(nrow(dat[dataset == \"Train\" & seq_length <= stat, ])/n_distinct(train_set$sequence), accuracy = 0.01),\")\" )\n\ncat(\"\\nSequences with 100 or fewer amino acids:\",\n nrow(dat[dataset == \"Train\" & seq_length <= 100, ]), \"(\", \n percent(nrow(dat[dataset == \"Train\" & seq_length <= 100, ])/n_distinct(train_set$sequence), accuracy = 0.01),\")\")\n\nstat <- dstats[dstats$stat == \"Median\", \"Train (n = 138924)\"]\ncat(\"\\nSequences with\", stat,\n \"(median) or fewer amino acids:\", nrow(dat[dataset == \"Train\" & seq_length <= stat, ]))\n\nstat <- dstats[dstats$stat == \"90%\", \"Train (n = 138924)\"]\ncat(\"\\nSequences with\", stat,\n \"(90th percentile) or fewer amino acids:\", nrow(dat[dataset == \"Train\" & seq_length <= stat, ]))\ncat(\"\\n\\nSequences with more than\", stat,\n \"(90th percentile) amino acids:\", nrow(dat[dataset == \"Train\" & seq_length > stat, ]))\ncat(\"\\nSequences with 10000 of amino acids or more:\",\n nrow(dat[dataset == \"Train\" & seq_length >= 10000, ]), \"(\", \n percent(nrow(dat[dataset == \"Train\" & seq_length >= 10000, ])/n_distinct(train_set$sequence), accuracy = 0.01),\")\")\n\n#-------------test set--------------------#\ncat(\"\\n\\n\\033[1mTest dataset\\033[0m\\n\")\n\nstat <- dstats[dstats$stat == \"Min\", \"Test (n = 139155)\"]\ncat(\"\\nSequences with only\", stat,\n \"(min) amino acids:\", nrow(dat[dataset == \"Test\" & seq_length <= stat, ]), \"(\", \n percent(nrow(dat[dataset == \"Test\" & seq_length <= stat, ])/n_distinct(test_set$sequence), accuracy = 0.01),\")\" )\n\ncat(\"\\nSequences with 100 or fewer amino acids:\",\n nrow(dat[dataset == \"Test\" & seq_length <= 100, ]), \"(\", \n percent(nrow(dat[dataset == \"Test\" & seq_length <= 100, ])/n_distinct(test_set$sequence), accuracy = 0.01),\")\")\n\nstat <- dstats[dstats$stat == \"Median\", \"Test (n = 139155)\"]\ncat(\"\\nSequences with\", stat,\n \"(median) or fewer amino acids:\", nrow(dat[dataset == \"Test\" & seq_length <= stat, ]))\n\nstat <- dstats[dstats$stat == \"90%\", \"Test (n = 139155)\"]\ncat(\"\\nSequences with\", stat,\n \"(90th percentile) or fewer amino acids:\", nrow(dat[dataset == \"Test\" & seq_length <= stat, ]))\ncat(\"\\n\\nSequences with more than\", stat,\n \"(90th percentile) amino acids:\", nrow(dat[dataset == \"Test\" & seq_length > stat, ]))\ncat(\"\\nSequences with 10000 of amino acids or more:\",\n nrow(dat[dataset == \"Test\" & seq_length >= 10000, ]), \"(\", \n percent(nrow(dat[dataset == \"Test\" & seq_length >= 10000, ])/n_distinct(test_set$EntryID), accuracy = 0.01),\")\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-05-23T12:09:42.319931Z","iopub.execute_input":"2023-05-23T12:09:42.322141Z","iopub.status.idle":"2023-05-23T12:09:42.622981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Taxons","metadata":{}},{"cell_type":"code","source":"cat(\"Number fo taxons in the train set:\", n_distinct(train_set$taxonomyID))\ncat(\"\\nNumber fo taxons in the test set:\", n_distinct(test_set$taxonomyID))","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:09:42.627197Z","iopub.execute_input":"2023-05-23T12:09:42.629296Z","iopub.status.idle":"2023-05-23T12:09:42.663593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 10)\nvenn_diag <- list(\n \"Train set\" = unique(train_set$taxonomyID), \n \"Test set\" = unique(test_set$taxonomyID))\n\nggvenn(venn_diag, show_elements = FALSE, label_sep = \"\\n\",\n text_size = 6, set_name_size = 8,\n fill_color = c(\"#cc7722\", \"#01796f\")) +\nggtitle(paste(\"Intersection of taxons in the datasets\")) +\ntheme(plot.title = element_text(size = 20, face = \"bold\", hjust = 0.5))","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:09:42.667556Z","iopub.execute_input":"2023-05-23T12:09:42.669595Z","iopub.status.idle":"2023-05-23T12:09:43.528937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dat <- rbind(train_set %>% select(\"taxonomyID\", \"prot_per_tax\") %>% unique %>% mutate(dataset = \"Train\"),\n test_set %>% select(\"taxonomyID\", \"prot_per_tax\") %>% unique %>% mutate(dataset = \"Test\"))\n\ncat(\"Distribution of the number of proteins per taxon\")\ndstats <- merge(stat_calcs(dat[dat$dataset == \"Train\", ], colname = \"prot_per_tax\") %>%\n dplyr::rename_with(.fn = ~paste0(\"Train (n = \",\n n_distinct(train_set$taxonomyID),\n \")\"),\n .cols = prot_per_tax),\n stat_calcs(dat[dat$dataset == \"Test\", ], colname = \"prot_per_tax\") %>%\n dplyr::rename_with(.fn = ~paste0(\"Test (n = \",\n n_distinct(test_set$taxonomyID),\n \")\"), .cols = prot_per_tax) ) %>%\n arrange(stat)\ndstats","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:09:43.532009Z","iopub.execute_input":"2023-05-23T12:09:43.533627Z","iopub.status.idle":"2023-05-23T12:09:43.729316Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#-------------train set--------------------#\ncat(\"\\033[1mTrain dataset (n =\", n_distinct(train_set$taxonomyID), \"taxons)\\033[0m\\n\")\n\nstat <- dstats[dstats$stat == \"Median\", \"Train (n = 3156)\"]\ncat(\"\\nTaxons with\", stat,\n \"(median) or fewer observations:\", nrow(dat[dataset == \"Train\" & prot_per_tax <= stat, ]))\n\nstat <- dstats[dstats$stat == \"90%\", \"Train (n = 3156)\"]\ncat(\"\\nTaxons with\", stat,\n \"(90th percentile) or fewer observations:\", nrow(dat[dataset == \"Train\" & prot_per_tax <= stat, ]))\ncat(\"\\n\\nTaxons with more than\", stat,\n \"(90th percentile) observations:\", nrow(dat[dataset == \"Train\" & prot_per_tax > stat, ]))\ncat(\"\\nTaxons with 1000 of observations or more:\",\n nrow(dat[dataset == \"Train\" & prot_per_tax >= 1000, ]), \"(\", \n percent(nrow(dat[dataset == \"Train\" & prot_per_tax >= 1000, ])/n_distinct(train_set$taxonomyID), accuracy = 0.01),\")\")\n\n#-------------test set--------------------#\ncat(\"\\n\\n\\033[1mTest dataset (n =\", n_distinct(test_set$taxonomyID), \"taxons)\\033[0m\\n\")\n\nstat <- dstats[dstats$stat == \"Min\", \"Test (n = 90)\"]\ncat(\"\\nTaxons with only\", stat,\n \"(min) observations:\", nrow(dat[dataset == \"Test\" & prot_per_tax <= stat, ]), \"(\", \n percent(nrow(dat[dataset == \"Test\" & prot_per_tax <= stat, ])/n_distinct(test_set$taxonomyID), accuracy = 0.01),\")\" )\n\nstat <- dstats[dstats$stat == \"Median\", \"Test (n = 90)\"]\ncat(\"\\nTaxons with\", stat,\n \"(median) or fewer observations:\", nrow(dat[dataset == \"Test\" & prot_per_tax <= stat, ]))\n\nstat <- dstats[dstats$stat == \"90%\", \"Test (n = 90)\"]\ncat(\"\\nTaxons with\", stat,\n \"(90th percentile) or fewer observations:\", nrow(dat[dataset == \"Test\" & prot_per_tax <= stat, ]))\ncat(\"\\n\\nTaxons with more than\", stat,\n \"(90th percentile) observations:\", nrow(dat[dataset == \"Test\" & prot_per_tax > stat, ]))\ncat(\"\\nTaxons with 1000 of observations or more:\",\n nrow(dat[dataset == \"Test\" & prot_per_tax >= 1000, ]), \"(\", \n percent(nrow(dat[dataset == \"Test\" & prot_per_tax >= 1000, ])/n_distinct(test_set$taxonomyID), accuracy = 0.01),\")\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-05-23T12:09:43.732342Z","iopub.execute_input":"2023-05-23T12:09:43.733922Z","iopub.status.idle":"2023-05-23T12:09:43.844482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(\"Taxons in the test set with only 1 (min) observarion:\")\ntest_set[prot_per_tax == 1, ] %>%\n select(-c('sequence', 'seq_length'))\n\ncat(\"\\nTaxons in the test with more than 4214 (90th percentile) observarions:\")\nsetorder(test_set[prot_per_tax > 4214], - prot_per_tax) %>%\n select(-c('EntryID', 'sequence', 'seq_length', 'prot_per_seq')) %>% unique","metadata":{"execution":{"iopub.status.busy":"2023-05-23T12:09:43.848434Z","iopub.execute_input":"2023-05-23T12:09:43.850533Z","iopub.status.idle":"2023-05-23T12:09:43.947599Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#paste(test_set[70001:nrow(test_set), \"EntryID\"], collapse = \", \")\n#paste(test_set[1:70000, \"EntryID\"], collapse = \", \")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-05-23T12:09:43.951493Z","iopub.execute_input":"2023-05-23T12:09:43.953584Z","iopub.status.idle":"2023-05-23T12:09:43.966557Z"},"trusted":true},"execution_count":null,"outputs":[]}]}