R/auc_thresh_kmeans_nextflow.R

In FloWuenne/scFunctions: Functions for single cell data analysis

#' This function is to analyze the nextflow implementation results of SCENIC!
#' This function takes the regulon AUC values and determines thresholds for binarizatio based on kmeans clustering.
#' Less specific than the SCENIC function AUCell_exploreThresholds but much faster.
#'
#'
#' @param regulonAUC The AUC values for all regulons as calculated by SCENIC (content of file:3.4_regulonAUC.Rds).
#' @keywords SCENIC, regulons, binary activity, kmeans, thresholds
#' @import SCENIC
#' @import svMisc
#' @import tidyverse
#' @export
#' @examples
#' \donttest{
#' regulon_thresholds <- auc_thresh_kmeans(regulonAUC)
#' #' }


auc_thresh_kmeans_nextflow <- function(regulonAUC){
    
  ## Iterate over each regulon in the AUC matrix
  regulons <- colnames(regulonAUC)
  
  kmeans_thresholds <- list()
  
  print("Processing regulon distributions...")
  
  for(regulon_no in 1:length(regulons)){
    
    progress(regulon_no)
    
    regulon <- regulons[regulon_no]
    
    df <- data.frame("auc" = regulonAUC[,regulon_no],
                     "cells"= rownames(regulonAUC),
                     "regulon" = regulon)
    
    ## Remove cells with 0 AUC as they interfere with kmeans clustering
    df <- df %>%
      subset(auc > 0)
    
    km <- kmeans(df$auc,centers=2)
    df$cluster <- as.factor(km$cluster)
    
    cluster1_max <- max(subset(df,cluster == 1)$auc)
    cluster2_max <- max(subset(df,cluster == 2)$auc)
    
    if(cluster1_max > cluster2_max){
      
      df <- df %>%
        mutate("cluster" = gsub(2,3,cluster)) %>%
        mutate("cluster" = gsub(1,2,cluster)) %>%
        mutate("cluster" = gsub(3,1,cluster))
    }
    
    df <- df %>%
      arrange(desc(auc))
    
    df_sub <- df %>%
      subset(cluster == 1)
    
    auc_thresholds <- df_sub[1,]$auc
    
    kmeans_thresholds[[regulon]] <- auc_thresholds
    
  }
  
  
  print("Done evaluating thresholds...")
  return(kmeans_thresholds)
  
}