johrpan/geposanui
1
0
Fork 0
mirror of https://github.com/johrpan/geposanui.git synced 2025-10-26 19:27:24 +01:00
Code Issues Projects Releases Packages Wiki Activity Actions

Update correlation method

Browse source 
This adds some performance optimizations and removes the penalty for
missing values.
This commit is contained in:
Elias Projahn 2021-10-11 10:10:28 +02:00
parent cd63486f0a
commit d87474e574
1 changed files with 29 additions and 27 deletions
Download patch file Download diff file Expand all files Collapse all files

56
correlation.R
View file

@ -1,6 +1,4 @@
library(data.table) library(data.table)
library(progress)
library(rlog)
#' Compute the mean correlation coefficient comparing gene distances with a set #' Compute the mean correlation coefficient comparing gene distances with a set
#' of reference genes. #' of reference genes.
@ -8,7 +6,7 @@ library(rlog)
#' The result will be a data.table with the following columns: #' The result will be a data.table with the following columns:
#' #'
#' - `gene` Gene ID of the processed gene. #' - `gene` Gene ID of the processed gene.
#' - `r_mean` Mean correlation coefficient. #' - `correlation` Mean correlation coefficient.
#' #'
#' @param distances Distance data to use. #' @param distances Distance data to use.
#' @param species_ids Species, whose data should be included. #' @param species_ids Species, whose data should be included.
@ -17,55 +15,59 @@ library(rlog)
process_correlation <- function(distances, species_ids, gene_ids, process_correlation <- function(distances, species_ids, gene_ids,
reference_gene_ids) { reference_gene_ids) {
results <- data.table(gene = gene_ids) results <- data.table(gene = gene_ids)
gene_count <- length(gene_ids)
reference_count <- length(reference_gene_ids) reference_count <- length(reference_gene_ids)
log_info(sprintf(
"Correlating %i genes from %i species with %i reference genes",
gene_count,
length(species_ids),
reference_count
))
progress <- progress_bar$new(
total = gene_count,
format = "Correlating genes [:bar] :percent (ETA :eta)"
)
# Prefilter distances by species. # Prefilter distances by species.
distances <- distances[species %chin% species_ids] distances <- distances[species %chin% species_ids]
for (i in 1:gene_count) { # Add an index for quickly accessing data per gene.
progress$tick() setkey(distances, gene)
gene_id <- gene_ids[i] # Prepare the reference genes' data.
gene_distances <- distances[gene == gene_id] reference_distances <- distances[gene %chin% reference_gene_ids]
if (nrow(gene_distances) < 10) { #' Perform the correlation for one gene.
next compute <- function(gene_id) {
gene_distances <- distances[gene_id]
gene_species_count <- nrow(gene_distances)
# Return a score of 0.0 if there is just one or no value at all.
if (gene_species_count <= 1) {
return(0.0)
} }
#' Buffer for the sum of correlation coefficients. #' Buffer for the sum of correlation coefficients.
r_sum <- 0 correlation_sum <- 0
# Correlate with all reference genes but not with the gene itself. # Correlate with all reference genes but not with the gene itself.
for (reference_gene_id in for (reference_gene_id in
reference_gene_ids[reference_gene_ids != gene_id]) { reference_gene_ids[reference_gene_ids != gene_id]) {
data <- merge( data <- merge(
gene_distances, gene_distances,
distances[gene == reference_gene_id], reference_distances[reference_gene_id],
by = "species" by = "species"
) )
# Skip this reference gene, if there are not enough value pairs.
# This will lessen the final score, because it effectively
# represents a correlation coefficient of 0.0.
if (nrow(data) <= 1) {
next
}
# Order data by the reference gene's distance to get a monotonic # Order data by the reference gene's distance to get a monotonic
# relation. # relation.
setorder(data, distance.y) setorder(data, distance.y)
r_sum <- r_sum + abs(cor(data[, distance.x], data[, distance.y])) correlation_sum <- correlation_sum + abs(cor(
data[, distance.x], data[, distance.y],
method = "spearman"
))
} }
results[gene == gene_id, r_mean := r_sum / reference_count] # Compute the score as the mean correlation coefficient.
score <- correlation_sum / reference_count
} }
results results[, correlation := compute(gene), by = 1:nrow(results)]
} }