speed up the deviance_residual_transform
mydeviance_residual_transform <- function(count_mat){ resid_mat <- matrix(nrow = nrow(count_mat), ncol = ncol(count_mat)) #row is cell,coll is gene row_sum <- rowSums(count_mat) #每个细胞总read col_sum <- colSums(count_mat) #每个基因总read total_sum <- sum(count_mat)
for (i in 1:nrow(count_mat)){ U_i <- row_sum[i]col_sum/total_sum cross_prod <- 2count_mat[i,]log((count_mat[i,]+1e-6)/(U_i+1e-6)) + 2(row_sum[i] - count_mat[i,])*log((row_sum[i]-count_mat[i,])/(row_sum[i]-U_i)) l <- sign(count_mat[i,]-U_i) * sqrt(cross_prod) resid_mat[i, ] <- l } return(resid_mat) }
dev_res <- mydeviance_residual_transform(t(as.matrix(combined_obj@assays$RNA@counts[1:1000,1:1000]))) dev_res2 <- deviance_residual_transform(t(as.matrix(combined_obj@assays$RNA@counts[1:1000,1:1000])))
max(dev_res -dev_res2) [1] 1.856293e-13