cell typing using negative binomial and f-test pvalue for error detection

NAMESPACE

@@ -19,17 +19,21 @@ export(getScoreCellType_gene)
 export(get_baselineCT)
 export(groupTranscripts_Delaunay)
 export(groupTranscripts_dbscan)
+export(lldist)
 export(neighborhood_for_resegment)
 export(neighborhood_for_resegment_spatstat)
+export(numCores)
 export(plotSpatialScoreMultiCells)
 export(prepSMI_for_fastReseg)
+export(quick_celltype)
 export(scoreGenesInRef)
 export(score_cell_segmentation_error)
 export(update_transDF_ResegActions)
 importFrom(Giotto,createSpatialNetwork)
 importFrom(Giotto,pDataDT)
 importFrom(Matrix,colMeans)
 importFrom(Matrix,rowSums)
+importFrom(RcppEigen,fastLmPure)
 importFrom(concaveman,concaveman)
 importFrom(data.table,as.data.table)
 importFrom(data.table,setDT)
@@ -61,3 +65,4 @@ importFrom(spatstat.geom,pp3)
 importFrom(spatstat.geom,ppp)
 importFrom(spatstat.geom,subset.pp3)
 importFrom(spatstat.geom,subset.ppp)
+importFrom(stats,dnbinom)

R/flag_errors.R

@@ -11,11 +11,12 @@
 #' \enumerate{
 #'    \item{cell_ID, cell id}
 #'    \item{transcript_num, number of transcripts in given cell}
-#'    \item{modAlt_rsq, summary(mod_alternative)$r.squared}
-#'    \item{lrtest_ChiSq, lrtest chi-squared value}
-#'    \item{lrtest_Pr, lrtest probability larger than chi-squared value, p-value}
+#'    \item{modAlt_rsq, the root mean square for residual of the alternative model }
+#'    \item{lm_Fstat, the F-test statistic of the alternative model against null model}
+#'    \item{lm_Pvalue, the p.value calculated from the F-test statstic}
 #' }
-#' @details For tLLRv2 score of transcripts within each cell,  run a quadratic model: mod_alternative = lm(tLLRv2 ~ x + y + x2 + y2 +xy) for 2D,  lm(tLLRv2 ~ x + y + z + x2 + y2 +z2 +xy + xz + yz) for 3D and a null model: mod_null = lm(tLLRv2 ~ 1); then run lmtest::lrtest(mod_alternative, mod_null). Return statistics for mod_alternative$fitted.values (standard deviation and minimal value), summary(mod_alternative)$r.squared and as well as lrtest chi-squared value.  
+#' @details For tLLRv2 score of transcripts within each cell,  run a quadratic model: mod_alternative = lm(tLLRv2 ~ x + y + x2 + y2 +xy) for 2D,  lm(tLLRv2 ~ x + y + z + x2 + y2 +z2 +xy + xz + yz) for 3D. Return the root mean square of residual after fitting, the F statistics and p.value of alternative model against null model. 
+#' @importFrom RcppEigen fastLmPure
 #' @export
 score_cell_segmentation_error <- function(chosen_cells, transcript_df, 
                                           cellID_coln = "CellId", 
@@ -66,10 +67,14 @@ score_cell_segmentation_error <- function(chosen_cells, transcript_df,
     # lm(tLLRv2 ~ x + y + x2 + y2 + xy) for 2D,  lm(tLLRv2 ~ x + y + z + x2 + y2 +z2 +xy + xz + yz) for 3D
     if(d2_or_d3 ==2){
       colnames(coord_df) <- c('cell_ID','score','x','y')
-      mod_formula <- 'score ~ x + y + x2 + y2 + xy'
+      # mod_formula <- 'score ~ x + y + x2 + y2 + xy'
+
+      colns_to_regress <- c('x','y','x2','y2','xy')
     } else {
       colnames(coord_df) <- c('cell_ID','score','x','y','z')
-      mod_formula <- 'score ~ x + y + z + x2 + y2 +z2 +xy + xz + yz'
+      # mod_formula <- 'score ~ x + y + z + x2 + y2 +z2 +xy + xz + yz'
+
+      colns_to_regress <- c('x','y','z','x2','y2','z2','xy','xz','yz')
     }
 
     coord_df[['x2']] <- coord_df[['x']]^2
@@ -83,22 +88,62 @@ score_cell_segmentation_error <- function(chosen_cells, transcript_df,
     }
 
 
-    # (3) perform lm and lrtest by group
-    my_fun <- function(data){
-      # null linear model, lm(tLLRv2 ~ 1)
-      mod_null <- lm(score~1, data = data)
+    # (3) perform lm by group
+    # # lrtest to evaluate spatial dependency 
+    # my_fun <- function(data){
+    #   # null linear model, lm(tLLRv2 ~ 1)
+    #   mod_null <- lm(score~1, data = data)
+    #   
+    #   # spatial quadratic model
+    #   # lm(tLLRv2 ~ x + y + x2 + y2 + xy) for 2D,  lm(tLLRv2 ~ x + y + z + x2 + y2 +z2 +xy + xz + yz) for 3D
+    #   mod_alternative <- lm(as.formula(mod_formula), data = data)
+    #   
+    #   #likelihood ratio test of nested model
+    #   lrtest_res <- lmtest::lrtest(mod_alternative, mod_null)
+    #   
+    #   outputs <- data.frame(transcript_num = nrow(data), 
+    #                         modAlt_rsq = summary(mod_alternative)$r.squared,
+    #                         lrtest_ChiSq = lrtest_res$Chisq[2],
+    #                         lrtest_Pr= lrtest_res$`Pr(>Chisq)`[2])
+    #   
+    #   return(outputs)
+    # }
+
+    # fastLmPure and F-statistics to evaluate spatial dependency 
+    fstat <- function(flmp,y){
+      n <- length(flmp$residuals)
+      sumsquares_residual_h0 <- var(y)*(n-1)
+      sumsquares_residual_h1 <- sum(flmp$residuals^2 )
+      p1 <- length(flmp$coefficients) 
+      p0 <- 1
+      fstat <- 
+        ((sumsquares_residual_h0 - sumsquares_residual_h1)/(p1-p0)) / 
+        ((sumsquares_residual_h1)/(n-p1))  
+      p.value <- pf(fstat, p1-p0, n-p1, lower.tail=FALSE)
 
-      # spatial quadratic model
+      return(list(fstat = c("value" = fstat, "numdf" = p1-p0, "dendf" = n-p1)
+                  ,p.value = p.value
+      ))
+    }
+
+    my_fun <- function(data){
+     # linear regression using spatial quadratic model
       # lm(tLLRv2 ~ x + y + x2 + y2 + xy) for 2D,  lm(tLLRv2 ~ x + y + z + x2 + y2 +z2 +xy + xz + yz) for 3D
-      mod_alternative <- lm(as.formula(mod_formula), data = data)
 
-      #likelihood ratio test of nested model
-      lrtest_res <- lmtest::lrtest(mod_alternative, mod_null)
+      ## note that have to add the column of 1's for the intercept
+      flmp <- RcppEigen::fastLmPure(y=data$score, 
+                                    X=as.matrix(cbind(rep(1, nrow(data)), 
+                                                      data[,.SD, .SDcols = colns_to_regress]))
+                                    ) 
+
+      # f-statistics and p.value
+      fstatistic <- fstat(flmp, data$score)
 
       outputs <- data.frame(transcript_num = nrow(data), 
-                            modAlt_rsq = summary(mod_alternative)$r.squared,
-                            lrtest_ChiSq = lrtest_res$Chisq[2],
-                            lrtest_Pr = lrtest_res$`Pr(>Chisq)`[2])
+                            modAlt_rsq = flmp$s,
+                            lm_Fstats = fstatistic[['fstat']][['value']],
+                            lm_Pvalue = fstatistic[['p.value']])
+
       return(outputs)
     }
 

R/get_baseline.R

@@ -171,6 +171,7 @@ choose_distance_cutoff <- function(transcript_df,
 #' @param refProfiles A matrix of cluster profiles, genes X clusters
 #' @param counts Counts matrix, cells X genes.
 #' @param clust Vector of cluster assignments for each cell in `counts`, default = NULL to automatically assign the cell cluster for each cell based on maximum transcript score  
+#' @param celltype_method use either `LogLikeRatio` or `NegBinomial` method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio)
 #' @return a list
 #' \enumerate{
 #'    \item{span_score, a matrix of average transcript tLLR score per molecule per cell for 22 distinct cell types in rows, percentile at (0%, 25%, 50%, 75%, 100%) in columns}
@@ -188,7 +189,11 @@ choose_distance_cutoff <- function(transcript_df,
 #' @export
 get_baselineCT <- function(refProfiles, 
                            counts, 
-                           clust = NULL){
+                           clust = NULL,
+                           celltype_method = 'LogLikeRatio'){
+
+  celltype_method <- match.arg(celltype_method, c('LogLikeRatio', 'NegBinomial'))
+
   # get common genes
   common_genes <- intersect(rownames(refProfiles), colnames(counts))
 
@@ -236,30 +241,56 @@ get_baselineCT <- function(refProfiles,
   counts <- as.matrix(counts)[, common_genes]
 
   # get score matrix based on refProfiles for each gene and cell ----
-  # replace zero in mean profiles with 1E-5
-  refProfiles <- pmax(refProfiles, 1e-5)
-  # tLL score
-  transcript_loglik <- scoreGenesInRef(genes = common_genes, ref_profiles = refProfiles)
-  # tLLR score, re-center on maximum per row/transcript
-  tmp_max <- apply(transcript_loglik, 1, max)
-  tLLRv2_geneMatrix <- sweep(transcript_loglik, 1, tmp_max, '-')
-  rm(tmp_max, transcript_loglik)
-
-  # get cell x cell-cluster score matrix = counts (cell x gene) %*% tLLR_score (gene x cell-cluster)
-  tLLRv2_cellMatrix <- counts %*% tLLRv2_geneMatrix
+  # replace zero in mean profiles with 1E-8
+  refProfiles <- pmax(refProfiles, 1e-8)
 
-  # assign cell type for each cell if not provided ----
-  if(is.null(clust)){
-    message('Perform cluster assignment based on maximum transcript score given the provided `refProfiles`.')
+  if(celltype_method == 'LogLikeRatio'){
+    # tLL score
+    transcript_loglik <- scoreGenesInRef(genes = common_genes, ref_profiles = refProfiles)
+    # tLLR score, re-center on maximum per row/transcript
+    tmp_max <- apply(transcript_loglik, 1, max)
+    tLLRv2_geneMatrix <- sweep(transcript_loglik, 1, tmp_max, '-')
+    rm(tmp_max, transcript_loglik)
+
+    # get cell x cell-cluster score matrix = counts (cell x gene) %*% tLLR_score (gene x cell-cluster)
+    tLLRv2_cellMatrix <- counts %*% tLLRv2_geneMatrix
+
+    # assign cell type for each cell if not provided
+    if(is.null(clust)){
+      message('Perform cluster assignment based on maximum transcript score given the provided `refProfiles`.')
+
+      # assign cell type based on max values
+      max_idx_1st <- max.col(tLLRv2_cellMatrix, ties.method="first")
+      clust <- colnames(tLLRv2_cellMatrix)[max_idx_1st]
 
-    # assign cell type based on max values
-    max_idx_1st <- max.col(tLLRv2_cellMatrix, ties.method="first")
-    clust <- colnames(tLLRv2_cellMatrix)[max_idx_1st]
+      rm(max_idx_1st)
+    }
+
+
+  } else if (celltype_method == 'NegBinomial'){
+    # get logliks for cell under all cell types
+    nb_res <- quick_celltype(counts, bg = 0.01, 
+                             reference_profiles = refProfiles, 
+                             align_genes = FALSE)
+
+    # per cell logliks for all cells, cell x cell-cluster score matrix, exclude cells of zero count
+    tLLRv2_cellMatrix <- nb_res[['logliks']][1: (length(nb_res[['clust']]) - length(nb_res[['zeroCells']])), ]
+
+    # assign cell type for each cell if not provided 
+    if(is.null(clust)){
+      message('Perform cluster assignment based on negative binomial model given the provided `refProfiles`.')
+      clust <- nb_res[['clust']][1: (length(nb_res[['clust']]) - length(nb_res[['zeroCells']]))]
+
+    }
+    rm(nb_res)
 
-    common_celltypes <- unique(clust)
-    rm(max_idx_1st)
+
+  } else {
+    stop(sprintf('The provided `celltype_method` = `%s` is not supported.', celltype_method))
   }
 
+  common_celltypes <- unique(clust)
+
   # get transcript number quantile profile ---
   all_transNum <- rowSums(counts)
   span_transNum_CellType <- tapply(all_transNum, 
@@ -274,6 +305,7 @@ get_baselineCT <- function(refProfiles,
     rowidx <- which(clust == each_celltype)
     all_tLLRv2[rowidx] <- tLLRv2_cellMatrix[rowidx, each_celltype]
   }
+
   # normalized by transcript number to get per molecule transcript score for each cell
   all_tLLRv2 <- all_tLLRv2/all_transNum
   span_tLLRv2_CellType <- tapply(all_tLLRv2,