From 9d6b2e4d50ac891217f32fd5db721e3c0062b04f Mon Sep 17 00:00:00 2001 From: Elias Projahn Date: Sat, 18 Sep 2021 23:10:52 +0200 Subject: [PATCH] Add new correlation method --- process.R => clustering.R | 21 ++++++--- correlation.R | 61 ++++++++++++++++++++++++ init.R | 97 +++++++++++++++++++++++++++++++++++++++ server.R | 60 ++---------------------- 4 files changed, 176 insertions(+), 63 deletions(-) rename process.R => clustering.R (65%) create mode 100644 correlation.R create mode 100644 init.R diff --git a/process.R b/clustering.R similarity index 65% rename from process.R rename to clustering.R index b3b26ff..e0b5a3e 100644 --- a/process.R +++ b/clustering.R @@ -1,21 +1,30 @@ library(data.table) library(rlog) -#' Process genes screening for a likely TPE-OLD. +#' Process genes clustering their distance to telomeres. #' -#' The return value will be a table containing genes and data to take in -#' account when regarding them as TPE-OLD candidates. +#' The return value will be a data.table with the following columns: +#' +#' - `gene` Gene ID of the processed gene. +#' - `cluster_length` Length of the largest cluster. +#' - `cluster_mean` Mean value of the largest cluster. +#' - `cluster_species` List of species contributing to the largest cluster. #' #' @param distances Gene distance data to use. #' @param species_ids IDs of species to include in the analysis. #' @param gene_ids Genes to include in the computation. -process_input <- function(distances, species_ids, gene_ids) { +process_clustering <- function(distances, species_ids, gene_ids) { results <- data.table(gene = gene_ids) gene_count <- length(gene_ids) - for (i in seq_along(gene_ids)) { + for (i in 1:gene_count) { gene_id <- gene_ids[i] - log_info(sprintf("Processing gene %i/%i (%s)", i, gene_count, gene_id)) + + log_info(sprintf( + "[%3i%%] Processing gene \"%s\"", + round(i / gene_count * 100), + gene_id + )) data <- distances[ species %chin% species_ids & gene == gene_id, diff --git a/correlation.R b/correlation.R new file mode 100644 index 0000000..f3bb536 --- /dev/null +++ b/correlation.R @@ -0,0 +1,61 @@ +library(data.table) +library(rlog) + +#' Compute the mean correlation coefficient comparing gene distances with a set +#' of reference genes. +#' +#' The result will be a data.table with the following columns: +#' +#' - `gene` Gene ID of the processed gene. +#' - `r_mean` Mean correlation coefficient. +#' +#' @param distances Distance data to use. +#' @param species_ids Species, whose data should be included. +#' @param gene_ids Genes to process. +#' @param reference_gene_ids Genes to compare to. +process_correlation <- function(distances, species_ids, gene_ids, + reference_gene_ids) { + log_info("Processing genes for correlation") + + results <- data.table(gene = gene_ids) + gene_count <- length(gene_ids) + reference_count <- length(reference_gene_ids) + + # Prefilter distances by species. + distances <- distances[species %chin% species_ids] + + for (i in 1:gene_count) { + gene_id <- gene_ids[i] + + log_info(sprintf( + "[%3i%%] Processing gene \"%s\"", + round(i / gene_count * 100), + gene_id + )) + + gene_distances <- distances[gene == gene_id] + + if (nrow(gene_distances) < 12) { + next + } + + #' Buffer for the sum of correlation coefficients. + r_sum <- 0 + + # Correlate with all reference genes but not with the gene itself. + for (reference_gene_id in + reference_gene_ids[reference_gene_ids != gene_id]) { + data <- merge( + gene_distances, + distances[gene == reference_gene_id], + by = "species" + ) + + r_sum <- r_sum + abs(cor(data[, distance.x], data[, distance.y])) + } + + results[gene == gene_id, r_mean := r_sum / reference_count] + } + + results +} \ No newline at end of file diff --git a/init.R b/init.R new file mode 100644 index 0000000..1c29603 --- /dev/null +++ b/init.R @@ -0,0 +1,97 @@ +source("clustering.R") +source("correlation.R") +source("input.R") +source("util.R") + +# Load input data + +species <- run_cached("input/species", retrieve_species) +genes <- run_cached("input/genes", retrieve_genes) + +distances <- run_cached( + "input/distances", + retrieve_distances, + species[, id], + genes[, id] +) + +# Load processed data + +all_species <- species[, id] +replicative_species <- species[replicative == TRUE, id] +all_genes <- genes[, id] +tpe_old_genes <- genes[suggested | verified == TRUE, id] + +clustering_all <- run_cached( + "all_species/clustering", + process_clustering, + distances, + all_species, + all_genes +) + +clustering_replicative <- run_cached( + "replicative_species/clustering", + process_clustering, + distances, + replicative_species, + all_genes +) + +correlation_all <- run_cached( + "all_species/correlation", + process_correlation, + distances, + all_species, + all_genes, + tpe_old_genes +) + +correlation_replicative <- run_cached( + "replicative_species/correlation", + process_correlation, + distances, + replicative_species, + all_genes, + tpe_old_genes +) + +# Merge processed data as well as gene information. + +results_all <- merge( + genes, + clustering_all, + by.x = "id", + by.y = "gene" +) + +results_all <- merge( + results_all, + correlation_all, + by.x = "id", + by.y = "gene" +) + +results_replicative <- merge( + genes, + clustering_replicative, + by.x = "id", + by.y = "gene" +) + +results_replicative <- merge( + results_replicative, + correlation_replicative, + by.x = "id", + by.y = "gene" +) + +# Rename `id` columns to `gene`. + +setnames(results_all, "id", "gene") +setnames(results_replicative, "id", "gene") + +# Order results by cluster length descendingly as a start. + +setorder(results_all, -cluster_length) +setorder(results_replicative, -cluster_length) \ No newline at end of file diff --git a/server.R b/server.R index f4de2c8..cc7cdd7 100644 --- a/server.R +++ b/server.R @@ -2,61 +2,8 @@ library(data.table) library(DT) library(shiny) -source("input.R") -source("process.R") +source("init.R") source("scatter_plot.R") -source("util.R") - -# Load input data - -species <- run_cached("species", retrieve_species) -genes <- run_cached("genes", retrieve_genes) - -distances <- run_cached( - "distances", - retrieve_distances, - species[, id], - genes[, id] -) - -#' Results computed for all species. -results_all <- run_cached( - "results_all", - process_input, - distances, - species[, id], - genes[, id] -) - -#' Results computed for known replicatively aging species. -results_replicative <- run_cached( - "results_replicative", - process_input, - distances, - species[replicative == TRUE, id], - genes[, id] -) - -# Add gene information to results for display. - -results_all <- merge( - results_all, - genes, - by.x = "gene", - by.y = "id" -) - -results_replicative <- merge( - results_replicative, - genes, - by.x = "gene", - by.y = "id" -) - -# Order results by cluster length descendingly. -# TODO: Once other methods have been added, this has to be dynamic. -setorder(results_all, -cluster_length) -setorder(results_replicative, -cluster_length) server <- function(input, output) { #' This expression applies all user defined filters to the available @@ -77,14 +24,13 @@ server <- function(input, output) { output$genes <- renderDT({ datatable( - results()[, .(.I, name, chromosome, cluster_length, cluster_mean)], + results()[, .(.I, name, cluster_length, r_mean)], rownames = FALSE, colnames = c( "Rank", "Gene", - "Chromosome", "Cluster length", - "Cluster mean" + "Correlation" ), style = "bootstrap" )