Need to have a faster implementation of getting deviances for splicing and stability events/features

[Arshammik]

Main

On the current pipeline, we use a glm fitting with a binomial family to model m1 as success and m2 as failure. Here is the original code we used previously.

# the glm function 
.glm_fit <- function(nz_indx, nz_val) {  
  
  temp_m1_m2 <- rep(0, ncol(temp_X))
  temp_m1_m2[nz_indx] <- nz_val
  
  #fit <- glm(cbind(temp_m1_m2[m1_indexes], temp_m1_m2[m2_indexes]) ~ 1, family = "binomial")
  fit <- tryCatch({
    glm(cbind(temp_m1_m2[m1_indexes], temp_m1_m2[m2_indexes]) ~ 1, family = "binomial")
  }, error = function(e) {
    message("Error in glm: ", e$message, ", due to not convargence, we put zero as deviance for this case.")
    return(NULL)
  }) 
  if (is.null(fit)) {
    return(0)  # Return NA if glm failed
  }
  return(fit$deviance)
 }

To do

1. For now, since we extract fit$deviance, need to do is to implement the deviance in a C++ script and benchmark with the previous results.
2. Then, create an R script and matching cpp function using Rcpp to make getting the high variable events easier. I suggest a function for splicing and stability.
3. (Optional) If the implementation is successful, we can try to build our native feature selection for gene expression and not use Seuart's Find VariableFeatures by using the vst method.

Thanks!

[Arshammik]

Mathematical concept

To solve this problem, we need insight into the deviance formula in binomial distribution. The formula for deviance in a binomial distribution can be defined as:

$$Deviance = 2[ylog(\frac{y}{n\hat{p}}) +(n-y)log(\frac{n-y}{n(1-\hat{p})})]$$

Where $y$ is the number of successes, and $n$ is the number of all trials (successes + failures). The MLE Estimate which is the $\hat{p}$ defines as $sum(y) / sum(n)$

Implementation

We can write the code in C++ for much faster implementation to calculate the binomial deviance contribution for each cell over columns (k) within each row (i). Here's the specific section of implemented code:

double dev_row = 0.0;
for (int k = 0; k < n_cols; k++){
  double y_val = M1(i, k); // Successes in cell (i,k)
  double f_val = M2(i, k); // Failures in cell (i,k)
  double n_i = y_val + f_val; // Total trials in cell (i,k)

The p_hat is defined as sum_y / sum_n in the following code. We can also skip the cells with zero trails.

if (n_i <= 0) continue;

double p_hat = sum_y / sum_n;
     // If p_hat is on the boundary, set deviance to 0.
     if (p_hat <= 0 || p_hat >= 1) {
       dev[i] = 0;
       continue;
     }

To compute the deviance contribution for Success Term, we should implement

$$ylog\frac{y}{n\hat{p}}$$

The related part in the cpp script is:

  double term = 0.0;
  // Term for successes (y_val):
  if (y_val > 0)
    term += y_val * std::log(y_val / (n_i * p_hat));

For the second term, to compute the deviance contribution for the Failure Term, which is

$$(n-y)log(\frac{n-y}{n(1-\hat{p})}$$

we used this script for C++ implementation

  // Term for failures (n_i - y_val):
  if ((n_i - y_val) > 0)
    term += (n_i - y_val) * std::log((n_i - y_val) / (n_i * (1 - p_hat)));

Lastly, we simply do the multiplication by two at the final lines of the script

dev_row += 2 * term; // 2 × [success term + failure term]

Initial demo tests indicate that the code is considerably faster and computes the same numbers with a minimal error rate for events with very small deviations (1e-10 <).

Benchmarking

Old method (glm in R)

We test this implementation on a test data set containing 5972 cells with 901198 events. Initially, we ran the previous method using tictoc code execution timing functions to monitor the timing. It is worth mentioning that the previous method used multithreading in data.table while the current one is in single-thread computation. Here is the original code for the previous method:

filter <- which(meta[, ID] == lib) 
m1_indexes <- 1:length(filter)
m2_indexes <- (length(filter) + 1):(2*length(filter))

# We convert it to matrix for getting vectors without names and faster processing
temp_X <- cbind(M1[,filter],M2[,filter])

# creating a summary data table
temp_X_summ <- summary(temp_X) |> as.data.table()

# doing a glm fit based on group by in data table
tic()
temp_dev <- temp_X_summ[, .(dev = .SD[, .glm_fit(j, x)]) ,i]
toc()

The total time taken for the previous method was around one hour.

> tic()
> temp_dev <- temp_X_summ[, .(dev = .SD[, .glm_fit(j, x)]) ,i]
Processed 309575 groups out of 309575. 100% done. Time elapsed: 3585s. ETA: 0s...
There were 50 or more warnings (use warnings() to see the first 50)
> toc()
3585.56 sec elapsed

New method

The new implementation demonstrates significantly enhanced performance metrics. Below is a benchmark analysis utilizing the same dataset with the newly introduced function called multigedi_get_deviance .

tic()
deviance_values <- multigedi_get_deviance(m1_matrix = temp_X[, 1:5972], m2_matrix = temp_X[,5973:11944])
toc()

With the console log

There are  1 libraries detected... 
Calculating the deviances for sample SRR26455875 has been completed!
24.509 sec elapsed

We conducted a Pearson correlation analysis between the new and old datasets to ensure consistency in our results.

> sum_deviances_cpp <- data.table(rn = deviance_values$events, deviance = deviance_values$sum_deviance)
> sum_deviances_dt <- as.data.table(sum_deviances, keep.rownames = TRUE)
> sum_deviances_dt <- sum_deviances_dt[rn %in% sum_deviances_cpp$rn, ]
> dt <- merge(sum_deviances_dt, sum_deviances_cpp, by = "rn")
> cor(dt$sum_deviances, dt$deviance)
[1] 1

Lastly, here is the smoothScatter(log(dt$sum_deviances + 1), log(dt$deviance + 1)) for a visual inspection as well.

[Arshammik]

Yet another example with much larger data

previous <- merge(previous, deviance_sum, by = "events")

cor(previous$R_sum_deviance, previous$sum_deviance)
[1] 0.9961936

and the related scatter plot