geposanui/clustering.R

library(data.table)

#' Perform a cluster analysis.
#'
#' This function will cluster the data using `hclust` and `cutree` (with the
#' specified height). Every cluster with at least two members qualifies for
#' further analysis. Clusters are then ranked based on their size in relation
#' to the total number of possible values (`n`). The return value is a final
#' score between zero and one. Lower ranking clusters contribute less to this
#' score.
clusteriness <- function(data, n, height = 1000000) {
    # Return a score of 0.0 if there is just one or no value at all.
    if (length(data) < 2) {
        return(0.0)
    }

    # Cluster the data and compute the cluster sizes.

    tree <- hclust(dist(data))
    clusters <- cutree(tree, h = height)
    cluster_sizes <- sort(tabulate(clusters), decreasing = TRUE)

    # Compute the "clusteriness" score.

    score <- 0.0

    for (i in seq_along(cluster_sizes)) {
        cluster_size <- cluster_sizes[i]

        if (cluster_size >= 2) {
            cluster_score <- cluster_size / n
            score <- score + cluster_score / i
        }
    }

    score
}

#' Process genes clustering their distance to telomeres.
#'
#' The return value will be a data.table with the following columns:
#'
#'  - `gene` Gene ID of the processed gene.
#'  - `clusteriness` Score quantidying the gene's clusters.
#'
#' @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_clustering <- function(distances, species_ids, gene_ids) {
    results <- data.table(gene = gene_ids)
    species_count <- length(species)

    # Prefilter the input data by species.
    distances <- distances[species %chin% species_ids]

    # Add an index for quickly accessing data per gene.
    setkey(distances, gene)

    #' Perform the cluster analysis for one gene.
    compute <- function(gene_id) {
        clusteriness(distances[gene_id, distance], species_count)
    }

    results[, clusteriness := compute(gene), by = 1:nrow(results)]
}