R/subFunctions.R

In MultiSEp: Predict Synthetic Lethality and Other Gene Dependency Relationships from Multiomics Data

kmeansplusplus <- function(x,k,iter.max = 10, algorithm=c("Hartigan-Wong","Lloyd","Forgy","MacQueen")) {
  if(k<1) stop("k must be >=1")
  if(is.null(dim(x))) dim(x) <- c(length(x),1)
  if(k>nrow(x)) stop("k must be >n")
  centres <- matrix(0,k,ncol(x))
  centres[1,] <- x[sample(1:nrow(x),1),]
  if(k>1) {
    for(cn in 2:k) {
      dist.nearest.centroid <- apply(x,1,function(xeach){min(apply(t(t(centres)-xeach)^2,1,sum))})
      # Avoid divide by zero
      dist.nearest.centroid[dist.nearest.centroid < .Machine$double.eps] <- .Machine$double.eps
      centres[cn,] <- x[sample(1:nrow(x),1,prob=dist.nearest.centroid/sum(dist.nearest.centroid)),]
    }
  } else {
    centres <- 1
  }
  if(any(duplicated(centres))) {
    warning('Unable to generate ',k,' distinct centres, falling back to ordinary kmeans...')
    centres <- k
  }
  return(kmeans(x,centres,iter.max=iter.max,algorithm=algorithm))
}

discretise <- function(dataset,method="quantile",quantiles=4,model=NULL,...) {
  if(method!="mog" && !is.null(model)) stop("Cannot use a model unless method is 'mog'")
  if(!is.null(model) && !inherits(model, "mixmodel")) stop("Model is not a valid 'mixmodel'")
  drop.dims <- FALSE
  if(is.null(dim(dataset))) {
    nm <- names(dataset)
    drop.dims <- TRUE
    dim(dataset) <- c(1,length(dataset))
    colnames(dataset) <- nm
  }
  dataset <- na.omit(dataset)
  res <- matrix(NA,nrow=nrow(dataset),ncol=ncol(dataset),dimnames=dimnames(dataset))
  for(i in 1:nrow(dataset)) {
    data.vec <- as.numeric(dataset[i,])
    switch(method,
           quantile={
             res[i,] <- quantcut(data.vec,q=seq(0,1,1/quantiles))
             levels(res[i,]) <- 1:quantiles
           },
           mog={
             if(is.null(model)) model <- EM.findk(data.vec,...)
             if(model$k==1) {
               res[i,] <- rep(1,ncol(dataset)) ##unimodal (set all to 1)
             } else {
               res[i,] <- mog.density(data.vec,model)$membership
             }
           },
           stop("recognised values for method are 'quantile' or 'mog' (mixture of Gaussians)")
    )
  }

  # Drop unimodal markers
  if(method=="mog") {
    unimodal <- apply(res==1,1,all)
    if(any(unimodal)) message("Dropping unimodal markers:",paste(names(unimodal)[unimodal],collapse=","))
    res <- res[!unimodal,]
    attr(res,"unimodal") <- ifelse(any(unimodal),which(unimodal),NA)
  }
  attr(res,"na.action") <- attr(dataset,"na.action")

  if(drop.dims && !is.null(nrow(res)) && nrow(res)==1) {
    res <- c(res)
    names(res) <- colnames(dataset)
  }
  if(!is.null(attr(dataset,"protein")))
    attr(res,"protein") <- attr(dataset,"protein")
  attr(res,"discretisation") <- ifelse(method=="mog","Mixture of Gaussians","Quantiles")

  return(res)
}

log_likelihood <- function(data.vec,mixmodel) {
  data.vec <- as.numeric(data.vec)
  probs <- repmat(t(mixmodel$mix.props),length(data.vec),1) * sapply(1:mixmodel$k,function(i) { dnorm(data.vec,mixmodel$means[i],mixmodel$std.dev[i]) })
  sum(log(rowSums(probs)))
}


init.em <- function(data.vec,k,model.type="V",init.method="kmeans++") {
  if(init.method=="quantiles")
    q <- discretise(data.vec,method="quantile",quantiles=k)
  if(init.method=="kmeans++")
    q <- kmeansplusplus(data.vec,k)$cluster
  mix.props <- as.numeric(repmat(1/k,1,k))
  means <- sapply(unique(q),function(i) mean(data.vec[i==q]))
  std.dev <- if(model.type=="V") sapply(unique(q),function(i) sd(data.vec[i==q]))
  else if(model.type=="E") rep(sd(data.vec),k)
  else stop("Unknown model type: ",model.type)
  mixmodel <- list(initialisation=init.method,k=k,means=means,std.dev=std.dev,mix.props=mix.props,model.type=model.type,n=length(data.vec))
  class(mixmodel) <- "mixmodel"
  return(mixmodel)
}

e.step <- function(data.vec,mixmodel) {
  if(!inherits(mixmodel, "mixmodel")) stop("Mixture parameters should be of mixmodel class")
  data.vec <- as.numeric(data.vec)
  probs <- repmat(t(mixmodel$mix.props),length(data.vec),1) * sapply(1:mixmodel$k,function(i) { dnorm(data.vec,mixmodel$means[i],mixmodel$std.dev[i]) })
  probs/rowSums(probs)
}

m.step <- function(data.vec,mem.probs,mixmodel) {
  data.vec <- as.numeric(data.vec)
  mixmodel$mix.props <- colSums(mem.probs)/length(data.vec)
  mixmodel$means <- colSums(data.vec*mem.probs)/colSums(mem.probs)
  data.vec.mean.centred <- data.vec - repmat(t(mixmodel$means),length(data.vec),1)
  mixmodel$std.dev <- if(mixmodel$model.type=="V") sqrt(colSums(mem.probs*(data.vec.mean.centred^2))/colSums(mem.probs)) else if(mixmodel$model.type=="E") rep(sqrt(sum(mem.probs*(data.vec.mean.centred^2))/sum(mem.probs)),mixmodel$k) else stop("Unknown model type: ",mixmodel$model.type)
  mixmodel
}

EM <- function(data.vec,k,model.type="V",it.max=1e3,lhood.tol=1e-4,init.method="kmeans++",max.restarts=20,always.restart=FALSE,min.mean.diff=0.01,min.sd=0.01,suppress.warnings=FALSE) {
  # Remove missing values
  if(any(is.na(data.vec))) {
    marker.name <- attr(data.vec,"protein")
    data.vec <- na.omit(data.vec)
    attr(data.vec,"protein") <- marker.name
  }

  if(length(unique(data.vec))<k) stop("Need at least k=",k,"distinct data points")

  mixmodel <- NULL
  n.restart <- 0
  n.unconverged <- 0
  while(n.restart==0 || ((always.restart || any(diff(sort(mixmodel.this$means))<min.mean.diff,na.rm=T) || any(mixmodel.this$std.dev<min.sd,na.rm=T)) && n.restart<max.restarts)) {
    # Initialise
    mixmodel.this <- init.em(data.vec,k,model.type=model.type,init.method=init.method)

    lhood.last <- log_likelihood(data.vec,mixmodel.this)

    i <- 1
    repeat {
      ##Read off the new membership probabilities
      mem.probs <- e.step(data.vec,mixmodel.this)
      ##Maximisation step - adapt parameters
      mixmodel.this <- m.step(data.vec,mem.probs,mixmodel.this)
      ##Compute the likelihood
      lhood.this <- log_likelihood(data.vec,mixmodel.this)
      ##Stop if unavailable to compute likelihood
      if(is.na(lhood.this)) break
      ##Stop on convergence
      if(abs((lhood.last-lhood.this)/lhood.last)<lhood.tol) break
      ## Stop if maximum number of iterations reached
      if(i>=it.max) {
        n.unconverged <- n.unconverged+1
        break
      }
      lhood.last <- lhood.this
      i <- i+1
    }
    if(is.null(mixmodel) || is.null(mixmodel$loglik) || is.na(mixmodel$loglik) || mixmodel$loglik>mixmodel.this$loglik) mixmodel <- mixmodel.this
    n.restart <- n.restart+1
  }

  # Sort the MoG modes in ascending order
  # Previous code ordered by means:  mode.order <- order(mixmodel$means)
  mode.order <- order(mixmodel$means - (mixmodel$std.dev)) # return indices of mixmodel according to value of one standard deviation below the mean
  mixmodel$means <- mixmodel$means[mode.order]
  mixmodel$std.dev <- mixmodel$std.dev[mode.order]
  mixmodel$mix.props <- mixmodel$mix.props[mode.order]

  if(!suppress.warnings) {
    if(any(is.na(mixmodel$mix.props))) warning("k=",k,": ",sum(is.na(mixmodel$mix.props))," empty model component(s) (decrease k or increase max.restarts)")
    if(n.unconverged==max.restarts) warning("k=",k,": Model did not converge (increase it.max, lhood.tol or max.restarts)")
    else if(any(diff(mixmodel$means)<min.mean.diff,na.rm=T)) warning("k=",k,": Some component means closer than requested (decrease k or increase max.restarts)")
    else if(any(mixmodel$std.dev<min.sd,na.rm=T)) warning("k=",k,": Some component standard deviations lower than requested (decrease k or increase max.restarts)")
  }

  mixmodel$loglik <- lhood.this
  mixmodel$bic <- tmap.bic(lhood.this,mixmodel,length(data.vec))
  mixmodel$aic <- tmap.aic(lhood.this,mixmodel)

  attr(mixmodel,"protein") <- attr(data.vec,"protein")
  class(mixmodel) <- "mixmodel"
  return(mixmodel)
}

tmap.bic <- function(loglik,mixmodel,n) {
  k <- ifelse(mixmodel$model.type=="V",2*mixmodel$k,mixmodel$k+1)
  # The mixture components count as k-1 free parameters
  if(mixmodel$k>1) k <- k+mixmodel$k-1
  -(2*loglik)+(k*log(n))
}

tmap.aic <- function(loglik,mixmodel) {
  k <- ifelse(mixmodel$model.type=="V",2*mixmodel$k,mixmodel$k+1)
  # The mixture components count as k-1 free parameters
  if(mixmodel$k>1) k <- k+mixmodel$k-1
  -(2*loglik)+(2*k)
}

mog.density <- function(data.vec,mixmodel,sort.data=FALSE) {
  if(any(is.na(data.vec))) stop("Please remove NA values from data.vec (use na.omit)")
  if(sort.data) data.vec <- sort(data.vec)
  y <- matrix(NA,mixmodel$k,length(data.vec))
  for(i in 1:mixmodel$k) y[i,] <- mixmodel$mix.props[i]*dnorm(as.numeric(data.vec),mean=mixmodel$means[i],sd=mixmodel$std.dev[i])
  prob <- apply(y,2,function(a){a/sum(a)})
  if(is.null(dim(prob))) dim(prob) <- c(1,length(prob))
  list(x=data.vec,y=apply(y,2,sum),prob=t(prob),membership=apply(y,2,which.max))
}

EM.findk <- function(data.vec,num.gaussians=1:9,model.types=c("E","V"),select.by="bic",suppress.warnings=TRUE,...) {
  stopifnot(select.by %in% c("aic","bic"))
  model.bic <- list()
  mixmodels <- list()

  best.k <- NA
  best.type <- NA
  best.bic <- Inf

  for(k in num.gaussians) {
    for(model.type in model.types) {
      mixid <- paste(k,model.type,sep="")
      tryCatch({
        mixmodels[[mixid]] <- EM(data.vec,k=k,model.type=model.type,suppress.warnings=suppress.warnings,...)
        model.bic[[mixid]] <- mixmodels[[mixid]][[select.by]]
        if(!is.na(model.bic[[mixid]]) && model.bic[[mixid]]<best.bic) {
          best.bic <- model.bic[[mixid]]
          best.k <- k
          best.type <- model.type
        }
      },error=function(e){if(!suppress.warnings) {
        warning("Mixture model k=",k,"type=",model.type,"failed to generate")
      }})
    }
  }

  if(is.na(best.k)) stop("Unable to generate a mixture model with any value of k in the supplied range")

  mixmodels[[paste(best.k,best.type,sep="")]]
}

prefilterMixModelClusters <- function(mixModelClusters) {
  mixModelClusters <- mixModelClusters[sapply(mixModelClusters, function(x) !all(is.atomic(x)))]
  return(mixModelClusters)
}

generateGenePairs <- function(mixModelClusters1,
                              mixModelClusters2 = NULL,
                              bidirectionalAnalysis = FALSE,
                              include_reverse_pairs = FALSE,
                              verbose = TRUE,
                              cores) {
  
  if (verbose) message("Generating genepairs for analysis")
  expression_genes <- names(mixModelClusters1)
  
  if (bidirectionalAnalysis) {
    # ___ 2-GMM analysis ___
    crispr_genes <- names(mixModelClusters2)
    
    # unique pairs with pbmclapply
    if (length(expression_genes) <= length(crispr_genes)) {
      genepairs_list <- pbmclapply(expression_genes, function(g1) {
        lapply(crispr_genes, function(g2) {
          if (g1 != g2) sort(c(g1, g2)) else NULL
        })
      }, mc.cores = cores)
    } else {
      genepairs_list <- pbmclapply(crispr_genes, function(g2) {
        lapply(expression_genes, function(g1) {
          if (g1 != g2) sort(c(g1, g2)) else NULL
        })
      }, mc.cores = cores)
    }
    
    # combine and deduplicate
    genepairs <- unique(as.data.frame(
      do.call(rbind, unlist(genepairs_list, recursive = FALSE)),
      stringsAsFactors = FALSE
    ))
    colnames(genepairs) <- c("Gene1", "Gene2")
    
    # if reverse pairs is true, add only valid reverse directions
    if (include_reverse_pairs) {
      genepairs_rev <- genepairs[
        genepairs$Gene1 %in% names(mixModelClusters2) &
          genepairs$Gene2 %in% names(mixModelClusters1),
      ]
      if (nrow(genepairs_rev) > 0) {
        genepairs_rev <- genepairs_rev[, c("Gene2", "Gene1")]
        colnames(genepairs_rev) <- c("Gene1", "Gene2")
        genepairs <- unique(rbind(genepairs, genepairs_rev))
      }
    }
    
  } else {
    # --- 1-GMM analysis ---
    if (length(expression_genes) < 2) {
      stop("Analysis requires at least two genes in mixModelClusters1 to generate gene pairs.")
    }
    n <- length(expression_genes)
    genepairs <- pbmclapply(1:(n - 1), function(i) {
      pairs <- lapply((i + 1):n, function(j) {
        c(expression_genes[i], expression_genes[j])
      })
      do.call(rbind, pairs)
    }, mc.cores = cores)
    
    genepairs <- as.data.frame(do.call(rbind, genepairs))
    colnames(genepairs) <- c("Gene1", "Gene2")
  }
  
  if (verbose) message("Generated ", nrow(genepairs), " gene pairs.")
  return(genepairs)
} 

align_gene_CCL <- function(gene1, gene2, mix1, mix2, verbose=TRUE) {
  df1 <- mix1[[gene1]]
  df2 <- mix2[[gene2]]

  if (is.null(df1) || is.null(df2)) {
    if (verbose) message("Skipping pair: gene not found in input mix lists: ", gene1, ", ", gene2)
    return(NULL)
  }
  
  # get unique cell-line IDs (column 1) in each
  lines1 <- unique(df1[[1]])
  lines2 <- unique(df2[[1]])
  # find the intersection
  shared <- intersect(lines1, lines2)
  if (length(shared) == 0) {
    if (verbose) message("No shared cell lines for ", gene1, " & ", gene2)
    return(NULL)
  }
  # subset both data frames to only those shared lines
  df1f <- df1[df1[[1]] %in% shared, , drop = FALSE]
  df2f <- df2[df2[[1]] %in% shared, , drop = FALSE]
  list(cluster1 = df1f, cluster2 = df2f,shared = shared)
}

OneGMMBinomial = function(gene1, gene2, mixModelClusters1, prop_cluster, user_tmp_map=NULL, directionality, effectsize =T) {
  tmp_map <- constructContingencyTable(mixModelClusters1[[gene1]][, 3], mixModelClusters1[[gene2]][, 3])
  n_count <- sum(tmp_map)
  # the subset of prop_cluster for the current genes
  prop_gene1 <- prop_cluster[prop_cluster$gene == gene1, ]
  prop_gene2 <- prop_cluster[prop_cluster$gene == gene2, ] 
  # one minus for 1-mid effect size calculation
  .binomialtest(gene1, gene2, prop_gene1, prop_gene2,
                tmp_map, n_count, directionality, effectsize, invertMidpoint=T)
}


TwoGMMBinomial = function(gene1, gene2, mixModelClusters1, mixModelClusters2, prop_cluster1, prop_cluster2, user_tmp_map=NULL, directionality, effectsize =T, verbose=TRUE) {
  filtered_lists = align_gene_CCL(gene1, gene2, mixModelClusters1, mixModelClusters2, verbose)
  if (is.null(filtered_lists)) return(NULL)
  v1 <- filtered_lists$cluster1[, 3]
  v2  <- filtered_lists$cluster2[, 3]
  tmp_map <- constructContingencyTable(v1, v2)
  n_count <- sum(tmp_map)
  # get the subset of prop_cluster dataframes for the current genes
  prop1_gene <- prop_cluster1[prop_cluster1$gene == gene1, ]
  prop2_gene <- prop_cluster2[prop_cluster2$gene == gene2, ]
  .binomialtest(gene1, gene2, prop1_gene, prop2_gene,
                tmp_map, n_count, directionality, effectsize)
}

.binomialtest <- function(gene1, gene2, prop_cluster1, prop_cluster2, tmp_map, n_count, directionality, effectsize,
                          invertMidpoint =F) {
  num_rows_tmp <- nrow(tmp_map)
  num_cols_tmp <- ncol(tmp_map)
  geneCluster_index <- data.frame()
  
  for (i in 1:num_rows_tmp) { # rows of tmp_map correspond to clusters of gene1
    for (j in 1:num_cols_tmp) { # columns of tmp_map correspond to clusters of gene2
      
      actual_count <- tmp_map[i, j]
      
      # Find the proportion for the current cluster combination
      prop_row_gene1 <- prop_cluster1[prop_cluster1$gene == gene1 & prop_cluster1$cluster == i, ]
      prop_row_gene2 <- prop_cluster2[prop_cluster2$gene == gene2 & prop_cluster2$cluster == j, ]
      
      if (nrow(prop_row_gene1) == 0 || nrow(prop_row_gene2) == 0)
        next
      
      MM1_proportion <- prop_row_gene1$prop
      MM2_proportion <- prop_row_gene2$prop
      exp_count <- (MM1_proportion * MM2_proportion) * n_count
      expected_cell_value = exp_count / n_count
      Cluster_Combination <- paste0(i, "_x_", j)
      TotalCluster_Combinations = num_rows_tmp * num_cols_tmp
      
      tmp <- binom.test(actual_count, n_count, p = expected_cell_value, alternative = directionality, conf.level = 0.95)
      
      num_clusters_gene1_data <- max(prop_cluster1$cluster[prop_cluster1$gene == gene1])
      num_clusters_gene2_data <- max(prop_cluster2$cluster[prop_cluster2$gene == gene2])
      
      if (effectsize) {
        mid = (tmp$conf.int[1] + tmp$conf.int[2]) / 2
        if (invertMidpoint) {
          mid = 1 - mid
        }
        tmp_results <- data.frame(Gene1 = gene1, Gene2 = gene2,
                                  Gene1_ClusterData = num_clusters_gene1_data, Gene2_ClusterData = num_clusters_gene2_data,
                                  Actual_Count = actual_count, Expected_Count = exp_count, TotalCluster_Combinations = TotalCluster_Combinations,
                                  Cluster_Combination = Cluster_Combination, Sample_Count = n_count, p_value = tmp$p.value, Effect_Size=mid)
      } else {
        tmp_results <- data.frame(Gene1 = gene1, Gene2 = gene2,
                                  Gene1_ClusterData = num_clusters_gene1_data, Gene2_ClusterData = num_clusters_gene2_data,
                                  Actual_Count = actual_count, Expected_Count = exp_count, TotalCluster_Combinations = TotalCluster_Combinations,
                                  Cluster_Combination = Cluster_Combination, Sample_Count = n_count, p_value = tmp$p.value)
      }
      geneCluster_index <- rbind(geneCluster_index, tmp_results)
    }
  }
  return(geneCluster_index)
}

constructContingencyTable <- function(mixModelCluster1, mixModelCluster2 = NULL) {
  # for 1D analysis, set mixModelCluster2 to mixModelCluster1
  if (is.null(mixModelCluster2)) {
    mixModelCluster2 <- mixModelCluster1
  }
  vals1 <- mixModelCluster1[!is.na(mixModelCluster1)]
  vals2 <- mixModelCluster2[!is.na(mixModelCluster2)]
  mixModelCluster1 = as.factor(mixModelCluster1)
  mixModelCluster2 = as.factor(mixModelCluster2)
  unique_values_cluster1 <- as.numeric(levels(mixModelCluster1))
  unique_values_cluster2 <-  as.numeric(levels(mixModelCluster2))
  unique_values_cluster1 <- max(as.numeric(unique_values_cluster1))
  unique_values_cluster2 <- max(as.numeric(unique_values_cluster2))
  tmp_map <- matrix(data = 0, nrow = unique_values_cluster1, ncol = unique_values_cluster2)
  for (j in 1:unique_values_cluster1) {
    for (k in 1:unique_values_cluster2) {
      tmp_map[j, k] <- sum(mixModelCluster1 == j & mixModelCluster2 == k)
    }
  }
  return(tmp_map)
}


GeneClusterProportions <- function(mixModel, gene_index, cores, all_patterns = TRUE) { 
  tryCatch({
    cl <- as.integer(mixModel[[gene_index]][, 3])
    cl <- cl[!is.na(cl)]
    n_count <- length(cl)
    if (n_count == 0) return(NULL)
    # the maximum cluster value 
    maxClusterVal <- max(cl, na.rm = TRUE)
    # proportion of samples, calculated for every cluster from 1 to maxClusterVal
    props <- tabulate(cl, nbins = maxClusterVal) / n_count
    
    df <- data.frame(gene = names(mixModel)[gene_index],
                     cluster = seq_len(maxClusterVal),
                     prop = props,
                     stringsAsFactors = FALSE)
    
    # for SL-only mode, keep cluster 1 row only
    if (!all_patterns) {
      df <- df[df$cluster == 1, , drop = FALSE]
    }
    df
  }, error = function(e) {
    NULL
  })
}

# to get common samples/cell lines when performing bidirectional analysis
GMM_CCL <- function(list1, list2 = NULL) {
  list1 <- prefilterMixModelClusters(list1)
  
  if (!is.null(list2)) {
    list2 <- prefilterMixModelClusters(list2)
    
    # unique cell lines for each gene in both lists
    cell_lines_list1 <- lapply(list1, function(df) unique(df[, 1]))
    cell_lines_list2 <- lapply(list2, function(df) unique(df[, 1]))
    
    # find common cell lines across all genes
    common_cell_lines <- Reduce(intersect, c(cell_lines_list1, cell_lines_list2))
    
    # filter each gene's data frame to include only the common cell lines
    list1_filtered <- lapply(list1, function(df) df[df[, 1] %in% common_cell_lines, , drop = FALSE])
    list2_filtered <- lapply(list2, function(df) df[df[, 1] %in% common_cell_lines, , drop = FALSE])
    
    # rename and return both
    names(list1_filtered) <- names(list1)
    names(list2_filtered) <- names(list2)
    return(list(list1 = list1_filtered, list2 = list2_filtered, samples = common_cell_lines))
    
  } else {
    # one-directional mode
    cell_lines_list <- lapply(list1, function(df) unique(df[, 1]))
    common_cell_lines <- Reduce(intersect, cell_lines_list)
    
    list1_filtered <- lapply(list1, function(df) df[df[, 1] %in% common_cell_lines, , drop = FALSE])
    names(list1_filtered) <- names(list1)
    
    return(list(list1 = list1_filtered, samples = common_cell_lines))
  }
}
mts_pickAxisTransform <- function(values,
                                  cuts = NULL,
                                  method = c("sqrt", "pseudo_log"),
                                  skew_threshold = 0.20,
                                  ratio_threshold = 12,
                                  use_skew  = TRUE,
                                  use_ratio = TRUE) {

  method <- match.arg(method)
  v <- values[is.finite(values)]

  if (length(v) < 3L) {
    out <- "linear"; attr(out, "diagnostics") <- list(reason = "too few finite values"); return(out)
  }

  mn <- min(v); mx <- max(v); rng <- mx - mn
  if (!is.finite(rng) || rng <= 0) {
    out <- "linear"; attr(out, "diagnostics") <- list(reason = "zero range"); return(out)
  }

  med          <- stats::median(v)
  has_negative <- mn < 0

  ## median position within the range
  med_pos <- (med - mn) / rng
  fired_skew  <- isTRUE(use_skew) && med_pos < skew_threshold

  ## backstop: dynamic-range ratio (only meaningful for non-negative data)
  ratio <- if (mn >= 0 && med > 0) mx / med else NA_real_
  fired_ratio <- isTRUE(use_ratio) && is.finite(ratio) && ratio > ratio_threshold

  squished <- fired_skew || fired_ratio

  if (squished) {
    decision <- if (method == "sqrt" && has_negative) "pseudo_log" else method
  } else {
    decision <- "linear"
  }

  attr(decision, "diagnostics") <- list(
    med_pos = med_pos, skew_fired = fired_skew, ratio = ratio,
    ratio_fired = fired_ratio, has_negative = has_negative,
    requested = method, decision = as.character(decision)
  )
  decision
}

## ggplot scale for one axis given a transform.
.mts_axis_scale <- function(axis = c("x", "y"),
                            transform = "linear",
                            limits = NULL,
                            breaks = ggplot2::waiver(),
                            expand = ggplot2::waiver(),
                            pseudo_log_sigma = 1) {
  axis <- match.arg(axis)
  oob  <- scales::oob_keep
  if (identical(transform, "sqrt")) {
    if (axis == "x") ggplot2::scale_x_sqrt(limits = limits, breaks = breaks, expand = expand, oob = oob)
    else             ggplot2::scale_y_sqrt(limits = limits, breaks = breaks, expand = expand, oob = oob)
  } else if (identical(transform, "pseudo_log")) {
    tr <- scales::pseudo_log_trans(sigma = pseudo_log_sigma, base = 10)
    if (axis == "x") ggplot2::scale_x_continuous(trans = tr, limits = limits, breaks = breaks, expand = expand, oob = oob)
    else             ggplot2::scale_y_continuous(trans = tr, limits = limits, breaks = breaks, expand = expand, oob = oob)
  } else {
    if (axis == "x") ggplot2::scale_x_continuous(limits = limits, breaks = breaks, expand = expand, oob = oob)
    else             ggplot2::scale_y_continuous(limits = limits, breaks = breaks, expand = expand, oob = oob)
  }
}

.mts_density_support <- function(values, modes) {
  v0 <- values[is.finite(values)]
  if (length(v0) < 2L) return(range(v0))
  lo <- min(v0); hi <- max(v0)
  for (m in unique(modes)) {
    vv <- values[modes == m]; vv <- vv[is.finite(vv)]
    if (length(vv) >= 2L && diff(range(vv)) > 0) {
      dd <- tryCatch(stats::density(vv), error = function(e) NULL)
      if (!is.null(dd)) { lo <- min(lo, min(dd$x)); hi <- max(hi, max(dd$x)) }
    }
  }
  c(lo, hi)
}

.mts_combine_breaks <- function(cuts, lims, transform, tol_frac = 0.03) {
  lo <- if (identical(transform, "sqrt")) max(0, lims[1]) else lims[1]
  base <- scales::breaks_extended()(c(lo, lims[2]))
  base <- base[is.finite(base) & base >= lims[1] & base <= lims[2]]
  cuts <- cuts[is.finite(cuts) & cuts >= lims[1] & cuts <= lims[2]]
  if (length(cuts)) {
    tol <- tol_frac * diff(lims)
    base <- base[!vapply(base, function(b) any(abs(b - cuts) < tol), logical(1))]
  }
  out <- sort(unique(round(c(base, cuts), 1)))
  if (length(out) == 0) ggplot2::waiver() else out
}

.mts_axis_limits <- function(values, modes, transform, override = NULL) {
  if (!is.null(override)) return(override)
  s <- .mts_density_support(values, modes)
  if (identical(transform, "sqrt")) c(0, s[2]) else s
}