From ab95e0a18b956d17b0ec90ba4ca8fab53e06bc0a Mon Sep 17 00:00:00 2001 From: lidanwu Date: Fri, 21 Oct 2022 00:25:26 +0000 Subject: [PATCH 01/10] use negbinomial for cell typing --- R/get_baseline.R | 79 ++++++--- R/neighborhood_spatstat.R | 85 +++++++--- R/preprocessing.R | 9 +- R/quick_celltyping.R | 147 +++++++++++++++++ R/update_segments.R | 77 ++++++--- R/wrapper_FastReseg.R | 201 ++++++++++++++++------- vignettes/0_flagErrorOnly_on_SMIobject.R | 9 +- vignettes/1_fastReseg_on_SMIobject.R | 10 +- 8 files changed, 490 insertions(+), 127 deletions(-) create mode 100644 R/quick_celltyping.R diff --git a/R/get_baseline.R b/R/get_baseline.R index 1379a11..ab34a67 100644 --- a/R/get_baseline.R +++ b/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,59 @@ 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) + # replace zero in mean profiles with 1E-8 + refProfiles <- pmax(refProfiles, 1e-8) - # 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] - - common_celltypes <- unique(clust) - rm(max_idx_1st) + 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] + + common_celltypes <- unique(clust) + rm(max_idx_1st) + } + + + } else if (celltype_method == ' NegBinomial'){ + + # 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`.') + + nb_res <- quick_celltype(counts, bg = 0.01, + reference_profiles = refProfiles, + align_genes = FALSE) + + + clust <- nb_res[['clust']] + common_celltypes <- unique(clust) + + } + + # per cell logliks for all cells, cell x cell-cluster score matrix + tLLRv2_cellMatrix <- nb_res[['logliks']] + + + } else { + stop(sprintf('The provided `celltype_method` = `%s` is not supported.', celltype_method)) } + + # get transcript number quantile profile --- all_transNum <- rowSums(counts) span_transNum_CellType <- tapply(all_transNum, @@ -274,6 +308,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, diff --git a/R/neighborhood_spatstat.R b/R/neighborhood_spatstat.R index 9939753..5e22176 100644 --- a/R/neighborhood_spatstat.R +++ b/R/neighborhood_spatstat.R @@ -2,8 +2,9 @@ #' @title neighborhood_for_resegment_spatstat #' @description find neighbor cells with transcripts that are direct neighbor of chosen_cell, check tLLRv2 score under neighbor cell type, return neighborhood information #' @param chosen_cells the cell_ID of chosen cells need to be evaluate for re-segmentation -#' @param score_GeneMatrix the gene x cell-type matrix of log-like score of gene in each cell type -#' @param score_baseline a named vector of score baseline for all cell type listed in score_GeneMatrix +#' @param score_GeneMatrix the gene x cell-type matrix of log-like score of gene in each cell type, needed if using `LogLikeRatio` cell typing method (default = NULL) +#' @param refProfiles A matrix of cluster profiles, genes X clusters, needed if using `NegBionomial` cell typing method (default = NULL) +#' @param score_baseline a named vector of score baseline for all cell type listed in `score_GeneMatrix` or `refProfiles` #' @param neighbor_distance_xy maximum cell-to-cell distance in x, y between the center of query cells to the center of neighbor cells with direct contact, same unit as input spatial coordinate. Default = NULL to use the 2 times of average 2D cell diameter. #' @param distance_cutoff maximum molecule-to-molecule distance within connected transcript group, same unit as input spatial coordinate (default = 2.7 micron). #' If set to NULL, the pipeline would first randomly choose no more than 2500 cells from up to 10 random picked ROIs with search radius to be 5 times of `neighbor_distance_xy`, and then calculate the minimal molecular distance between picked cells. The pipeline would further use the 5 times of 90% quantile of minimal molecular distance as `distance_cutoff`. This calculation is slow and is not recommended for large transcript data.frame. @@ -13,6 +14,7 @@ #' @param transID_coln the column name of transcript_ID in transcript_df #' @param transGene_coln the column name of target or gene name in transcript_df #' @param transSpatLocs_coln the column name of 1st, 2nd, optional 3rd spatial dimension of each transcript in transcript_df +#' @param celltype_method use either `LogLikeRatio` or `NegBinomial` method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio) #' @importFrom spatstat.geom ppp subset.ppp nncross pp3 subset.pp3 #' @return a data.frame #' #' \enumerate{ @@ -28,7 +30,8 @@ #' @details Locate neighbor cells of each query cell firstly via cell-to-cell distance in 2D plane within neighbor_distance_xy, then via molecule-to-molecule 3D distance within distance_cutoff. If no neighbor cells found for query cell, use the cell id and cell type of query cell to fill in the columns for neighbor cells in returned data.frame #' @export neighborhood_for_resegment_spatstat <- function(chosen_cells = NULL, - score_GeneMatrix, + score_GeneMatrix = NULL, + refProfiles = NULL, score_baseline = NULL, neighbor_distance_xy = NULL, distance_cutoff = 2.7, @@ -37,7 +40,18 @@ neighborhood_for_resegment_spatstat <- function(chosen_cells = NULL, celltype_coln = "cell_type", transID_coln = "transcript_id", transGene_coln = "target", - transSpatLocs_coln = c('x','y','z')){ + transSpatLocs_coln = c('x','y','z'), + celltype_method = 'LogLikeRatio'){ + + celltype_method <- match.arg(celltype_method, c('LogLikeRatio', 'NegBinomial')) + + if(celltype_method == 'LogLikeRatio' & is.null(score_GeneMatrix)){ + stop("Must provided `score_GeneMatrix` when using log-likelihood ratio based cell typing method.") + } + + if(celltype_method == 'NegBinomial' & is.null(refProfiles)){ + stop("Must provided `refProfiles` when using negative binomial cell typing method.") + } if(is.null(chosen_cells)){ stop("Must define chosen_cells to start resegmentation evaluation in neighborhood of each chosen cell.") @@ -45,12 +59,23 @@ neighborhood_for_resegment_spatstat <- function(chosen_cells = NULL, if(!is.null(score_baseline)){ if(!any(class(score_baseline) %in% c('numeric'))){ - stop("The provided score_baseline must be either a named numeric vector or NUll which would disable comparision with score_baseline. ") + stop("The provided `score_baseline` must be either a named numeric vector or NUll which would disable comparision with `score_baseline`. ") } - if(length(setdiff(colnames(score_GeneMatrix), names(score_baseline)))>0){ - stop(sprintf("The provided score_baseline is missing for the following cell types used in score_GeneMatrix: `%s`.", - paste0(setdiff(colnames(score_GeneMatrix), names(score_baseline)), collapse ="`, `"))) + + + if(celltype_method == 'LogLikeRatio'){ + if(length(setdiff(colnames(score_GeneMatrix), names(score_baseline)))>0){ + stop(sprintf("The provided `score_baseline` is missing for the following cell types used in `score_GeneMatrix`: `%s`.", + paste0(setdiff(colnames(score_GeneMatrix), names(score_baseline)), collapse ="`, `"))) + } + } else if (celltype_method =='NegBinomial'){ + if(length(setdiff(colnames(refProfiles), names(score_baseline)))>0){ + stop(sprintf("The provided `score_baseline` is missing for the following cell types used in `refProfiles`: `%s`.", + paste0(setdiff(colnames(refProfiles), names(score_baseline)), collapse ="`, `"))) + } } + + } @@ -100,8 +125,14 @@ neighborhood_for_resegment_spatstat <- function(chosen_cells = NULL, common_cells <- unique(transcript_df[[cellID_coln]]) # get common genes - common_genes <- intersect(rownames(score_GeneMatrix), - unique(transcript_df[[transGene_coln]])) + if(celltype_method == 'LogLikeRatio'){ + common_genes <- intersect(rownames(score_GeneMatrix), + unique(transcript_df[[transGene_coln]])) + } else if (celltype_method =='NegBinomial'){ + common_genes <- intersect(rownames(refProfiles), + unique(transcript_df[[transGene_coln]])) + } + message(sprintf("Found %d common cells and %d common genes among transcript_df, cell_networkDT, and score_GeneMatrix. ", length(common_cells), length(common_genes))) @@ -121,9 +152,13 @@ neighborhood_for_resegment_spatstat <- function(chosen_cells = NULL, } } - score_GeneMatrix <- score_GeneMatrix[common_genes, ] + transcript_df <- transcript_df[which(transcript_df[[cellID_coln]] %in% common_cells & transcript_df[[transGene_coln]] %in% common_genes), ] + if(celltype_method == 'LogLikeRatio'){ + score_GeneMatrix <- score_GeneMatrix[common_genes, ] + } + # get per cell dataframe and search range if neighbor_distance_xy = NULL perCell_coordM <- transcript_df[, list(CenterX = mean(get(transSpatLocs_coln[1])), CenterY = mean(get(transSpatLocs_coln[2])), @@ -163,6 +198,13 @@ neighborhood_for_resegment_spatstat <- function(chosen_cells = NULL, # data for chosen cells only chosen_transDF <- transcript_df[which(transcript_df[[cellID_coln]] %in% chosen_cells), ] + # get cell type and logliks for chosen cells only + if(celltype_method == 'NegBinomial'){ + exprMat <- reshape2::acast(chosen_transDF, as.formula(paste(cellID_coln, '~', transGene_coln)), length) + nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = TRUE) + rm(exprMat) + } + ## get molecular_distance_cutoff between neighbor cells from 10 randomly selected ROIs with 5* neighbor_distance_xy if(is.null(distance_cutoff)){ if(nrow(perCell_coordM)> 2500){ @@ -360,15 +402,19 @@ neighborhood_for_resegment_spatstat <- function(chosen_cells = NULL, } - # get score matrix for each transcript in query cell - cell_score <- score_GeneMatrix[query_df[[transGene_coln]], ] - if(nrow(query_df) >1 ){ - cell_score <- colSums(cell_score) + # get sum of score matrix for each transcript in query cell + if(celltype_method == 'LogLikeRatio'){ + cell_score <- score_GeneMatrix[query_df[[transGene_coln]], ] + if(nrow(query_df) >1 ){ + cell_score <- colSums(cell_score) + } + cell_score <- matrix(cell_score, nrow = 1, dimnames = list(each_cell, colnames(score_GeneMatrix))) + } else if (celltype_method =='NegBinomial'){ + cell_score <- nb_res[['logliks']][each_cell, , drop = FALSE] + } - cell_score <- matrix(cell_score, nrow = 1, dimnames = list(each_cell, colnames(score_GeneMatrix))) + max_idx_1st <- max.col(cell_score, ties.method="first") - - queryPerCell_df <- data.frame(CellId = each_cell, cell_type = query_df[[celltype_coln]][1], transcript_num = nrow(query_df), @@ -382,8 +428,9 @@ neighborhood_for_resegment_spatstat <- function(chosen_cells = NULL, neighbors_df <- df_subset[which(df_subset[[cellID_coln]] %in% directCell_neighbors),] neighbors_df <- as.data.frame(neighbors_df)[, c(cellID_coln, celltype_coln)] neighbors_df <- unique(neighbors_df) - neighbors_df[['score_under_neighbor']] <- cell_score[each_cell, neighbors_df[[celltype_coln]]]/nrow(query_df) + neighbors_df[['score_under_neighbor']] <- cell_score[each_cell, neighbors_df[[celltype_coln]]]/nrow(query_df) + # if score baseline is provided, compare the net score above baseline to choose the consistent neighbor cell types if(!is.null(score_baseline)){ neighbors_df[['baseline']] <- score_baseline[neighbors_df[[celltype_coln]]] diff --git a/R/preprocessing.R b/R/preprocessing.R index 512ffc3..4bb2761 100644 --- a/R/preprocessing.R +++ b/R/preprocessing.R @@ -8,6 +8,7 @@ #' @param removeUnpaired flag to remove FOVs with unpaired target call files and fov position information; default = FALSE, to stop processing when missing target call files #' @param blacklist_genes a vector of genes to be excluded from reference profile estimation (default = NULL) #' @param pixel_size the micrometer size of image pixel listed in `Width` and `Height` dimension of each cell stored in `cell_metadata` of the existing SMI object (default = 0.18) +#' @param celltype_method use either `LogLikeRatio` or `NegBinomial` method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio) #' @return a list #' \describe{ #' \item{counts}{a cells X genes count matrix for entire dataset, stored in SMI object.} @@ -42,7 +43,11 @@ prepSMI_for_fastReseg <- function(path_to_SMIobject, cellClus_to_exclude = NULL, removeUnpaired = FALSE, blacklist_genes = NULL, - pixel_size = 0.18){ + pixel_size = 0.18, + celltype_method = 'LogLikeRatio'){ + + celltype_method <- match.arg(celltype_method, c('LogLikeRatio', 'NegBinomial')) + #### (1) prepare sample annotation file ---- ## check config_loading colns_to_use <- c('folderpathColumn','slidefoldersColumn','slidenameColumn','votedfoldersColumn','versionColumn') @@ -198,7 +203,7 @@ prepSMI_for_fastReseg <- function(path_to_SMIobject, # lowerCutoff_transNum, a named vector of 25% quantile of cluster-specific per molecule per cell transcript number, to be used as transcript number cutoff such that higher than the cutoff is required to keep query cell as it is # higherCutoff_transNum, a named vector of median value of cluster-specific per molecule per cell transcript number, to be used as transcript number cutoff such that lower than the cutoff is required to keep query cell as it is when there is neighbor cell of consistent cell type. # clust_used, a named vector of cluster assignments for each cell used in baseline calculation, cell_ID in `counts` as name - baselineData <- get_baselineCT(refProfiles = refProfiles, counts = counts, clust = as.character(clust)) + baselineData <- get_baselineCT(refProfiles = refProfiles, counts = counts, clust = as.character(clust), celltype_method = celltype_method) score_baseline <- baselineData[['score_baseline']] lowerCutoff_transNum <- baselineData[['lowerCutoff_transNum']] diff --git a/R/quick_celltyping.R b/R/quick_celltyping.R new file mode 100644 index 0000000..8fee5a1 --- /dev/null +++ b/R/quick_celltyping.R @@ -0,0 +1,147 @@ +#' @title lldist +#' @description calculate the log-likelihood of cell belong to certain cell cluster given the reference profile using negative binomial model +#' @param x a vector of a reference cell cluster +#' @param mat a cells x genes matrix of expression levels in all cells +#' @param bg a vector for background level of each cell (default: 0.01) +#' @param size the parameters for dnbinom function (default: 10) +#' @param digits the number of digits for rounding +#' +#' @importFrom Matrix rowSums +#' @importFrom stats dnbinom +#' @return a cell x cell_type matrix of the log-likelihood +lldist <- function(x, mat, bg = 0.01, size = 10, digits = 2) { + # convert to matrix form if only a vector was input: + if (is.vector(mat)) { + mat <- as(matrix(mat, nrow = 1), "dgCMatrix") + } else if (is.matrix(mat)) { + mat <- as(mat, "dgCMatrix") + } else if (!is(mat, "dgCMatrix")) { + errorMessage <- + sprintf( + "The `type` of parameter `mat` needs to be of one of dgCMatrix, vector, matrix, array, but is found to be of type %s", + class(mat) + ) + stop(errorMessage) + } + + # accept a single value of bg if input by user: + if (length(bg) == 1) { + bg <- rep(bg, nrow(mat)) + } + + # Check dimensions on bg and stop with informative error if not conformant + if (is.vector(bg)) { + if (!identical(length(bg), nrow(mat))) { + errorMessage <- sprintf("Dimensions of count matrix and background are not conformant.\nCount matrix rows: %d, length of bg: %d", + nrow(mat), length(bg)) + stop(errorMessage) + } + } + + # calc scaling factor to put y on the scale of x: + if (is.vector(bg)) { + bgsub <- apply(mat, 2, function(xx) xx - bg) + bgsub <- pmax(bgsub, 0) + } else { + bgsub <- pmax(mat - bg, 0) + } + + sum_of_x <- sum(x) + s <- Matrix::rowSums(bgsub) / sum_of_x + # override it if s is negative: + s[s <= 0] <- Matrix::rowSums(mat[s <= 0, , drop = FALSE]) / sum_of_x + + # yhat is a n_cells x n_genes matrix of expected values + yhat <- s %*% t(x) + + # log-likelihood, cell x gene matrix + res <- stats::dnbinom(x = as.matrix(mat), size = size, mu = yhat, log = TRUE) + # get sum of loglike per cell under current cell cluster + res <- rowSums(res) + + names(res) <- rownames(mat) + return(round(res, digits)) +} + + +#' Get number of cores for parallelized operations +#' +#' @return number of cores to use for mclapply +#' @export +numCores <- function() { + num_cores <- 1 + if (.Platform$OS.type == "unix") { + if (is.null(getOption("mc.cores"))) { + num_cores <- parallel::detectCores() - 2 + } else { + num_cores <- getOption("mc.cores") + } + + } + return(num_cores) +} + +#' @title quick_celltype +#' @description Classify cells based on reference profiles given the maximum log-likelihood calculated via negative binomial model +#' @param x Counts matrix (or dgCMatrix), cells x genes. +#' @param bg a vector for background level of each cell (default = 0.01) +#' @param reference_profiles Matrix of expression profiles of pre-defined clusters, genes x clusters. +#' @param nb_size The size parameter to assume for the NB distribution. +#' @param align_genes Logical, for whether to align the counts matrix and the reference_profiles by gene ID. +#' @return A list, with the following elements: +#' \enumerate{ +#' \item clust: a vector given cells' cluster assignments +#' \item logliks: Matrix of cells' log-likelihoods under each cluster. Cells in rows, clusters in columns. +#' } +#' +quick_celltype <- function(x, bg = 0.01, reference_profiles, nb_size = 10, align_genes = TRUE) { + + if (any(rowSums(x) == 0)) { + stop("Cells with 0 counts were found. Please remove.") + } + + # accept a single value of bg if input by user: + if (length(bg) == 1) { + bg <- rep(bg, nrow(x)) + names(bg) <- rownames(x) + } + + # align genes: + if (align_genes) { + sharedgenes <- intersect(rownames(reference_profiles), colnames(x)) + lostgenes <- setdiff(colnames(x), rownames(reference_profiles)) + + # subset: + x <- x[, sharedgenes] + reference_profiles <- reference_profiles[sharedgenes, ] + + # warn about genes being lost: + if ((length(lostgenes) > 0) && length(lostgenes) < 50) { + message(paste0("The following genes in the count data are missing from reference_profiles and will be omitted from cell typing: ", + paste0(lostgenes, collapse = ","))) + } + if (length(lostgenes) > 50) { + message(paste0(length(lostgenes), " genes in the count data are missing from reference_profiles and will be omitted from cell typing")) + } + } + + + # get logliks + logliks <- parallel::mclapply(asplit(reference_profiles, 2), + lldist, + mat = x, + bg = bg, + size = nb_size, + mc.cores = numCores()) + logliks <- do.call(cbind, logliks) + + + + # get remaining outputs + clust <- colnames(logliks)[apply(logliks, 1, which.max)] + names(clust) <- rownames(logliks) + + out <- list(clust = clust, + logliks = round(logliks, 4)) + return(out) +} diff --git a/R/update_segments.R b/R/update_segments.R index 67c6a8d..4e64e72 100644 --- a/R/update_segments.R +++ b/R/update_segments.R @@ -2,12 +2,14 @@ #' @description Update transcript data.frame based on resegmentation action, calculate the new cell type and mean per cell spatial coordinates #' @param transcript_df the data.frame of transcript to be updated #' @param reseg_full_converter a named converter to update the cell ID in `transcript_df`, cell_ID in name would be converted to cell_ID in value; discard cell_ID with value = NA -#' @param score_GeneMatrix a gene x cell-type score matrix +#' @param score_GeneMatrix the gene x cell-type matrix of log-like score of gene in each cell type, needed if using `LogLikeRatio` cell typing method (default = NULL) +#' @param refProfiles A matrix of cluster profiles, genes X clusters, needed if using `NegBionomial` cell typing method (default = NULL) #' @param transGene_coln the column name of target or gene name in `transcript_df` #' @param cellID_coln the column name of cell_ID in `transcript_df` #' @param celltype_coln the column name of cell type in `transcript_df` #' @param spatLocs_colns column names for 1st, 2nd and optional 3rd dimension of spatial coordinates in transcript_df #' @param return_perCellDF flag to return gene x cell count matrix and per cell DF with updated mean spatial coordinates and new cell type +#' @param celltype_method use either `LogLikeRatio` or `NegBinomial` method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio) #' @return a list #' \describe{ #' \item{updated_transDF}{the updated transcript_df with `updated_cellID` and `updated_celltype` column based on reseg_full_converter} @@ -18,12 +20,25 @@ #' @export update_transDF_ResegActions <- function(transcript_df, reseg_full_converter, - score_GeneMatrix, + score_GeneMatrix = NULL, + refProfiles = NULL, transGene_coln = 'target', cellID_coln = 'cell_ID', celltype_coln = 'cell_type', spatLocs_colns = c("x","y","z"), - return_perCellDF = TRUE){ + return_perCellDF = TRUE, + celltype_method = 'LogLikeRatio'){ + + celltype_method <- match.arg(celltype_method, c('LogLikeRatio', 'NegBinomial')) + + if(celltype_method == 'LogLikeRatio' & is.null(score_GeneMatrix)){ + stop("Must provided `score_GeneMatrix` when using log-likelihood ratio based cell typing method.") + } + + if(celltype_method == 'NegBinomial' & is.null(refProfiles)){ + stop("Must provided `refProfiles` when using negative binomial cell typing method.") + } + # check format of transcript_df if(any(!c(cellID_coln, transGene_coln, celltype_coln, spatLocs_colns) %in% colnames(transcript_df))){ @@ -36,16 +51,25 @@ update_transDF_ResegActions <- function(transcript_df, common_cells <- intersect(unique(names(reseg_full_converter)), unique(transcript_df[[cellID_coln]])) # get common genes - common_genes <- intersect(rownames(score_GeneMatrix), - unique(transcript_df[[transGene_coln]])) - message(sprintf("Found %d common cells and %d common genes among `names(reseg_full_converter)`, `transcript_df`, and `score_GeneMatrix`. ", + if(celltype_method == 'LogLikeRatio'){ + common_genes <- intersect(rownames(score_GeneMatrix), + unique(transcript_df[[transGene_coln]])) + score_GeneMatrix <- score_GeneMatrix[common_genes, ] + } else if (celltype_method =='NegBinomial'){ + common_genes <- intersect(rownames(refProfiles), + unique(transcript_df[[transGene_coln]])) + refProfiles <- refProfiles[common_genes, ] + } + + + message(sprintf("Found %d common cells and %d common genes among `names(reseg_full_converter)`, `transcript_df`, and `score_GeneMatrix` or `refProfiles. ", length(common_cells), length(common_genes))) if(any(length(common_cells) <1, length(common_genes)<1)){ - stop("Too few common cells or genes to proceed. Check if score_GeneMatrix is a gene x cell-type matrix.") + stop("Too few common cells or genes to proceed. Check if `score_GeneMatrix` or `refProfiles` is a gene x cell-type matrix.") } - score_GeneMatrix <- score_GeneMatrix[common_genes, ] + # split reseg_full_converter into different types of cells # get idx @@ -68,19 +92,30 @@ update_transDF_ResegActions <- function(transcript_df, ## get new cell types for cells being updated ---- subTransDF <- transcript_df[which(transcript_df[['updated_cellID']] %in% unique(cells_to_update) & transcript_df[[transGene_coln]] %in% common_genes), ] - # get score for each transcripts - transcriptGeneScore <- score_GeneMatrix[subTransDF[[transGene_coln]], ] - - tmp_score <- as.data.frame(transcriptGeneScore) - tmp_score[['updated_cellID']] <- subTransDF[['updated_cellID']] - - tmp_score <-data.table::setDT(tmp_score)[, lapply(.SD, sum), by = 'updated_cellID'] - tmp_cellID <- tmp_score[['updated_cellID']] - tmp_score[['updated_cellID']] <- NULL - # assign cell type based on max values - max_idx_1st <- max.col(tmp_score,ties.method="first") - newCellTypes <- colnames(tmp_score)[max_idx_1st] - names(newCellTypes) <- tmp_cellID + + if(celltype_method == 'LogLikeRatio'){ + # get score for each transcripts + transcriptGeneScore <- score_GeneMatrix[subTransDF[[transGene_coln]], ] + + tmp_score <- as.data.frame(transcriptGeneScore) + tmp_score[['updated_cellID']] <- subTransDF[['updated_cellID']] + + tmp_score <- data.table::setDT(tmp_score)[, lapply(.SD, sum), by = 'updated_cellID'] + tmp_cellID <- tmp_score[['updated_cellID']] + tmp_score[['updated_cellID']] <- NULL + # assign cell type based on max values + max_idx_1st <- max.col(tmp_score,ties.method="first") + newCellTypes <- colnames(tmp_score)[max_idx_1st] + names(newCellTypes) <- tmp_cellID + + } else if (celltype_method =='NegBinomial'){ + exprMat <- reshape2::acast(subTransDF, as.formula(paste('updated_cellID', '~', transGene_coln)), length) + newCellTypes <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = TRUE)[['clust']] + rm(exprMat) + + } + + # update cell type transcript_df[['updated_celltype']] <- transcript_df[[celltype_coln]] diff --git a/R/wrapper_FastReseg.R b/R/wrapper_FastReseg.R index a994669..80134fa 100644 --- a/R/wrapper_FastReseg.R +++ b/R/wrapper_FastReseg.R @@ -12,7 +12,7 @@ #' @param flagCell_lrtest_cutoff the cutoff of lrtest_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile #' @param svmClass_score_cutoff the cutoff of transcript score to separate between high and low score transcripts in SVM (default = -2) #' @param svm_args a list of arguments to pass to svm function for identifying low-score transcript groups in space, typically involve kernel, gamma, scale -#' @param groupTranscripts_method use either 'delaunay' or 'dbscan' method to group transcripts in space (default = 'delaunay') +#' @param groupTranscripts_method use either `delaunay` or `dbscan` method to group transcripts in space (default = 'delaunay') #' @param cellular_distance_cutoff maximum cell-to-cell distance in x, y between the center of query cells to the center of neighbor cells with direct contact, same unit as input spatial coordinate. Default = NULL to use the 2 times of average 2D cell diameter. #' @param molecular_distance_cutoff maximum molecule-to-molecule distance within connected transcript group, same unit as input spatial coordinate (default = 2.7 micron). #' If set to NULL, the pipeline would first randomly choose no more than 2500 cells from up to 10 random picked ROIs with search radius to be 5 times of `cellular_distance_cutoff`, and then calculate the minimal molecular distance between picked cells. The pipeline would further use the 5 times of 90% quantile of minimal molecular distance as `molecular_distance_cutoff`. This calculation is slow and is not recommended for large transcript data.frame. @@ -25,6 +25,7 @@ #' @param return_perCellData flag to return gene x cell count matrix and per cell DF with updated mean spatial coordinates and new cell type #' @param includeAllRefGenes flag to include all genes in `refProfiles` in the returned `updated_perCellExprs` with missing genes of value 0 (default = FALSE) #' @param ctrl_genes a vector of control genes that are present in input transcript data.frame but not present in `counts` or `refProfiles`; the `ctrl_genes` would be included in FastReseg analysis. (default = NULL) +#' @param celltype_method use either `LogLikeRatio` or `NegBinomial` method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio) #' @return a list #' \describe{ #' \item{modStats_ToFlagCells}{a data.frame for spatial modeling statistics of each cell, output of `score_cell_segmentation_error` function, return when `return_intermediates` = TRUE} @@ -84,9 +85,11 @@ fastReseg_core_externalRef <- function(refProfiles, return_intermediates = TRUE, return_perCellData = TRUE, includeAllRefGenes = FALSE, - ctrl_genes = NULL){ + ctrl_genes = NULL, + celltype_method = 'LogLikeRatio'){ groupTranscripts_method <- match.arg(groupTranscripts_method, c('delaunay', 'dbscan')) + celltype_method <- match.arg(celltype_method, c('LogLikeRatio', 'NegBinomial')) # final results final_res <- list() @@ -174,7 +177,7 @@ fastReseg_core_externalRef <- function(refProfiles, } ## get tLL score matrix - meanCelltype_profiles <- pmax(refProfiles, 1e-5) + meanCelltype_profiles <- pmax(refProfiles, 1e-8) transcript_loglik <- scoreGenesInRef(genes = common_genes, ref_profiles = meanCelltype_profiles) # tLLRv2 score, re-center on maximum per row/transcript @@ -274,17 +277,30 @@ fastReseg_core_externalRef <- function(refProfiles, ## (0.4) for each cell, get new cell type based on maximum score ---- - # `getCellType_maxScore` function returns a list contains element `cellType_DF`, a data.frame with cell in row, cell_ID and cell_type in column. - tmp_df <- getCellType_maxScore(score_GeneMatrix = tLLRv2_geneMatrix, - transcript_df = transcript_df, - transID_coln = transID_coln, - transGene_coln = transGene_coln, - cellID_coln = cellID_coln, - return_transMatrix = FALSE) - - select_cellmeta <- tmp_df[['cellType_DF']] + if(celltype_method == 'LogLikeRatio'){ + # `getCellType_maxScore` function returns a list contains element `cellType_DF`, a data.frame with cell in row, cell_ID and cell_type in column. + tmp_df <- getCellType_maxScore(score_GeneMatrix = tLLRv2_geneMatrix, + transcript_df = transcript_df, + transID_coln = transID_coln, + transGene_coln = transGene_coln, + cellID_coln = cellID_coln, + return_transMatrix = FALSE) + + select_cellmeta <- tmp_df[['cellType_DF']] + rm(tmp_df) + + } else if (celltype_method =='NegBinomial'){ + # `quick_celltype` function returns a list contains element `clust`, a vector given cells' cluster assignments. + exprMat <- reshape2::acast(transcript_df, as.formula(paste0(cellID_coln, "~", transGene_coln)), length) + nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = TRUE) + select_cellmeta <- data.frame(cellID = names(nb_res[['clust']]), + celltype = nb_res[['clust']]) + rm(exprMat, nb_res) + } + + colnames(select_cellmeta) <- c(cellID_coln,'tLLRv2_maxCellType') - rm(tmp_df) + transcript_df <- merge(transcript_df, select_cellmeta, by = cellID_coln) message(sprintf("Found %d cells and assigned cell type based on the provided 'refProfiles` cluster profiles.", nrow(select_cellmeta))) @@ -530,13 +546,25 @@ fastReseg_core_externalRef <- function(refProfiles, flagged_transDF_SVM3[, tmp_cellID := ifelse(connect_group == 0, get(cellID_coln), paste0(get(cellID_coln),'_g', connect_group))] # get new cell type of each group based on maximum - # `getCellType_maxScore` function returns a list contains element `cellType_DF`, a data.frame with cell in row, cell_ID and cell_type in column. - tmp_df <- getCellType_maxScore(score_GeneMatrix = tLLRv2_geneMatrix, - transcript_df = flagged_transDF_SVM3, - transID_coln = transID_coln, - transGene_coln = transGene_coln, - cellID_coln = "tmp_cellID", - return_transMatrix = FALSE) + if(celltype_method == 'LogLikeRatio'){ + # `getCellType_maxScore` function returns a list contains element `cellType_DF`, a data.frame with cell in row, cell_ID and cell_type in column. + tmp_df <- getCellType_maxScore(score_GeneMatrix = tLLRv2_geneMatrix, + transcript_df = flagged_transDF_SVM3, + transID_coln = transID_coln, + transGene_coln = transGene_coln, + cellID_coln = "tmp_cellID", + return_transMatrix = FALSE) + } else if (celltype_method =='NegBinomial'){ + # `quick_celltype` function returns a list contains element `clust`, a vector given cells' cluster assignments. + exprMat <- reshape2::acast(flagged_transDF_SVM3, as.formula(paste0("tmp_cellID", "~", transGene_coln)), length) + nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = TRUE) + tmp_df <- list(cellType_DF = data.frame(cellID = names(nb_res[['clust']]), + celltype = nb_res[['clust']])) + rm(exprMat, nb_res) + } + + + colnames(tmp_df[['cellType_DF']]) <- c('tmp_cellID','group_maxCellType') flagged_transDF_SVM3 <- merge(flagged_transDF_SVM3, tmp_df[['cellType_DF']], by = 'tmp_cellID', all.x = TRUE) flagged_transDF_SVM3 <- as.data.frame(flagged_transDF_SVM3) @@ -571,17 +599,36 @@ fastReseg_core_externalRef <- function(refProfiles, ### search within absolute distance, consider 25um in xy for cell level search and 15 pixel = 2.7um to be direct neighbor at transcript level. # using spatstat to locate neighbor cells and rank them by minimal molecular distance to query cell # `neighborhood_for_resegment_spatstat` function returns a data.frame with each cell in row and its neighborhood information in columns - neighborReSeg_df <- neighborhood_for_resegment_spatstat(chosen_cells = cells_to_use, - score_GeneMatrix = tLLRv2_geneMatrix, - score_baseline = score_baseline, - neighbor_distance_xy = cellular_distance_cutoff, - distance_cutoff = molecular_distance_cutoff, - transcript_df = reseg_transcript_df, - cellID_coln = "tmp_cellID", - celltype_coln = "group_maxCellType", - transID_coln = transID_coln, - transGene_coln = transGene_coln, - transSpatLocs_coln = spatLocs_colns) + if(celltype_method == 'LogLikeRatio'){ + neighborReSeg_df <- neighborhood_for_resegment_spatstat(chosen_cells = cells_to_use, + score_GeneMatrix = tLLRv2_geneMatrix, + score_baseline = score_baseline, + neighbor_distance_xy = cellular_distance_cutoff, + distance_cutoff = molecular_distance_cutoff, + transcript_df = reseg_transcript_df, + cellID_coln = "tmp_cellID", + celltype_coln = "group_maxCellType", + transID_coln = transID_coln, + transGene_coln = transGene_coln, + transSpatLocs_coln = spatLocs_colns, + celltype_method = celltype_method) + } else if (celltype_method =='NegBinomial'){ + neighborReSeg_df <- neighborhood_for_resegment_spatstat(chosen_cells = cells_to_use, + refProfiles = refProfiles, + score_baseline = score_baseline, + neighbor_distance_xy = cellular_distance_cutoff, + distance_cutoff = molecular_distance_cutoff, + transcript_df = reseg_transcript_df, + cellID_coln = "tmp_cellID", + celltype_coln = "group_maxCellType", + transID_coln = transID_coln, + transGene_coln = transGene_coln, + transSpatLocs_coln = spatLocs_colns, + celltype_method = celltype_method) + } + + + #### (4.4) decide resegmentation operation: merge, new cell, or discard ---- # # `decide_ReSegment_Operations_leidenCut` function returns a list containing the following 4 elements: @@ -613,15 +660,28 @@ fastReseg_core_externalRef <- function(refProfiles, # `updated_transDF`, the updated transcript_df with `updated_cellID` and `updated_celltype` column based on reseg_full_converter. # `perCell_DT`, a per cell data.table with mean spatial coordinates and new cell type when return_perCellDF = TRUE. # `perCell_expression`, a gene x cell count sparse matrix for updated transcript data.frame when return_perCellDF = TRUE. + if(celltype_method == 'LogLikeRatio'){ + post_reseg_results <- update_transDF_ResegActions(transcript_df = reseg_transcript_df, + reseg_full_converter = reseg_actions$reseg_full_converter, + score_GeneMatrix = tLLRv2_geneMatrix, + transGene_coln = transGene_coln, + cellID_coln = 'tmp_cellID', + celltype_coln = 'group_maxCellType', + spatLocs_colns = spatLocs_colns, + return_perCellDF = return_perCellData, + celltype_method = celltype_method) + } else if (celltype_method =='NegBinomial'){ + post_reseg_results <- update_transDF_ResegActions(transcript_df = reseg_transcript_df, + reseg_full_converter = reseg_actions$reseg_full_converter, + refProfiles = refProfiles, + transGene_coln = transGene_coln, + cellID_coln = 'tmp_cellID', + celltype_coln = 'group_maxCellType', + spatLocs_colns = spatLocs_colns, + return_perCellDF = return_perCellData, + celltype_method = celltype_method) + } - post_reseg_results <- update_transDF_ResegActions(transcript_df = reseg_transcript_df, - reseg_full_converter = reseg_actions$reseg_full_converter, - score_GeneMatrix = tLLRv2_geneMatrix, - transGene_coln = transGene_coln, - cellID_coln = 'tmp_cellID', - celltype_coln = 'group_maxCellType', - spatLocs_colns = spatLocs_colns, - return_perCellDF = return_perCellData) # get tLLRv2 score under updated cell type @@ -732,6 +792,7 @@ fastReseg_core_externalRef <- function(refProfiles, #' @param return_perCellData flag to return and save to output folder for gene x cell count matrix and per cell DF with updated mean spatial coordinates and new cell type #' @param combine_extra flag to combine original extracellular transcripts and trimmed transcripts back to the updated transcript data.frame, slow process if many transcript in each FOV file. (default = FALSE) #' @param ctrl_genes a vector of control genes that are present in input transcript data.frame but not present in `counts` or `refProfiles`; the `ctrl_genes` would be included in FastReseg analysis. (default = NULL) +#' @param celltype_method use either `LogLikeRatio` or `NegBinomial` method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio) #' @return a list #' \describe{ #' \item{refProfiles}{a genes X clusters matrix of cluster-specific reference profiles used in resegmenation pipeline} @@ -901,10 +962,14 @@ fastReseg_internalRef <- function(counts, save_intermediates = TRUE, return_perCellData = TRUE, combine_extra = FALSE, - ctrl_genes = NULL){ + ctrl_genes = NULL, + celltype_method = 'LogLikeRatio'){ groupTranscripts_method <- match.arg(groupTranscripts_method, c('delaunay', 'dbscan')) - message(sprintf("Use %s for grouping low-score transcripts within each cell in space. ", groupTranscripts_method)) + message(sprintf("Use `%s` for grouping low-score transcripts within each cell in space. ", groupTranscripts_method)) + + celltype_method <- match.arg(celltype_method, c('LogLikeRatio', 'NegBinomial')) + message(sprintf("Use `%s` method for quick cell typing in resegmentation. ", celltype_method)) # spatial dimension d2_or_d3 <- length(spatLocs_colns) @@ -1005,11 +1070,11 @@ fastReseg_internalRef <- function(counts, # lowerCutoff_transNum, a named vector of 25% quantile of cluster-specific per molecule per cell transcript number, to be used as transcript number cutoff such that higher than the cutoff is required to keep query cell as it is # higherCutoff_transNum, a named vector of median value of cluster-specific per molecule per cell transcript number, to be used as transcript number cutoff such that lower than the cutoff is required to keep query cell as it is when there is neighbor cell of consistent cell type. # clust_used, a named vector of cluster assignments for each cell used in baseline calculation, cell_ID in `counts` as name - baselineData <- get_baselineCT(refProfiles = refProfiles, counts = counts, clust = as.character(clust)) + baselineData <- get_baselineCT(refProfiles = refProfiles, counts = counts, clust = as.character(clust), celltype_method = celltype_method) } else { # reference profiles exists, but no cluster assignment - baselineData <- get_baselineCT(refProfiles = refProfiles, counts = counts, clust = NULL) + baselineData <- get_baselineCT(refProfiles = refProfiles, counts = counts, clust = NULL, celltype_method = celltype_method) clust = baselineData[['clust_used']] } @@ -1246,7 +1311,8 @@ fastReseg_internalRef <- function(counts, return_intermediates = save_intermediates, return_perCellData = return_perCellData, includeAllRefGenes = TRUE, - ctrl_genes = ctrl_genes) + ctrl_genes = ctrl_genes, + celltype_method = celltype_method) # intracellular in original and updated segmentation each_segRes[['updated_transDF']][['transComp']] <- 'intraC' @@ -1400,6 +1466,7 @@ fastReseg_internalRef <- function(counts, #' @param path_to_output the file path to output folder; directory would be created by function if not exists; `flagged_transDF`, the reformatted transcript data.frame with transcripts of low goodness-of-fit flagged by` SVM_class = 0`, and `modStats_ToFlagCells`, the per cell evaluation output of segmentation error, and `classDF_ToFlagTrans`, the class assignment of transcripts within each flagged cells are saved as individual csv files for each FOV, respectively. #' @param combine_extra flag to combine original extracellular transcripts back to the flagged transcript data.frame. (default = FALSE) #' @param ctrl_genes a vector of control genes that are present in input transcript data.frame but not present in `counts` or `refProfiles`; the `ctrl_genes` would be included in FastReseg analysis. (default = NULL) +#' @param celltype_method use either `LogLikeRatio` or `NegBinomial` method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio) #' @return a list #' \describe{ #' \item{refProfiles}{a genes * clusters matrix of cluster-specific reference profiles used in resegmenation pipeline} @@ -1502,7 +1569,11 @@ findSegmentError_allFiles <- function(counts, gamma = 0.4), path_to_output = "reSeg_res", combine_extra = FALSE, - ctrl_genes = NULL){ + ctrl_genes = NULL, + celltype_method = 'LogLikeRatio'){ + + celltype_method <- match.arg(celltype_method, c('LogLikeRatio', 'NegBinomial')) + message(sprintf("Use `%s` method for quick cell typing in resegmentation. ", celltype_method)) # spatial dimension d2_or_d3 <- length(spatLocs_colns) @@ -1603,11 +1674,11 @@ findSegmentError_allFiles <- function(counts, # lowerCutoff_transNum, a named vector of 25% quantile of cluster-specific per molecule per cell transcript number, to be used as transcript number cutoff such that higher than the cutoff is required to keep query cell as it is # higherCutoff_transNum, a named vector of median value of cluster-specific per molecule per cell transcript number, to be used as transcript number cutoff such that lower than the cutoff is required to keep query cell as it is when there is neighbor cell of consistent cell type. # clust_used, a named vector of cluster assignments for each cell used in baseline calculation, cell_ID in `counts` as name - baselineData <- get_baselineCT(refProfiles = refProfiles, counts = counts, clust = as.character(clust)) + baselineData <- get_baselineCT(refProfiles = refProfiles, counts = counts, clust = as.character(clust), celltype_method = celltype_method) } else { # reference profiles exists, but no cluster assignment - baselineData <- get_baselineCT(refProfiles = refProfiles, counts = counts, clust = NULL) + baselineData <- get_baselineCT(refProfiles = refProfiles, counts = counts, clust = NULL, celltype_method = celltype_method) clust = baselineData[['clust_used']] } @@ -1633,7 +1704,7 @@ findSegmentError_allFiles <- function(counts, # but also combine perCell data from all FOVs to return ## (0) get transcript score matrix for each gene based on reference profile - transcript_loglik <- scoreGenesInRef(genes = common_genes, ref_profiles = pmax(refProfiles, 1e-5)) + transcript_loglik <- scoreGenesInRef(genes = common_genes, ref_profiles = pmax(refProfiles, 1e-8)) # tLLRv2 score, re-center on maximum per row/transcript tmp_max <- apply(transcript_loglik, 1, max) @@ -1700,18 +1771,30 @@ findSegmentError_allFiles <- function(counts, ## (2) for each cell, get new cell type based on maximum score ---- - # `getCellType_maxScore` function returns a list contains element `cellType_DF`, a data.frame with cell in row, cell_ID and cell_type in column. - tmp_df <- getCellType_maxScore(score_GeneMatrix = tLLRv2_geneMatrix, - transcript_df = transcript_df[['intraC']], - transID_coln = 'UMI_transID', - transGene_coln = 'target', - cellID_coln = 'UMI_cellID', - return_transMatrix = FALSE) + if(celltype_method == 'LogLikeRatio'){ + # `getCellType_maxScore` function returns a list contains element `cellType_DF`, a data.frame with cell in row, cell_ID and cell_type in column. + tmp_df <- getCellType_maxScore(score_GeneMatrix = tLLRv2_geneMatrix, + transcript_df = transcript_df[['intraC']], + transID_coln = 'UMI_transID', + transGene_coln = 'target', + cellID_coln = 'UMI_cellID', + return_transMatrix = FALSE) + + select_cellmeta <- tmp_df[['cellType_DF']] + rm(tmp_df) + + } else if (celltype_method =='NegBinomial'){ + # `quick_celltype` function returns a list contains element `clust`, a vector given cells' cluster assignments. + exprMat <- reshape2::acast(transcript_df[['intraC']], as.formula(paste0("UMI_cellID", "~", "target")), length) + nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = TRUE) + select_cellmeta <- data.frame(cellID = names(nb_res[['clust']]), + celltype = nb_res[['clust']]) + rm(exprMat, nb_res) + + } - select_cellmeta <- tmp_df[['cellType_DF']] - colnames(select_cellmeta) <- c('UMI_cellID','tLLRv2_maxCellType') - rm(tmp_df) + colnames(select_cellmeta) <- c('UMI_cellID','tLLRv2_maxCellType') all_cells <- select_cellmeta[['UMI_cellID']] transcript_df[['intraC']] <- merge(transcript_df[['intraC']], select_cellmeta, by = 'UMI_cellID') @@ -1730,7 +1813,7 @@ findSegmentError_allFiles <- function(counts, - ## (4.1) spatial modeling of tLLR score profile within each cell to identify cells with strong spatail dependency + ## (4.1) spatial modeling of tLLR score profile within each cell to identify cells with strong spatial dependency # `score_cell_segmentation_error` function returns a data.frame with cell in row and spatial modeling outcomes in columns tmp_df <- score_cell_segmentation_error(chosen_cells = all_cells, transcript_df = transcript_df[['intraC']], diff --git a/vignettes/0_flagErrorOnly_on_SMIobject.R b/vignettes/0_flagErrorOnly_on_SMIobject.R index be940e7..effb4e6 100644 --- a/vignettes/0_flagErrorOnly_on_SMIobject.R +++ b/vignettes/0_flagErrorOnly_on_SMIobject.R @@ -38,6 +38,9 @@ removeUnpaired <- FALSE # flag to include ctrl_genes in analysis include_ctrlgenes <- FALSE +# use either `LogLikeRatio` or `NegBinomial` method for quick cell typing and corresponding score_baseline calculation +celltype_method <- 'LogLikeRatio' + # SMI TAP output folder scAnalysis_dir <- "/home/rstudio/NAS_data/lwu/testRun/SpatialTest/giotto_test/melanoma/Run4104_melanoma_980plx/giotto_output/Run4104_cellpose_vs_oldDASH" # get config files to get path to sample annotation file @@ -68,7 +71,8 @@ smi_inputs <- prepSMI_for_fastReseg(path_to_SMIobject = path_to_SMIobject, cellClus_to_exclude = cellClus_to_exclude, removeUnpaired = removeUnpaired, blacklist_genes = blacklist_genes, - pixel_size = pixel_size) + pixel_size = pixel_size, + celltype_method = celltype_method) # write `transDF_fov_fileInfo` into csv file write.csv(smi_inputs[['transDF_fov_fileInfo']], file = fs::path(sub_out_dir, 'transDF_fov_fileInfo.csv')) @@ -130,7 +134,8 @@ reseg_outputs <- findSegmentError_allFiles( extracellular_cellID = c(0), # CellId = 0 means extracelluar transcripts in raw data path_to_output = sub_out_dir, combine_extra = TRUE, # if TRUE, extracellular and trimmed transcripts are included in the updated transcript data.frame - ctrl_genes = ctrl_genes + ctrl_genes = ctrl_genes, + celltype_method = celltype_method ) diff --git a/vignettes/1_fastReseg_on_SMIobject.R b/vignettes/1_fastReseg_on_SMIobject.R index ea78d54..c8617a8 100644 --- a/vignettes/1_fastReseg_on_SMIobject.R +++ b/vignettes/1_fastReseg_on_SMIobject.R @@ -35,6 +35,10 @@ removeUnpaired <- TRUE # flag to include ctrl_genes in analysis include_ctrlgenes <- FALSE +# use either `LogLikeRatio` or `NegBinomial` method for quick cell typing and corresponding score_baseline calculation +celltype_method <- 'LogLikeRatio' + + # SMI TAP output folder scAnalysis_dir <- "/home/rstudio/smiqumulo/01 SMI TAP project/SMI-0003_DavidTing_MGH/6.4 Analysis combined" # get config files to get path to sample annotation file @@ -63,7 +67,8 @@ smi_inputs <- prepSMI_for_fastReseg(path_to_SMIobject = path_to_SMIobject, cellClus_to_exclude = cellClus_to_exclude, removeUnpaired = removeUnpaired, blacklist_genes = blacklist_genes, - pixel_size = pixel_size) + pixel_size = pixel_size, + celltype_method = celltype_method) # write `transDF_fov_fileInfo` into csv file write.csv(smi_inputs[['transDF_fov_fileInfo']], file = fs::path(sub_out_dir, 'transDF_fov_fileInfo.csv')) @@ -150,7 +155,8 @@ reseg_outputs <- fastReseg_internalRef( imputeFlag_missingCTs = TRUE, path_to_output = sub_out_dir, combine_extra = FALSE, # if TRUE, extracellular and trimmed transcripts are included in the updated transcript data.frame - ctrl_genes = ctrl_genes + ctrl_genes = ctrl_genes, + celltype_method = celltype_method ) From 1fb3bda6f2b99e5c06d13753fbb2f691e3f4c70c Mon Sep 17 00:00:00 2001 From: lidanwu Date: Fri, 21 Oct 2022 01:36:12 +0000 Subject: [PATCH 02/10] fastLm f-test for segment error detection; change var name --- R/flag_errors.R | 79 +++++++++++++++++++----- R/wrapper_FastReseg.R | 58 ++++++++--------- vignettes/0_flagErrorOnly_on_SMIobject.R | 14 ++--- vignettes/1_SMI_data_input_cleanup.R | 4 +- vignettes/1_fastReseg_on_SMIobject.R | 4 +- vignettes/2_resegment_and_visualize.R | 72 ++++++++++----------- 6 files changed, 138 insertions(+), 93 deletions(-) diff --git a/R/flag_errors.R b/R/flag_errors.R index 3af90e2..7173e7c 100644 --- a/R/flag_errors.R +++ b/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) } diff --git a/R/wrapper_FastReseg.R b/R/wrapper_FastReseg.R index 80134fa..4f1dcdb 100644 --- a/R/wrapper_FastReseg.R +++ b/R/wrapper_FastReseg.R @@ -9,7 +9,7 @@ #' @param spatLocs_colns column names for 1st, 2nd and optional 3rd dimension of spatial coordinates in `transcript_df` #' @param extracellular_cellID a vector of cell_ID for extracellular transcripts which would be removed from the resegmention pipeline (default = NULL) #' @param flagModel_TransNum_cutoff the cutoff of transcript number to do spatial modeling for identification of wrongly segmented cells (default = 50) -#' @param flagCell_lrtest_cutoff the cutoff of lrtest_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile +#' @param flagCell_lm_cutoff the cutoff of lm_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile #' @param svmClass_score_cutoff the cutoff of transcript score to separate between high and low score transcripts in SVM (default = -2) #' @param svm_args a list of arguments to pass to svm function for identifying low-score transcript groups in space, typically involve kernel, gamma, scale #' @param groupTranscripts_method use either `delaunay` or `dbscan` method to group transcripts in space (default = 'delaunay') @@ -66,7 +66,7 @@ fastReseg_core_externalRef <- function(refProfiles, spatLocs_colns = c('x','y','z'), extracellular_cellID = NULL, flagModel_TransNum_cutoff = 50, - flagCell_lrtest_cutoff = 5, + flagCell_lm_cutoff = 5, svmClass_score_cutoff = -2, svm_args = list(kernel = "radial", scale = FALSE, @@ -299,7 +299,7 @@ fastReseg_core_externalRef <- function(refProfiles, } - colnames(select_cellmeta) <- c(cellID_coln,'tLLRv2_maxCellType') + colnames(select_cellmeta) <- c(cellID_coln,'tSum_maxCellType') transcript_df <- merge(transcript_df, select_cellmeta, by = cellID_coln) @@ -312,7 +312,7 @@ fastReseg_core_externalRef <- function(refProfiles, transcript_df = transcript_df, transID_coln = transID_coln, transGene_coln = transGene_coln, - celltype_coln = 'tLLRv2_maxCellType') + celltype_coln = 'tSum_maxCellType') transcript_df <- merge(transcript_df, tmp_df, by = transID_coln) rm(tmp_df) @@ -323,7 +323,7 @@ fastReseg_core_externalRef <- function(refProfiles, transcript_df = transcript_df, cellID_coln = cellID_coln, transID_coln = transID_coln, - score_coln = 'score_tLLRv2_maxCellType', + score_coln = 'score_tSum_maxCellType', spatLocs_colns = spatLocs_colns, model_cutoff = flagModel_TransNum_cutoff) @@ -338,7 +338,7 @@ fastReseg_core_externalRef <- function(refProfiles, } else{ #-log10(P) - tmp_df[['lrtest_-log10P']] <- -log10(tmp_df[['lrtest_Pr']]) + tmp_df[['lm_-log10P']] <- -log10(tmp_df[['lm_Pvalue']]) modStats_tLLRv2_3D <- merge(select_cellmeta, tmp_df, by.x = cellID_coln, by.y = 'cell_ID') rm(tmp_df) @@ -347,10 +347,10 @@ fastReseg_core_externalRef <- function(refProfiles, } - ## (1.2) flag cells based on linear regression of tLLRv2, lrtest_-log10P - flagged_cells <- modStats_tLLRv2_3D[[cellID_coln]][which(modStats_tLLRv2_3D[['lrtest_-log10P']] > flagCell_lrtest_cutoff )] - message(sprintf("%d cells, %.4f of all evaluated cells, are flagged for resegmentation with lrtest_-log10P > %.1f.", - length(flagged_cells), length(flagged_cells)/nrow(modStats_tLLRv2_3D), flagCell_lrtest_cutoff)) + ## (1.2) flag cells based on linear regression of tLLRv2, lm_-log10P + flagged_cells <- modStats_tLLRv2_3D[[cellID_coln]][which(modStats_tLLRv2_3D[['lm_-log10P']] > flagCell_lm_cutoff )] + message(sprintf("%d cells, %.4f of all evaluated cells, are flagged for resegmentation with lm_-log10P > %.1f.", + length(flagged_cells), length(flagged_cells)/nrow(modStats_tLLRv2_3D), flagCell_lm_cutoff)) } @@ -376,12 +376,12 @@ fastReseg_core_externalRef <- function(refProfiles, reseg_transcript_df <- transcript_df reseg_transcript_df[['connect_group']] <- 0 reseg_transcript_df[['tmp_cellID']] <- reseg_transcript_df[[cellID_coln]] - reseg_transcript_df[['group_maxCellType']] <- reseg_transcript_df[['tLLRv2_maxCellType']] + reseg_transcript_df[['group_maxCellType']] <- reseg_transcript_df[['tSum_maxCellType']] # update transcript df with resegmentation outcomes reseg_transcript_df[['updated_cellID']] <- reseg_transcript_df[[cellID_coln]] - reseg_transcript_df[['updated_celltype']] <- reseg_transcript_df[['tLLRv2_maxCellType']] - reseg_transcript_df[['score_updated_celltype']] <- reseg_transcript_df[['score_tLLRv2_maxCellType']] + reseg_transcript_df[['updated_celltype']] <- reseg_transcript_df[['tSum_maxCellType']] + reseg_transcript_df[['score_updated_celltype']] <- reseg_transcript_df[['score_tSum_maxCellType']] final_res[['updated_transDF']] <- reseg_transcript_df @@ -447,7 +447,7 @@ fastReseg_core_externalRef <- function(refProfiles, transcript_df = flagged_transDF3d, cellID_coln = cellID_coln, transID_coln = transID_coln, - score_coln = 'score_tLLRv2_maxCellType', + score_coln = 'score_tSum_maxCellType', spatLocs_colns = spatLocs_colns, model_cutoff = flagModel_TransNum_cutoff, score_cutoff = svmClass_score_cutoff, @@ -587,7 +587,7 @@ fastReseg_core_externalRef <- function(refProfiles, tmp_idx <- which(is.na(reseg_transcript_df[['connect_group']])) reseg_transcript_df[['connect_group']][tmp_idx]<-rep(0, length(tmp_idx)) reseg_transcript_df[['tmp_cellID']][tmp_idx] <- reseg_transcript_df[[cellID_coln]][tmp_idx] - reseg_transcript_df[['group_maxCellType']][tmp_idx] <- reseg_transcript_df[['tLLRv2_maxCellType']][tmp_idx] + reseg_transcript_df[['group_maxCellType']][tmp_idx] <- reseg_transcript_df[['tSum_maxCellType']][tmp_idx] rm(tmp_idx) ## (4.3) evaluate the neighborhood of each group for re-segmentation ---- @@ -774,7 +774,7 @@ fastReseg_core_externalRef <- function(refProfiles, #' @param spatLocs_colns column names for 1st, 2nd and optional 3rd dimension of spatial coordinates in `transcript_df` #' @param extracellular_cellID a vector of cell_ID for extracellular transcripts which would be removed from the resegmention pipeline (default = NULL) #' @param flagModel_TransNum_cutoff the cutoff of transcript number to do spatial modeling for identification of wrongly segmented cells (default = 50) -#' @param flagCell_lrtest_cutoff the cutoff of lrtest_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile +#' @param flagCell_lm_cutoff the cutoff of lm_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile #' @param svmClass_score_cutoff the cutoff of transcript score to separate between high and low score transcripts in SVM (default = -2) #' @param svm_args a list of arguments to pass to svm function for identifying low-score transcript groups in space, typically involve kernel, gamma, scale #' @param groupTranscripts_method use either 'delaunay' or 'dbscan' method to group transcripts in space (default = 'delaunay') @@ -941,7 +941,7 @@ fastReseg_internalRef <- function(counts, spatLocs_colns = c('x','y','z'), extracellular_cellID = NULL, flagModel_TransNum_cutoff = 50, - flagCell_lrtest_cutoff = 5, + flagCell_lm_cutoff = 5, svmClass_score_cutoff = -2, svm_args = list(kernel = "radial", scale = FALSE, @@ -1302,7 +1302,7 @@ fastReseg_internalRef <- function(counts, cellID_coln = 'UMI_cellID', spatLocs_colns = c('x','y','z')[1:d2_or_d3], flagModel_TransNum_cutoff = flagModel_TransNum_cutoff, - flagCell_lrtest_cutoff = flagCell_lrtest_cutoff, + flagCell_lm_cutoff = flagCell_lm_cutoff, svmClass_score_cutoff = svmClass_score_cutoff, svm_args = svm_args, groupTranscripts_method = groupTranscripts_method, @@ -1460,7 +1460,7 @@ fastReseg_internalRef <- function(counts, #' @param spatLocs_colns column names for 1st, 2nd and optional 3rd dimension of spatial coordinates in `transcript_df` #' @param extracellular_cellID a vector of cell_ID for extracellular transcripts which would be removed from the resegmention pipeline (default = NULL) #' @param flagModel_TransNum_cutoff the cutoff of transcript number to do spatial modeling for identification of wrongly segmented cells (default = 50) -#' @param flagCell_lrtest_cutoff the cutoff of lrtest_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile +#' @param flagCell_lm_cutoff the cutoff of lm_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile #' @param svmClass_score_cutoff the cutoff of transcript score to separate between high and low score transcripts in SVM (default = -2) #' @param svm_args a list of arguments to pass to svm function for identifying low-score transcript groups in space, typically involve kernel, gamma, scale #' @param path_to_output the file path to output folder; directory would be created by function if not exists; `flagged_transDF`, the reformatted transcript data.frame with transcripts of low goodness-of-fit flagged by` SVM_class = 0`, and `modStats_ToFlagCells`, the per cell evaluation output of segmentation error, and `classDF_ToFlagTrans`, the class assignment of transcripts within each flagged cells are saved as individual csv files for each FOV, respectively. @@ -1562,7 +1562,7 @@ findSegmentError_allFiles <- function(counts, spatLocs_colns = c('x','y','z'), extracellular_cellID = NULL, flagModel_TransNum_cutoff = 50, - flagCell_lrtest_cutoff = 5, + flagCell_lm_cutoff = 5, svmClass_score_cutoff = -2, svm_args = list(kernel = "radial", scale = FALSE, @@ -1794,7 +1794,7 @@ findSegmentError_allFiles <- function(counts, } - colnames(select_cellmeta) <- c('UMI_cellID','tLLRv2_maxCellType') + colnames(select_cellmeta) <- c('UMI_cellID','tSum_maxCellType') all_cells <- select_cellmeta[['UMI_cellID']] transcript_df[['intraC']] <- merge(transcript_df[['intraC']], select_cellmeta, by = 'UMI_cellID') @@ -1807,7 +1807,7 @@ findSegmentError_allFiles <- function(counts, transcript_df = transcript_df[['intraC']], transID_coln = 'UMI_transID', transGene_coln = 'target', - celltype_coln = 'tLLRv2_maxCellType') + celltype_coln = 'tSum_maxCellType') transcript_df[['intraC']] <- merge(transcript_df[['intraC']], tmp_df, by = 'UMI_transID') rm(tmp_df) @@ -1819,7 +1819,7 @@ findSegmentError_allFiles <- function(counts, transcript_df = transcript_df[['intraC']], cellID_coln = 'UMI_cellID', transID_coln = 'UMI_transID', - score_coln = 'score_tLLRv2_maxCellType', + score_coln = 'score_tSum_maxCellType', spatLocs_colns = c('x','y','z')[1:d2_or_d3], model_cutoff = flagModel_TransNum_cutoff) @@ -1830,16 +1830,16 @@ findSegmentError_allFiles <- function(counts, } else { #-log10(P) - tmp_df[['lrtest_-log10P']] <- -log10(tmp_df[['lrtest_Pr']]) + tmp_df[['lm_-log10P']] <- -log10(tmp_df[['lm_Pvalue']]) modStats_ToFlagCells <- merge(select_cellmeta, tmp_df, by.x = 'UMI_cellID', by.y = 'cell_ID') rm(tmp_df) - ## (4.2) flag cells based on linear regression of tLLRv2, lrtest_-log10P - modStats_ToFlagCells[['flagged']] <- (modStats_ToFlagCells[['lrtest_-log10P']] > flagCell_lrtest_cutoff ) + ## (4.2) flag cells based on linear regression of tLLRv2, lm_-log10P + modStats_ToFlagCells[['flagged']] <- (modStats_ToFlagCells[['lm_-log10P']] > flagCell_lm_cutoff ) flagged_cells <- modStats_ToFlagCells[['UMI_cellID']][modStats_ToFlagCells[['flagged']]] - message(sprintf("%d cells, %.4f of all evaluated cells, are flagged for resegmentation with lrtest_-log10P > %.1f.", - length(flagged_cells), length(flagged_cells)/nrow(modStats_ToFlagCells), flagCell_lrtest_cutoff)) + message(sprintf("%d cells, %.4f of all evaluated cells, are flagged for resegmentation with lm_-log10P > %.1f.", + length(flagged_cells), length(flagged_cells)/nrow(modStats_ToFlagCells), flagCell_lm_cutoff)) # write into disk # add idx as file idx @@ -1864,7 +1864,7 @@ findSegmentError_allFiles <- function(counts, transcript_df = classDF_ToFlagTrans, cellID_coln = 'UMI_cellID', transID_coln = 'UMI_transID', - score_coln = 'score_tLLRv2_maxCellType', + score_coln = 'score_tSum_maxCellType', spatLocs_colns = c('x','y','z')[1:d2_or_d3], model_cutoff = flagModel_TransNum_cutoff, score_cutoff = svmClass_score_cutoff, diff --git a/vignettes/0_flagErrorOnly_on_SMIobject.R b/vignettes/0_flagErrorOnly_on_SMIobject.R index effb4e6..eec76d3 100644 --- a/vignettes/0_flagErrorOnly_on_SMIobject.R +++ b/vignettes/0_flagErrorOnly_on_SMIobject.R @@ -82,8 +82,8 @@ write.csv(smi_inputs[['transDF_fov_fileInfo']], # # cutoff of transcript number to do spatial modeling for identification of wrongly segmented cells (default = 50) # flagModel_TransNum_cutoff = 50 # -# # cutoff of lrtest_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile -# flagCell_lrtest_cutoff = 5 +# # cutoff of lm_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile +# flagCell_lm_cutoff = 5 # # # cutoff of transcript score to separate between high and low score transcripts in SVM (default = -2) # svmClass_score_cutoff = -2 @@ -235,16 +235,16 @@ message(sprintf("%d cells, %.2f%% of all cells, are flagged for potential cell s combined_modStats_ToFlagCells <- reseg_outputs$combined_modStats_ToFlagCells -# cutoff of lrtest_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile (default =5) +# cutoff of lm_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile (default =5) # lower values would flag more cells with potential segmentation error -flagCell_lrtest_cutoff = 3 +flagCell_lm_cutoff = 3 # cells with potential segmentation errors, flagged by new cutoff -combined_flaggedCells <- combined_modStats_ToFlagCells[combined_modStats_ToFlagCells['lrtest_-log10P'] > flagCell_lrtest_cutoff, 'UMI_cellID'] -message(sprintf("%d cells, %.2f%% of all cells, are flagged for potential cell segemntation error based on provided `flagCell_lrtest_cutoff` = %.2f. ", +combined_flaggedCells <- combined_modStats_ToFlagCells[combined_modStats_ToFlagCells['lm_-log10P'] > flagCell_lm_cutoff, 'UMI_cellID'] +message(sprintf("%d cells, %.2f%% of all cells, are flagged for potential cell segemntation error based on provided `flagCell_lm_cutoff` = %.2f. ", length(combined_flaggedCells), length(combined_flaggedCells)/nrow(reseg_outputs$combined_modStats_ToFlagCells), - flagCell_lrtest_cutoff)) + flagCell_lm_cutoff)) #### (6) redo identification of low goodness-of-fit transcript groups as needed ---- diff --git a/vignettes/1_SMI_data_input_cleanup.R b/vignettes/1_SMI_data_input_cleanup.R index d4beb08..2ee1724 100644 --- a/vignettes/1_SMI_data_input_cleanup.R +++ b/vignettes/1_SMI_data_input_cleanup.R @@ -41,8 +41,8 @@ config_dimension[['CellNeighbor_z']] = config_dimension[['zstep_um']]*4 config_dimension[['CellNeighbor_xy_in_transDF']] = config_dimension[['CellNeighbor_xy']]/3 -# lrtest_-log10P cutoff for flagging wronlgy segmented cells -config_dimension[['flagCell_lrtestCutoff']] = 10 +# lm_-log10P cutoff for flagging wronlgy segmented cells +config_dimension[['flagCell_lmCutoff']] = 10 # high and low tLLRv2 score cutoff for SVM config_dimension[['tLLRv2_SVMcutoff']] = -2 diff --git a/vignettes/1_fastReseg_on_SMIobject.R b/vignettes/1_fastReseg_on_SMIobject.R index c8617a8..c818b8b 100644 --- a/vignettes/1_fastReseg_on_SMIobject.R +++ b/vignettes/1_fastReseg_on_SMIobject.R @@ -82,8 +82,8 @@ save(smi_inputs, file = fs::path(sub_out_dir, "smi_inputs.RData")) # # cutoff of transcript number to do spatial modeling for identification of wrongly segmented cells (default = 50) # flagModel_TransNum_cutoff = 50 # -# # cutoff of lrtest_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile -# flagCell_lrtest_cutoff = 5 +# # cutoff of lm_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile +# flagCell_lm_cutoff = 5 # # # cutoff of transcript score to separate between high and low score transcripts in SVM (default = -2) # svmClass_score_cutoff = -2 diff --git a/vignettes/2_resegment_and_visualize.R b/vignettes/2_resegment_and_visualize.R index ec52c0d..d36cf98 100644 --- a/vignettes/2_resegment_and_visualize.R +++ b/vignettes/2_resegment_and_visualize.R @@ -1,7 +1,7 @@ #### pipeline to process multi-FOV multi-slot data ## part 2: ## (2.0) load exisitng RData containing data.frame for all transcripts and cells and the reference cell type profiles -## (2.1) flag cells based on linear regression of transcript score using lrtest_-log10P cutoff +## (2.1) flag cells based on linear regression of transcript score using lm_-log10P cutoff ## (2.2) use SVM~hyperplane to identify the connected transcripts group based on transcript score ## (2.3) do network analysis on flagged transcript to split flagged transcript groups in space ## (2.4) re-segmentation of each flagged transcript groups based on their neighborhood @@ -24,7 +24,7 @@ config_dimension[['CellNeighbor_xy']] = 25 config_dimension[['CellNeighbor_xy_in_transDF']] = 25 # change flaging cutoff to 5 -config_dimension[['flagCell_lrtestCutoff']] = 5 +config_dimension[['flagCell_lmCutoff']] = 5 # change svm configuraiton config_dimension[['svm_config']] <- list(kernel = "radial", @@ -96,18 +96,18 @@ system.time(tmp_df <- getCellType_maxScore(score_GeneMatrix = tLLRv2_geneMatrix_ # user system elapsed # 431.748 61.822 96.555 -colnames(tmp_df[['cellType_DF']]) <- c('cell_ID','cleaned_tLLRv2_maxCellType') +colnames(tmp_df[['cellType_DF']]) <- c('cell_ID','cleaned_tSum_maxCellType') ## cells with no prior cell typing info in gem would be removed, due to failed QC select_cellmeta <- merge(select_cellmeta, tmp_df[['cellType_DF']], by = 'cell_ID') # get score on assigned cell types at transcript level ## option 1: not used -#### keep all cells in all_transDF, update the cell_type of missing cells based on cleaned_tLLRv2_maxCellType +#### keep all cells in all_transDF, update the cell_type of missing cells based on cleaned_tSum_maxCellType #### remove c_0 from data.frame ## option 2: only keep the transcripts have cell_type information in gem all_transDF <- merge(all_transDF, - select_cellmeta[, c('cell_ID','cell_type','cleaned_tLLRv2_maxCellType')], + select_cellmeta[, c('cell_ID','cell_type','cleaned_tSum_maxCellType')], by = 'cell_ID') # need to have NotDet in reference if calculate tLLRv2 score for orignla cell type @@ -117,7 +117,7 @@ tmp_df <- getScoreCellType_gene(score_GeneMatrix = tLLRv2_geneMatrix_cleaned, transcript_df = all_transDF, transID_coln = "transcript_id", transGene_coln = "target", - celltype_coln = 'cleaned_tLLRv2_maxCellType') + celltype_coln = 'cleaned_tSum_maxCellType') all_transDF <- merge(all_transDF, tmp_df, by = 'transcript_id') reseg_logInfo[['common_data']][['feat_count']][['cellNum']] <- nrow(select_cellmeta) @@ -144,26 +144,26 @@ if(TRUE){ transcript_df = all_transDF, cellID_coln = 'cell_ID', transID_coln = 'transcript_id', - score_coln = 'score_cleaned_tLLRv2_maxCellType', + score_coln = 'score_cleaned_tSum_maxCellType', spatLocs_colns = c('x','y','z'), model_cutoff = 50)) #-log10(P) - tmp_df[['lrtest_-log10P']] <- -log10(tmp_df[['lrtest_Pr']]) - modStats_cleaned_tLLRv2_3D <- merge(tmp_df, select_cellmeta[, c('cell_ID','cleaned_tLLRv2_maxCellType','slide')], by = 'cell_ID') + tmp_df[['lm_-log10P']] <- -log10(tmp_df[['lm_Pvalue']]) + modStats_cleaned_tLLRv2_3D <- merge(tmp_df, select_cellmeta[, c('cell_ID','cleaned_tSum_maxCellType','slide')], by = 'cell_ID') # visualize extreme cells, 2D plots if(TRUE){ tmp_df <- data.table::copy(modStats_cleaned_tLLRv2_3D) - tmp_df <- tmp_df[order(tmp_df[['lrtest_-log10P']]),] - tmp_df[['labels']] <- paste0(tmp_df[['slide']],'_', tmp_df[['cleaned_tLLRv2_maxCellType']], - ', -log10P=', round(tmp_df[['lrtest_-log10P']],2)) + tmp_df <- tmp_df[order(tmp_df[['lm_-log10P']]),] + tmp_df[['labels']] <- paste0(tmp_df[['slide']],'_', tmp_df[['cleaned_tSum_maxCellType']], + ', -log10P=', round(tmp_df[['lm_-log10P']],2)) chosen_cells <- c(tmp_df$cell_ID[1:9], tmp_df$cell_ID[(nrow(tmp_df)-9):nrow(tmp_df)]) fig1 <- plotSpatialScoreMultiCells(chosen_cells = chosen_cells, cell_labels = tmp_df[match(chosen_cells, tmp_df$cell_ID), 'labels'], transcript_df = all_transDF, cellID_coln = "cell_ID", transID_coln = "transcript_id", - score_coln = "score_cleaned_tLLRv2_maxCellType", + score_coln = "score_cleaned_tSum_maxCellType", spatLocs_colns = c("x","y")) pdf(fs::path(sub_out_dir3, paste0(blockID,"_SpatialPlot_SpatModel2_tLLRv2_-log10P_cleanNBclust6_flagExtreme.pdf")), @@ -172,13 +172,13 @@ if(TRUE){ dev.off() # histogram for linear regression -log10P values - fig <- ggplot(tmp_df, aes(x = get('lrtest_-log10P'))) + + fig <- ggplot(tmp_df, aes(x = get('lm_-log10P'))) + geom_histogram(aes(y=..density..), fill = 'blue',color = 'black')+ - geom_vline(xintercept= quantile(tmp_df[['lrtest_-log10P']], 0.9), + geom_vline(xintercept= quantile(tmp_df[['lm_-log10P']], 0.9), linetype="dashed", color = 'red')+ - labs(y = 'density', x = 'lrtest_-log10P', + labs(y = 'density', x = 'lm_-log10P', title = paste0(nrow(tmp_df),' cells above model_cutoff, skip ', length(common_cells) - nrow(tmp_df), ' cells')) - ggsave(plot = fig, filename = fs::path(sub_out_dir3, paste0(blockID,"_histogram_lrtest_-log10P_allCells.jpeg"))) + ggsave(plot = fig, filename = fs::path(sub_out_dir3, paste0(blockID,"_histogram_lm_-log10P_allCells.jpeg"))) @@ -186,11 +186,11 @@ if(TRUE){ } ####re-segmentation based on tLLRv2 score (1) flag cells, identify transcript groups ---- -## (1) flag cells based on linear regression of tLLRv2, lrtest_-log10P -#5640 cells, 0.0292 of all evaluated cells, are flagged for resegmentation with lrtest_-log10P > 5.0. -flagged_cells_cleaned <- modStats_cleaned_tLLRv2_3D[['cell_ID']][which(modStats_cleaned_tLLRv2_3D[['lrtest_-log10P']] > config_dimension[['flagCell_lrtestCutoff']])] -message(sprintf("%d cells, %.4f of all evaluated cells, are flagged for resegmentation with lrtest_-log10P > %.1f.", - length(flagged_cells_cleaned), length(flagged_cells_cleaned)/nrow(modStats_cleaned_tLLRv2_3D),config_dimension[['flagCell_lrtestCutoff']])) +## (1) flag cells based on linear regression of tLLRv2, lm_-log10P +#5640 cells, 0.0292 of all evaluated cells, are flagged for resegmentation with lm_-log10P > 5.0. +flagged_cells_cleaned <- modStats_cleaned_tLLRv2_3D[['cell_ID']][which(modStats_cleaned_tLLRv2_3D[['lm_-log10P']] > config_dimension[['flagCell_lmCutoff']])] +message(sprintf("%d cells, %.4f of all evaluated cells, are flagged for resegmentation with lm_-log10P > %.1f.", + length(flagged_cells_cleaned), length(flagged_cells_cleaned)/nrow(modStats_cleaned_tLLRv2_3D),config_dimension[['flagCell_lmCutoff']])) ## (2) use SVM~hyperplane to identify the connected transcripts group based on tLLRv2 score ---- @@ -210,7 +210,7 @@ system.time(tmp_df <- flagTranscripts_SVM(chosen_cells = flagged_cells_cleaned, transcript_df = flagged_transDF3d_cleaned, cellID_coln = 'cell_ID', transID_coln = 'transcript_id', - score_coln = 'score_cleaned_tLLRv2_maxCellType', + score_coln = 'score_cleaned_tSum_maxCellType', spatLocs_colns = c('x','y','z'), model_cutoff = 50, score_cutoff = flag_tLLRv2_cutoff, @@ -237,14 +237,14 @@ reseg_logInfo[['flagging_cleaned_SVM']][['flagged_transID']] <- flaggedSVM_trans # visualize the spatial plot of score # choose 500 cells to plot, space across -log10P values tmp_df <- modStats_cleaned_tLLRv2_3D[modStats_cleaned_tLLRv2_3D[['cell_ID']] %in% flagged_cells_cleaned, ] -tmp_df <- tmp_df[order(tmp_df[['lrtest_-log10P']]),] +tmp_df <- tmp_df[order(tmp_df[['lm_-log10P']]),] cells_for_plots <- tmp_df[['cell_ID']][seq(1, length(flagged_cells_cleaned), by = round(length(flagged_cells_cleaned)/500))] cells_for_plots <- unique(cells_for_plots) if(TRUE){ tmp_df <- flagged_transDF_SVM3 - for(score_coln in c('score_cleaned_tLLRv2_maxCellType','DecVal')){ + for(score_coln in c('score_cleaned_tSum_maxCellType','DecVal')){ # visualize the flagged transcripts in chosen_cells if(score_coln == 'DecVal') { score_mid = 0 @@ -798,7 +798,7 @@ if(TRUE){ model_cutoff = 50) #-log10(P) - tmp_df[['lrtest_-log10P']] <- -log10(tmp_df[['lrtest_Pr']]) + tmp_df[['lm_-log10P']] <- -log10(tmp_df[['lm_Pr']]) alteredOnly_modStats_resegcleanSVM_leiden_tLLRv2_3D <- merge(tmp_df, post_reseg_results_cleanSVM_leiden$perCell_DT[, .SD, .SDcols = c('updated_cellID','updated_celltype')], by.x = 'cell_ID', by.y = 'updated_cellID') @@ -809,9 +809,9 @@ if(TRUE){ if(TRUE){ tmp_df <- data.table::copy(alteredOnly_modStats_resegcleanSVM_leiden_tLLRv2_3D) colnames(tmp_df)[1] <- 'cell_ID' - tmp_df <- tmp_df[order(tmp_df[['lrtest_-log10P']]),] + tmp_df <- tmp_df[order(tmp_df[['lm_-log10P']]),] tmp_df[['labels']] <- paste0(tmp_df[['slide']],'_', tmp_df[['updated_celltype']], - ', -log10P=', round(tmp_df[['lrtest_-log10P']],2)) + ', -log10P=', round(tmp_df[['lm_-log10P']],2)) chosen_cells <- c(tmp_df$cell_ID[1:9], tmp_df$cell_ID[(nrow(tmp_df)-9):nrow(tmp_df)]) fig1 <- plotSpatialScoreMultiCells(chosen_cells = chosen_cells, cell_labels = tmp_df[match(chosen_cells, tmp_df$cell_ID), 'labels'], @@ -827,13 +827,13 @@ if(TRUE){ dev.off() # histogram for linear regression -log10P values - fig <- ggplot(tmp_df, aes(x = get('lrtest_-log10P'))) + + fig <- ggplot(tmp_df, aes(x = get('lm_-log10P'))) + geom_histogram(aes(y=..density..), fill = 'blue',color = 'black')+ - geom_vline(xintercept= quantile(tmp_df[['lrtest_-log10P']], 0.9), + geom_vline(xintercept= quantile(tmp_df[['lm_-log10P']], 0.9), linetype="dashed", color = 'red')+ - labs(y = 'density', x = 'lrtest_-log10P', + labs(y = 'density', x = 'lm_-log10P', title = paste0(nrow(tmp_df),' cells above model_cutoff, skip ', length(tmp_cellID) - nrow(tmp_df), ' cells')) - ggsave(plot = fig, filename = fs::path(sub_out_dir3, paste0(blockID,"_histogram_lrtest_-log10P_resegmented_AlteredOnly_cleanNBclust6_SVMleiden.jpeg"))) + ggsave(plot = fig, filename = fs::path(sub_out_dir3, paste0(blockID,"_histogram_lm_-log10P_resegmented_AlteredOnly_cleanNBclust6_SVMleiden.jpeg"))) } } @@ -866,15 +866,15 @@ if(TRUE){ colnames(tmp2_df) <- c('slide', 'flaggedCells') stats_df <- merge(stats_df, tmp2_df, by = 'slide') - fig <- ggplot(tmp_df, aes(x = get('lrtest_-log10P'), group = as.factor(slide), + fig <- ggplot(tmp_df, aes(x = get('lm_-log10P'), group = as.factor(slide), color = as.factor(slide), fill = as.factor(slide))) + geom_histogram(position="dodge2", alpha = 0.5)+ geom_vline(xintercept= 10, linetype="dashed", color = 'black')+ - labs(x = 'lrtest_-log10P', color = "slide", fill = "slide", + labs(x = 'lm_-log10P', color = "slide", fill = "slide", title = paste0("slide ", paste0(tmp2_df$slide, collapse = ":"), "w/ flagged cell #, ", paste0(tmp2_df$flaggedCells, collapse = ":")))+ theme_linedraw() - ggsave(plot = fig, filename = fs::path(sub_out_dir3, paste0(blockID,"_histogram_lrtest_-log10P_slide-comparison_cleanNBclust6_SVMleiden.jpeg"))) + ggsave(plot = fig, filename = fs::path(sub_out_dir3, paste0(blockID,"_histogram_lm_-log10P_slide-comparison_cleanNBclust6_SVMleiden.jpeg"))) @@ -944,7 +944,7 @@ if(TRUE){ tmp_df[['slide']] <- sapply(strsplit(tmp_df$cell_ID, "_"),"[[", 2) tmp_df[['fov']] <- sapply(strsplit(tmp_df$cell_ID, "_"),"[[", 3) tmp_df <- tmp_df[which(tmp_df$fov ==1), ] - tmp_df <- tmp_df[order(tmp_df[['lrtest_-log10P']]),] + tmp_df <- tmp_df[order(tmp_df[['lm_-log10P']]),] cells_for_plots <- tmp_df[['cell_ID']] if(length(cells_for_plots) >510){ From e61e4cb79339b02441351e9f74eddd6fa5b12412 Mon Sep 17 00:00:00 2001 From: lidanwu Date: Fri, 21 Oct 2022 01:51:04 +0000 Subject: [PATCH 03/10] update mans --- NAMESPACE | 5 ++++ R/quick_celltyping.R | 3 +- man/fastReseg_core_externalRef.Rd | 11 ++++--- man/fastReseg_internalRef.Rd | 9 ++++-- man/findSegmentError_allFiles.Rd | 9 ++++-- man/get_baselineCT.Rd | 9 +++++- man/lldist.Rd | 25 ++++++++++++++++ man/neighborhood_for_resegment_spatstat.Rd | 14 ++++++--- man/numCores.Rd | 14 +++++++++ man/prepSMI_for_fastReseg.Rd | 5 +++- man/quick_celltype.Rd | 35 ++++++++++++++++++++++ man/score_cell_segmentation_error.Rd | 8 ++--- man/update_transDF_ResegActions.Rd | 12 ++++++-- 13 files changed, 135 insertions(+), 24 deletions(-) create mode 100644 man/lldist.Rd create mode 100644 man/numCores.Rd create mode 100644 man/quick_celltype.Rd diff --git a/NAMESPACE b/NAMESPACE index 7a1d160..9ea1adc 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -19,10 +19,13 @@ 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) @@ -30,6 +33,7 @@ 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) diff --git a/R/quick_celltyping.R b/R/quick_celltyping.R index 8fee5a1..af290ec 100644 --- a/R/quick_celltyping.R +++ b/R/quick_celltyping.R @@ -9,6 +9,7 @@ #' @importFrom Matrix rowSums #' @importFrom stats dnbinom #' @return a cell x cell_type matrix of the log-likelihood +#' @export lldist <- function(x, mat, bg = 0.01, size = 10, digits = 2) { # convert to matrix form if only a vector was input: if (is.vector(mat)) { @@ -93,7 +94,7 @@ numCores <- function() { #' \item clust: a vector given cells' cluster assignments #' \item logliks: Matrix of cells' log-likelihoods under each cluster. Cells in rows, clusters in columns. #' } -#' +#' @export quick_celltype <- function(x, bg = 0.01, reference_profiles, nb_size = 10, align_genes = TRUE) { if (any(rowSums(x) == 0)) { diff --git a/man/fastReseg_core_externalRef.Rd b/man/fastReseg_core_externalRef.Rd index c872333..28c05ab 100644 --- a/man/fastReseg_core_externalRef.Rd +++ b/man/fastReseg_core_externalRef.Rd @@ -13,7 +13,7 @@ fastReseg_core_externalRef( spatLocs_colns = c("x", "y", "z"), extracellular_cellID = NULL, flagModel_TransNum_cutoff = 50, - flagCell_lrtest_cutoff = 5, + flagCell_lm_cutoff = 5, svmClass_score_cutoff = -2, svm_args = list(kernel = "radial", scale = FALSE, gamma = 0.4), groupTranscripts_method = c("delaunay", "dbscan"), @@ -28,7 +28,8 @@ fastReseg_core_externalRef( return_intermediates = TRUE, return_perCellData = TRUE, includeAllRefGenes = FALSE, - ctrl_genes = NULL + ctrl_genes = NULL, + celltype_method = "LogLikeRatio" ) } \arguments{ @@ -48,13 +49,13 @@ fastReseg_core_externalRef( \item{flagModel_TransNum_cutoff}{the cutoff of transcript number to do spatial modeling for identification of wrongly segmented cells (default = 50)} -\item{flagCell_lrtest_cutoff}{the cutoff of lrtest_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile} +\item{flagCell_lm_cutoff}{the cutoff of lm_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile} \item{svmClass_score_cutoff}{the cutoff of transcript score to separate between high and low score transcripts in SVM (default = -2)} \item{svm_args}{a list of arguments to pass to svm function for identifying low-score transcript groups in space, typically involve kernel, gamma, scale} -\item{groupTranscripts_method}{use either 'delaunay' or 'dbscan' method to group transcripts in space (default = 'delaunay')} +\item{groupTranscripts_method}{use either \code{delaunay} or \code{dbscan} method to group transcripts in space (default = 'delaunay')} \item{molecular_distance_cutoff}{maximum molecule-to-molecule distance within connected transcript group, same unit as input spatial coordinate (default = 2.7 micron). If set to NULL, the pipeline would first randomly choose no more than 2500 cells from up to 10 random picked ROIs with search radius to be 5 times of \code{cellular_distance_cutoff}, and then calculate the minimal molecular distance between picked cells. The pipeline would further use the 5 times of 90\% quantile of minimal molecular distance as \code{molecular_distance_cutoff}. This calculation is slow and is not recommended for large transcript data.frame.} @@ -78,6 +79,8 @@ If set to NULL, the pipeline would first randomly choose no more than 2500 cells \item{includeAllRefGenes}{flag to include all genes in \code{refProfiles} in the returned \code{updated_perCellExprs} with missing genes of value 0 (default = FALSE)} \item{ctrl_genes}{a vector of control genes that are present in input transcript data.frame but not present in \code{counts} or \code{refProfiles}; the \code{ctrl_genes} would be included in FastReseg analysis. (default = NULL)} + +\item{celltype_method}{use either \code{LogLikeRatio} or \code{NegBinomial} method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio)} } \value{ a list diff --git a/man/fastReseg_internalRef.Rd b/man/fastReseg_internalRef.Rd index 88b0f00..23a1a58 100644 --- a/man/fastReseg_internalRef.Rd +++ b/man/fastReseg_internalRef.Rd @@ -21,7 +21,7 @@ fastReseg_internalRef( spatLocs_colns = c("x", "y", "z"), extracellular_cellID = NULL, flagModel_TransNum_cutoff = 50, - flagCell_lrtest_cutoff = 5, + flagCell_lm_cutoff = 5, svmClass_score_cutoff = -2, svm_args = list(kernel = "radial", scale = FALSE, gamma = 0.4), groupTranscripts_method = c("delaunay", "dbscan"), @@ -38,7 +38,8 @@ fastReseg_internalRef( save_intermediates = TRUE, return_perCellData = TRUE, combine_extra = FALSE, - ctrl_genes = NULL + ctrl_genes = NULL, + celltype_method = "LogLikeRatio" ) } \arguments{ @@ -73,7 +74,7 @@ Notice that some assays like SMI has XY axes swapped between stage and each FOV \item{flagModel_TransNum_cutoff}{the cutoff of transcript number to do spatial modeling for identification of wrongly segmented cells (default = 50)} -\item{flagCell_lrtest_cutoff}{the cutoff of lrtest_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile} +\item{flagCell_lm_cutoff}{the cutoff of lm_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile} \item{svmClass_score_cutoff}{the cutoff of transcript score to separate between high and low score transcripts in SVM (default = -2)} @@ -108,6 +109,8 @@ If set to NULL, the pipeline would first randomly choose no more than 2500 cells \item{ctrl_genes}{a vector of control genes that are present in input transcript data.frame but not present in \code{counts} or \code{refProfiles}; the \code{ctrl_genes} would be included in FastReseg analysis. (default = NULL)} +\item{celltype_method}{use either \code{LogLikeRatio} or \code{NegBinomial} method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio)} + \item{filepath_fov_coln}{the column name of each individual file of per FOV transcript data.frame in \code{transDF_fileInfo}} } \value{ diff --git a/man/findSegmentError_allFiles.Rd b/man/findSegmentError_allFiles.Rd index 2e9de7b..bdcbe59 100644 --- a/man/findSegmentError_allFiles.Rd +++ b/man/findSegmentError_allFiles.Rd @@ -21,12 +21,13 @@ findSegmentError_allFiles( spatLocs_colns = c("x", "y", "z"), extracellular_cellID = NULL, flagModel_TransNum_cutoff = 50, - flagCell_lrtest_cutoff = 5, + flagCell_lm_cutoff = 5, svmClass_score_cutoff = -2, svm_args = list(kernel = "radial", scale = FALSE, gamma = 0.4), path_to_output = "reSeg_res", combine_extra = FALSE, - ctrl_genes = NULL + ctrl_genes = NULL, + celltype_method = "LogLikeRatio" ) } \arguments{ @@ -61,7 +62,7 @@ Notice that some assays like SMI has XY axes swapped between stage and each FOV \item{flagModel_TransNum_cutoff}{the cutoff of transcript number to do spatial modeling for identification of wrongly segmented cells (default = 50)} -\item{flagCell_lrtest_cutoff}{the cutoff of lrtest_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile} +\item{flagCell_lm_cutoff}{the cutoff of lm_-log10P to identify putative wrongly segemented cells with strong spatial dependency in transcript score profile} \item{svmClass_score_cutoff}{the cutoff of transcript score to separate between high and low score transcripts in SVM (default = -2)} @@ -73,6 +74,8 @@ Notice that some assays like SMI has XY axes swapped between stage and each FOV \item{ctrl_genes}{a vector of control genes that are present in input transcript data.frame but not present in \code{counts} or \code{refProfiles}; the \code{ctrl_genes} would be included in FastReseg analysis. (default = NULL)} +\item{celltype_method}{use either \code{LogLikeRatio} or \code{NegBinomial} method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio)} + \item{filepath_fov_coln}{the column name of each individual file of per FOV transcript data.frame in \code{transDF_fileInfo}} } \value{ diff --git a/man/get_baselineCT.Rd b/man/get_baselineCT.Rd index 227b3a3..0eea1e0 100644 --- a/man/get_baselineCT.Rd +++ b/man/get_baselineCT.Rd @@ -4,7 +4,12 @@ \alias{get_baselineCT} \title{get_baselineCT} \usage{ -get_baselineCT(refProfiles, counts, clust = NULL) +get_baselineCT( + refProfiles, + counts, + clust = NULL, + celltype_method = "LogLikeRatio" +) } \arguments{ \item{refProfiles}{A matrix of cluster profiles, genes X clusters} @@ -12,6 +17,8 @@ get_baselineCT(refProfiles, counts, clust = NULL) \item{counts}{Counts matrix, cells X genes.} \item{clust}{Vector of cluster assignments for each cell in \code{counts}, default = NULL to automatically assign the cell cluster for each cell based on maximum transcript score} + +\item{celltype_method}{use either \code{LogLikeRatio} or \code{NegBinomial} method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio)} } \value{ a list diff --git a/man/lldist.Rd b/man/lldist.Rd new file mode 100644 index 0000000..9b77b9a --- /dev/null +++ b/man/lldist.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/quick_celltyping.R +\name{lldist} +\alias{lldist} +\title{lldist} +\usage{ +lldist(x, mat, bg = 0.01, size = 10, digits = 2) +} +\arguments{ +\item{x}{a vector of a reference cell cluster} + +\item{mat}{a cells x genes matrix of expression levels in all cells} + +\item{bg}{a vector for background level of each cell (default: 0.01)} + +\item{size}{the parameters for dnbinom function (default: 10)} + +\item{digits}{the number of digits for rounding} +} +\value{ +a cell x cell_type matrix of the log-likelihood +} +\description{ +calculate the log-likelihood of cell belong to certain cell cluster given the reference profile using negative binomial model +} diff --git a/man/neighborhood_for_resegment_spatstat.Rd b/man/neighborhood_for_resegment_spatstat.Rd index b81c5ad..588e69a 100644 --- a/man/neighborhood_for_resegment_spatstat.Rd +++ b/man/neighborhood_for_resegment_spatstat.Rd @@ -6,7 +6,8 @@ \usage{ neighborhood_for_resegment_spatstat( chosen_cells = NULL, - score_GeneMatrix, + score_GeneMatrix = NULL, + refProfiles = NULL, score_baseline = NULL, neighbor_distance_xy = NULL, distance_cutoff = 2.7, @@ -15,15 +16,18 @@ neighborhood_for_resegment_spatstat( celltype_coln = "cell_type", transID_coln = "transcript_id", transGene_coln = "target", - transSpatLocs_coln = c("x", "y", "z") + transSpatLocs_coln = c("x", "y", "z"), + celltype_method = "LogLikeRatio" ) } \arguments{ \item{chosen_cells}{the cell_ID of chosen cells need to be evaluate for re-segmentation} -\item{score_GeneMatrix}{the gene x cell-type matrix of log-like score of gene in each cell type} +\item{score_GeneMatrix}{the gene x cell-type matrix of log-like score of gene in each cell type, needed if using \code{LogLikeRatio} cell typing method (default = NULL)} -\item{score_baseline}{a named vector of score baseline for all cell type listed in score_GeneMatrix} +\item{refProfiles}{A matrix of cluster profiles, genes X clusters, needed if using \code{NegBionomial} cell typing method (default = NULL)} + +\item{score_baseline}{a named vector of score baseline for all cell type listed in \code{score_GeneMatrix} or \code{refProfiles}} \item{neighbor_distance_xy}{maximum cell-to-cell distance in x, y between the center of query cells to the center of neighbor cells with direct contact, same unit as input spatial coordinate. Default = NULL to use the 2 times of average 2D cell diameter.} @@ -41,6 +45,8 @@ If set to NULL, the pipeline would first randomly choose no more than 2500 cells \item{transGene_coln}{the column name of target or gene name in transcript_df} \item{transSpatLocs_coln}{the column name of 1st, 2nd, optional 3rd spatial dimension of each transcript in transcript_df} + +\item{celltype_method}{use either \code{LogLikeRatio} or \code{NegBinomial} method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio)} } \value{ a data.frame diff --git a/man/numCores.Rd b/man/numCores.Rd new file mode 100644 index 0000000..3706905 --- /dev/null +++ b/man/numCores.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/quick_celltyping.R +\name{numCores} +\alias{numCores} +\title{Get number of cores for parallelized operations} +\usage{ +numCores() +} +\value{ +number of cores to use for mclapply +} +\description{ +Get number of cores for parallelized operations +} diff --git a/man/prepSMI_for_fastReseg.Rd b/man/prepSMI_for_fastReseg.Rd index 1bf5608..e20afc4 100644 --- a/man/prepSMI_for_fastReseg.Rd +++ b/man/prepSMI_for_fastReseg.Rd @@ -11,7 +11,8 @@ prepSMI_for_fastReseg( cellClus_to_exclude = NULL, removeUnpaired = FALSE, blacklist_genes = NULL, - pixel_size = 0.18 + pixel_size = 0.18, + celltype_method = "LogLikeRatio" ) } \arguments{ @@ -28,6 +29,8 @@ prepSMI_for_fastReseg( \item{blacklist_genes}{a vector of genes to be excluded from reference profile estimation (default = NULL)} \item{pixel_size}{the micrometer size of image pixel listed in \code{Width} and \code{Height} dimension of each cell stored in \code{cell_metadata} of the existing SMI object (default = 0.18)} + +\item{celltype_method}{use either \code{LogLikeRatio} or \code{NegBinomial} method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio)} } \value{ a list diff --git a/man/quick_celltype.Rd b/man/quick_celltype.Rd new file mode 100644 index 0000000..64fe80f --- /dev/null +++ b/man/quick_celltype.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/quick_celltyping.R +\name{quick_celltype} +\alias{quick_celltype} +\title{quick_celltype} +\usage{ +quick_celltype( + x, + bg = 0.01, + reference_profiles, + nb_size = 10, + align_genes = TRUE +) +} +\arguments{ +\item{x}{Counts matrix (or dgCMatrix), cells x genes.} + +\item{bg}{a vector for background level of each cell (default = 0.01)} + +\item{reference_profiles}{Matrix of expression profiles of pre-defined clusters, genes x clusters.} + +\item{nb_size}{The size parameter to assume for the NB distribution.} + +\item{align_genes}{Logical, for whether to align the counts matrix and the reference_profiles by gene ID.} +} +\value{ +A list, with the following elements: +\enumerate{ +\item clust: a vector given cells' cluster assignments +\item logliks: Matrix of cells' log-likelihoods under each cluster. Cells in rows, clusters in columns. +} +} +\description{ +Classify cells based on reference profiles given the maximum log-likelihood calculated via negative binomial model +} diff --git a/man/score_cell_segmentation_error.Rd b/man/score_cell_segmentation_error.Rd index d68f8ad..8b62650 100644 --- a/man/score_cell_segmentation_error.Rd +++ b/man/score_cell_segmentation_error.Rd @@ -34,14 +34,14 @@ data.frame with columns for \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} } } \description{ Score each cell for how much their transcripts change their goodness-of-fit over space. } \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. +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. } diff --git a/man/update_transDF_ResegActions.Rd b/man/update_transDF_ResegActions.Rd index 6f1a49b..5cc2dcd 100644 --- a/man/update_transDF_ResegActions.Rd +++ b/man/update_transDF_ResegActions.Rd @@ -7,12 +7,14 @@ update_transDF_ResegActions( transcript_df, reseg_full_converter, - score_GeneMatrix, + score_GeneMatrix = NULL, + refProfiles = NULL, transGene_coln = "target", cellID_coln = "cell_ID", celltype_coln = "cell_type", spatLocs_colns = c("x", "y", "z"), - return_perCellDF = TRUE + return_perCellDF = TRUE, + celltype_method = "LogLikeRatio" ) } \arguments{ @@ -20,7 +22,9 @@ update_transDF_ResegActions( \item{reseg_full_converter}{a named converter to update the cell ID in \code{transcript_df}, cell_ID in name would be converted to cell_ID in value; discard cell_ID with value = NA} -\item{score_GeneMatrix}{a gene x cell-type score matrix} +\item{score_GeneMatrix}{the gene x cell-type matrix of log-like score of gene in each cell type, needed if using \code{LogLikeRatio} cell typing method (default = NULL)} + +\item{refProfiles}{A matrix of cluster profiles, genes X clusters, needed if using \code{NegBionomial} cell typing method (default = NULL)} \item{transGene_coln}{the column name of target or gene name in \code{transcript_df}} @@ -31,6 +35,8 @@ update_transDF_ResegActions( \item{spatLocs_colns}{column names for 1st, 2nd and optional 3rd dimension of spatial coordinates in transcript_df} \item{return_perCellDF}{flag to return gene x cell count matrix and per cell DF with updated mean spatial coordinates and new cell type} + +\item{celltype_method}{use either \code{LogLikeRatio} or \code{NegBinomial} method for quick cell typing and corresponding score_baseline calculation (default = LogLikeRatio)} } \value{ a list From cebdc7997c22f2f6f8f96bb85f8aeabb21d80766 Mon Sep 17 00:00:00 2001 From: lidanwu Date: Fri, 21 Oct 2022 02:27:55 +0000 Subject: [PATCH 04/10] fix typo --- R/get_baseline.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/get_baseline.R b/R/get_baseline.R index ab34a67..635270d 100644 --- a/R/get_baseline.R +++ b/R/get_baseline.R @@ -268,7 +268,7 @@ get_baselineCT <- function(refProfiles, } - } else if (celltype_method == ' NegBinomial'){ + } else if (celltype_method == 'NegBinomial'){ # assign cell type for each cell if not provided ---- if(is.null(clust)){ From 88ebca6a06c77f0d1abd03c4ed2dc4b8e40a1e1f Mon Sep 17 00:00:00 2001 From: lidanwu Date: Fri, 21 Oct 2022 02:57:39 +0000 Subject: [PATCH 05/10] update baseline calc --- R/get_baseline.R | 29 +++++++++++++---------------- 1 file changed, 13 insertions(+), 16 deletions(-) diff --git a/R/get_baseline.R b/R/get_baseline.R index 635270d..73683e0 100644 --- a/R/get_baseline.R +++ b/R/get_baseline.R @@ -255,44 +255,41 @@ get_baselineCT <- function(refProfiles, # 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 ---- + # 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] - - common_celltypes <- unique(clust) + 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) - # assign cell type for each cell if not provided ---- + # per cell logliks for all cells, cell x cell-cluster score matrix + tLLRv2_cellMatrix <- nb_res[['logliks']] + + # 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`.') - - nb_res <- quick_celltype(counts, bg = 0.01, - reference_profiles = refProfiles, - align_genes = FALSE) - - clust <- nb_res[['clust']] - common_celltypes <- unique(clust) } - - # per cell logliks for all cells, cell x cell-cluster score matrix - tLLRv2_cellMatrix <- nb_res[['logliks']] - + rm(nb_res) + } 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) From de7a0788259a6afd93ddd5ec7b389089d869783c Mon Sep 17 00:00:00 2001 From: lidanwu Date: Fri, 21 Oct 2022 17:22:45 +0000 Subject: [PATCH 06/10] deal with no altered cell --- R/update_segments.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/R/update_segments.R b/R/update_segments.R index 4e64e72..0ad02f9 100644 --- a/R/update_segments.R +++ b/R/update_segments.R @@ -92,8 +92,9 @@ update_transDF_ResegActions <- function(transcript_df, ## get new cell types for cells being updated ---- subTransDF <- transcript_df[which(transcript_df[['updated_cellID']] %in% unique(cells_to_update) & transcript_df[[transGene_coln]] %in% common_genes), ] - - if(celltype_method == 'LogLikeRatio'){ + if(nrow(subTransDF)<1){ + newCellTypes <- NULL + } else if(celltype_method == 'LogLikeRatio'){ # get score for each transcripts transcriptGeneScore <- score_GeneMatrix[subTransDF[[transGene_coln]], ] From 7171cc6b4630992592b5055e7950fd0f369da8f1 Mon Sep 17 00:00:00 2001 From: lidanwu Date: Sat, 22 Oct 2022 02:03:10 +0000 Subject: [PATCH 07/10] same gene set for all negBinomial cell typing --- R/neighborhood_spatstat.R | 11 +++++++++-- R/update_segments.R | 11 +++++++++-- R/wrapper_FastReseg.R | 35 +++++++++++++++++++++++++++++------ 3 files changed, 47 insertions(+), 10 deletions(-) diff --git a/R/neighborhood_spatstat.R b/R/neighborhood_spatstat.R index 5e22176..bbd07f5 100644 --- a/R/neighborhood_spatstat.R +++ b/R/neighborhood_spatstat.R @@ -201,8 +201,15 @@ neighborhood_for_resegment_spatstat <- function(chosen_cells = NULL, # get cell type and logliks for chosen cells only if(celltype_method == 'NegBinomial'){ exprMat <- reshape2::acast(chosen_transDF, as.formula(paste(cellID_coln, '~', transGene_coln)), length) - nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = TRUE) - rm(exprMat) + # fill missing genes that in refProfiles but not in current data as 0 + missingGenes <- setdiff(rownames(refProfiles), colnames(exprMat)) + exprMat <- cbind(exprMat, + matrix(0, nrow = nrow(exprMat), ncol = length(missingGenes), + dimnames = list(rownames(exprMat), missingGenes))) + exprMat <- exprMat[, rownames(refProfiles), drop = FALSE] + + nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = FALSE) + rm(exprMat, missingGenes) } ## get molecular_distance_cutoff between neighbor cells from 10 randomly selected ROIs with 5* neighbor_distance_xy diff --git a/R/update_segments.R b/R/update_segments.R index 0ad02f9..517367f 100644 --- a/R/update_segments.R +++ b/R/update_segments.R @@ -111,8 +111,15 @@ update_transDF_ResegActions <- function(transcript_df, } else if (celltype_method =='NegBinomial'){ exprMat <- reshape2::acast(subTransDF, as.formula(paste('updated_cellID', '~', transGene_coln)), length) - newCellTypes <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = TRUE)[['clust']] - rm(exprMat) + # fill missing genes that in refProfiles but not in current data as 0 + missingGenes <- setdiff(rownames(refProfiles), colnames(exprMat)) + exprMat <- cbind(exprMat, + matrix(0, nrow = nrow(exprMat), ncol = length(missingGenes), + dimnames = list(rownames(exprMat), missingGenes))) + exprMat <- exprMat[, rownames(refProfiles), drop = FALSE] + + newCellTypes <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = FALSE)[['clust']] + rm(exprMat, missingGenes) } diff --git a/R/wrapper_FastReseg.R b/R/wrapper_FastReseg.R index 4f1dcdb..369392d 100644 --- a/R/wrapper_FastReseg.R +++ b/R/wrapper_FastReseg.R @@ -292,10 +292,17 @@ fastReseg_core_externalRef <- function(refProfiles, } else if (celltype_method =='NegBinomial'){ # `quick_celltype` function returns a list contains element `clust`, a vector given cells' cluster assignments. exprMat <- reshape2::acast(transcript_df, as.formula(paste0(cellID_coln, "~", transGene_coln)), length) - nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = TRUE) + # fill missing genes that in refProfiles but not in current data as 0 + missingGenes <- setdiff(rownames(refProfiles), colnames(exprMat)) + exprMat <- cbind(exprMat, + matrix(0, nrow = nrow(exprMat), ncol = length(missingGenes), + dimnames = list(rownames(exprMat), missingGenes))) + exprMat <- exprMat[, rownames(refProfiles), drop = FALSE] + + nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = FALSE) select_cellmeta <- data.frame(cellID = names(nb_res[['clust']]), celltype = nb_res[['clust']]) - rm(exprMat, nb_res) + rm(exprMat, nb_res, missingGenes) } @@ -557,10 +564,17 @@ fastReseg_core_externalRef <- function(refProfiles, } else if (celltype_method =='NegBinomial'){ # `quick_celltype` function returns a list contains element `clust`, a vector given cells' cluster assignments. exprMat <- reshape2::acast(flagged_transDF_SVM3, as.formula(paste0("tmp_cellID", "~", transGene_coln)), length) - nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = TRUE) + # fill missing genes that in refProfiles but not in current data as 0 + missingGenes <- setdiff(rownames(refProfiles), colnames(exprMat)) + exprMat <- cbind(exprMat, + matrix(0, nrow = nrow(exprMat), ncol = length(missingGenes), + dimnames = list(rownames(exprMat), missingGenes))) + exprMat <- exprMat[, rownames(refProfiles), drop = FALSE] + + nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = FALSE) tmp_df <- list(cellType_DF = data.frame(cellID = names(nb_res[['clust']]), celltype = nb_res[['clust']])) - rm(exprMat, nb_res) + rm(exprMat, nb_res, missingGenes) } @@ -1688,6 +1702,8 @@ findSegmentError_allFiles <- function(counts, stop("Too few common genes between the `refProfiles` (genes X clusters) and `counts` (cells X genes), check if correct format. ") } + refProfiles <- refProfiles[common_genes, ] + ## initialize list to collect each FOV outputs ---- all_segRes <- list() @@ -1786,10 +1802,17 @@ findSegmentError_allFiles <- function(counts, } else if (celltype_method =='NegBinomial'){ # `quick_celltype` function returns a list contains element `clust`, a vector given cells' cluster assignments. exprMat <- reshape2::acast(transcript_df[['intraC']], as.formula(paste0("UMI_cellID", "~", "target")), length) - nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = TRUE) + # fill missing genes that in refProfiles but not in current data as 0 + missingGenes <- setdiff(rownames(refProfiles), colnames(exprMat)) + exprMat <- cbind(exprMat, + matrix(0, nrow = nrow(exprMat), ncol = length(missingGenes), + dimnames = list(rownames(exprMat), missingGenes))) + exprMat <- exprMat[, rownames(refProfiles), drop = FALSE] + + nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = FALSE) select_cellmeta <- data.frame(cellID = names(nb_res[['clust']]), celltype = nb_res[['clust']]) - rm(exprMat, nb_res) + rm(exprMat, nb_res, missingGenes) } From 6dec36b246930a6f5b5bd118810f656f8b9ffad5 Mon Sep 17 00:00:00 2001 From: lidanwu Date: Sat, 22 Oct 2022 02:07:17 +0000 Subject: [PATCH 08/10] fix 1 cell negbinomal type --- R/quick_celltyping.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/R/quick_celltyping.R b/R/quick_celltyping.R index af290ec..24674e1 100644 --- a/R/quick_celltyping.R +++ b/R/quick_celltyping.R @@ -40,10 +40,11 @@ lldist <- function(x, mat, bg = 0.01, size = 10, digits = 2) { } # calc scaling factor to put y on the scale of x: - if (is.vector(bg)) { + if (is.vector(bg) & nrow(mat) >1) { bgsub <- apply(mat, 2, function(xx) xx - bg) bgsub <- pmax(bgsub, 0) - } else { + } else { + # bg is a matrix or mat only has 1 row bgsub <- pmax(mat - bg, 0) } @@ -135,8 +136,7 @@ quick_celltype <- function(x, bg = 0.01, reference_profiles, nb_size = 10, align size = nb_size, mc.cores = numCores()) logliks <- do.call(cbind, logliks) - - + # get remaining outputs clust <- colnames(logliks)[apply(logliks, 1, which.max)] From 2215c7a5cf7b72aa8257cff632156fbf79efdcac Mon Sep 17 00:00:00 2001 From: lidanwu Date: Sat, 22 Oct 2022 19:35:36 +0000 Subject: [PATCH 09/10] fix cells of 0 counts --- R/get_baseline.R | 6 +++--- R/neighborhood_spatstat.R | 23 ++++++++++++++++---- R/quick_celltyping.R | 45 ++++++++++++++++++++++++++++++++++----- R/update_segments.R | 12 +++++++++-- R/wrapper_FastReseg.R | 41 +++++++++++++++++++++++++++++++---- man/quick_celltype.Rd | 7 +++--- 6 files changed, 113 insertions(+), 21 deletions(-) diff --git a/R/get_baseline.R b/R/get_baseline.R index 73683e0..aebaa1d 100644 --- a/R/get_baseline.R +++ b/R/get_baseline.R @@ -273,13 +273,13 @@ get_baselineCT <- function(refProfiles, reference_profiles = refProfiles, align_genes = FALSE) - # per cell logliks for all cells, cell x cell-cluster score matrix - tLLRv2_cellMatrix <- nb_res[['logliks']] + # 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']] + clust <- nb_res[['clust']][1: (length(nb_res[['clust']]) - length(nb_res[['zeroCells']]))] } rm(nb_res) diff --git a/R/neighborhood_spatstat.R b/R/neighborhood_spatstat.R index bbd07f5..2c1e428 100644 --- a/R/neighborhood_spatstat.R +++ b/R/neighborhood_spatstat.R @@ -209,6 +209,19 @@ neighborhood_for_resegment_spatstat <- function(chosen_cells = NULL, exprMat <- exprMat[, rownames(refProfiles), drop = FALSE] nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = FALSE) + + # transcript groups without informative genes would use the original cluster assignment + if(!is.null(nb_res[['zeroCells']])){ + oldCT_df <- unique(chosen_transDF[get(cellID_coln) %in% nb_res[['zeroCells']], + .SD, .SDcols = c(cellID_coln, celltype_coln)]) + oldCT_vector <- oldCT_df[[celltype_coln]] + names(oldCT_vector) <- oldCT_df[[cellID_coln]] + nb_res[['clust']] <- c(nb_res[['clust']][1: (length(nb_res[['clust']]) - length(nb_res[['zeroCells']]))], + oldCT_vector) + rm(oldCT_df, oldCT_vector) + + } + rm(exprMat, missingGenes) } @@ -416,17 +429,19 @@ neighborhood_for_resegment_spatstat <- function(chosen_cells = NULL, cell_score <- colSums(cell_score) } cell_score <- matrix(cell_score, nrow = 1, dimnames = list(each_cell, colnames(score_GeneMatrix))) + maxCT_1st <- colnames(cell_score)[max.col(cell_score, ties.method="first")] + } else if (celltype_method =='NegBinomial'){ cell_score <- nb_res[['logliks']][each_cell, , drop = FALSE] - + maxCT_1st <- nb_res[['clust']][each_cell] } - max_idx_1st <- max.col(cell_score, ties.method="first") + queryPerCell_df <- data.frame(CellId = each_cell, cell_type = query_df[[celltype_coln]][1], transcript_num = nrow(query_df), - self_celltype = colnames(cell_score)[max_idx_1st], - score_under_self = cell_score[each_cell, max_idx_1st]/nrow(query_df)) + self_celltype = maxCT_1st, + score_under_self = cell_score[each_cell, maxCT_1st]/nrow(query_df)) # get neighbor cell type and check the query transcript cell type against it diff --git a/R/quick_celltyping.R b/R/quick_celltyping.R index 24674e1..3f8aad5 100644 --- a/R/quick_celltyping.R +++ b/R/quick_celltyping.R @@ -91,15 +91,35 @@ numCores <- function() { #' @param nb_size The size parameter to assume for the NB distribution. #' @param align_genes Logical, for whether to align the counts matrix and the reference_profiles by gene ID. #' @return A list, with the following elements: -#' \enumerate{ -#' \item clust: a vector given cells' cluster assignments -#' \item logliks: Matrix of cells' log-likelihoods under each cluster. Cells in rows, clusters in columns. +#' \describe{ +#' \item{clust}{a vector given cells' cluster assignments, return NA for cells of zero counts. } +#' \item{logliks}{a cells x clusters matrix of cells' log-likelihoods under each cluster, return -Inf for cells of zero counts. } +#' \item{zeroCells}{a vector of cells of zero count, return NULL if none} #' } #' @export quick_celltype <- function(x, bg = 0.01, reference_profiles, nb_size = 10, align_genes = TRUE) { if (any(rowSums(x) == 0)) { - stop("Cells with 0 counts were found. Please remove.") + zeroCells <- rownames(x)[rowSums(x)==0] + message(sprintf("%d cells with 0 counts are found. Return `clust = NA`, `logliks = -Inf` for those cells: `%s`.", + length(zeroCells), paste0(zeroCells, collapse = "`, `"))) + + # bg is a vector of same length as x + if(is.vector(bg)){ + if(identical(length(bg), nrow(x))){ + bg <- bg[rowSums(x)!=0] + } + + } else { + # bg is a cell x gene matrix + bg <- bg[rowSums(x)!=0, ] + } + + # remove zero cells from x + x <- x[rowSums(x)!=0, ] + + } else { + zeroCells <- NULL } # accept a single value of bg if input by user: @@ -108,6 +128,7 @@ quick_celltype <- function(x, bg = 0.01, reference_profiles, nb_size = 10, align names(bg) <- rownames(x) } + # align genes: if (align_genes) { sharedgenes <- intersect(rownames(reference_profiles), colnames(x)) @@ -143,6 +164,20 @@ quick_celltype <- function(x, bg = 0.01, reference_profiles, nb_size = 10, align names(clust) <- rownames(logliks) out <- list(clust = clust, - logliks = round(logliks, 4)) + logliks = round(logliks, 4), + zeroCells = zeroCells) + + # add in zeroCell data + if(!is.null(zeroCells)){ + clust <- rep(NA, length(zeroCells)) + names(clust) <- zeroCells + + logliks <- matrix(-Inf, nrow = length(zeroCells), ncol = ncol(out[['logliks']]), + dimnames = list(zeroCells, colnames(out[['logliks']]))) + + out[['clust']] <- c(out[['clust']], clust) + out[['logliks']] <- rbind(out[['logliks']], logliks) + } + return(out) } diff --git a/R/update_segments.R b/R/update_segments.R index 517367f..62fd71a 100644 --- a/R/update_segments.R +++ b/R/update_segments.R @@ -118,8 +118,16 @@ update_transDF_ResegActions <- function(transcript_df, dimnames = list(rownames(exprMat), missingGenes))) exprMat <- exprMat[, rownames(refProfiles), drop = FALSE] - newCellTypes <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = FALSE)[['clust']] - rm(exprMat, missingGenes) + nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = FALSE) + + # transcript groups without informative genes would use the original cluster assignment + if(!is.null(nb_res[['zeroCells']])){ + cells_to_update <- setdiff(cells_to_update, nb_res[['zeroCells']]) + } + + newCellTypes <- nb_res[['clust']] + + rm(exprMat, missingGenes, nb_res) } diff --git a/R/wrapper_FastReseg.R b/R/wrapper_FastReseg.R index 369392d..5975dfc 100644 --- a/R/wrapper_FastReseg.R +++ b/R/wrapper_FastReseg.R @@ -302,6 +302,16 @@ fastReseg_core_externalRef <- function(refProfiles, nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = FALSE) select_cellmeta <- data.frame(cellID = names(nb_res[['clust']]), celltype = nb_res[['clust']]) + + # transcript groups without informative genes would be assigned with 1st cell type in refProfiles + if(!is.null(nb_res[['zeroCells']])){ + message(sprintf("Found %d cells of zero informative counts in original transcript_df, assign initial cell type = `%s`: `%s`.", + length(nb_res[['zeroCells']]), colnames(refProfiles)[1], + paste0(nb_res[['zeroCells']], collapse = "`, `"))) + + select_cellmeta[['celltype']][is.na(select_cellmeta[['celltype']])] <- colnames(refProfiles)[1] + } + rm(exprMat, nb_res, missingGenes) } @@ -561,6 +571,8 @@ fastReseg_core_externalRef <- function(refProfiles, transGene_coln = transGene_coln, cellID_coln = "tmp_cellID", return_transMatrix = FALSE) + tmp_df <- tmp_df[['cellType_DF']] + } else if (celltype_method =='NegBinomial'){ # `quick_celltype` function returns a list contains element `clust`, a vector given cells' cluster assignments. exprMat <- reshape2::acast(flagged_transDF_SVM3, as.formula(paste0("tmp_cellID", "~", transGene_coln)), length) @@ -572,15 +584,26 @@ fastReseg_core_externalRef <- function(refProfiles, exprMat <- exprMat[, rownames(refProfiles), drop = FALSE] nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = FALSE) - tmp_df <- list(cellType_DF = data.frame(cellID = names(nb_res[['clust']]), - celltype = nb_res[['clust']])) + + tmp_df <- data.frame(tmp_cellID = names(nb_res[['clust']]), + SVM_cell_type = nb_res[['clust']]) + + # transcript groups without informative genes would use the original cluster assignment + if(!is.null(nb_res[['zeroCells']])){ + oldCT_df <- unique(flagged_transDF_SVM3[tmp_cellID %in% nb_res[['zeroCells']], + .SD, .SDcols = c('tmp_cellID','SVM_cell_type')]) + + tmp_df <- rbind(tmp_df[!is.na(tmp_df[['SVM_cell_type']]), ], oldCT_df) + rm(oldCT_df) + } + rm(exprMat, nb_res, missingGenes) } - colnames(tmp_df[['cellType_DF']]) <- c('tmp_cellID','group_maxCellType') - flagged_transDF_SVM3 <- merge(flagged_transDF_SVM3, tmp_df[['cellType_DF']], by = 'tmp_cellID', all.x = TRUE) + colnames(tmp_df) <- c('tmp_cellID','group_maxCellType') + flagged_transDF_SVM3 <- merge(flagged_transDF_SVM3, tmp_df, by = 'tmp_cellID', all.x = TRUE) flagged_transDF_SVM3 <- as.data.frame(flagged_transDF_SVM3) rm(tmp_df) @@ -1812,6 +1835,16 @@ findSegmentError_allFiles <- function(counts, nb_res <- quick_celltype(exprMat, bg = 0, reference_profiles = refProfiles, align_genes = FALSE) select_cellmeta <- data.frame(cellID = names(nb_res[['clust']]), celltype = nb_res[['clust']]) + + # transcript groups without informative genes would be assigned with 1st cell type in refProfiles + if(!is.null(nb_res[['zeroCells']])){ + message(sprintf("Found %d cells of zero informative counts in original transcript_df, assign initial cell type = `%s`: `%s`.", + length(nb_res[['zeroCells']]), colnames(refProfiles)[1], + paste0(nb_res[['zeroCells']], collapse = "`, `"))) + + select_cellmeta[['celltype']][is.na(select_cellmeta[['celltype']])] <- colnames(refProfiles)[1] + } + rm(exprMat, nb_res, missingGenes) } diff --git a/man/quick_celltype.Rd b/man/quick_celltype.Rd index 64fe80f..c196c3b 100644 --- a/man/quick_celltype.Rd +++ b/man/quick_celltype.Rd @@ -25,9 +25,10 @@ quick_celltype( } \value{ A list, with the following elements: -\enumerate{ -\item clust: a vector given cells' cluster assignments -\item logliks: Matrix of cells' log-likelihoods under each cluster. Cells in rows, clusters in columns. +\describe{ +\item{clust}{a vector given cells' cluster assignments, return NA for cells of zero counts. } +\item{logliks}{a cells x clusters matrix of cells' log-likelihoods under each cluster, return -Inf for cells of zero counts. } +\item{zeroCells}{a vector of cells of zero count, return NULL if none} } } \description{ From 97fc7a236c4893d4716c1f86b3b17f1fc129329c Mon Sep 17 00:00:00 2001 From: lidanwu Date: Sat, 22 Oct 2022 20:03:54 +0000 Subject: [PATCH 10/10] cleanup refProfiles --- R/wrapper_FastReseg.R | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/R/wrapper_FastReseg.R b/R/wrapper_FastReseg.R index 5975dfc..73e4622 100644 --- a/R/wrapper_FastReseg.R +++ b/R/wrapper_FastReseg.R @@ -963,7 +963,7 @@ fastReseg_core_externalRef <- function(refProfiles, #' @export #' fastReseg_internalRef <- function(counts, - clust, + clust = NULL, refProfiles = NULL, transDF_fileInfo = NULL, filepath_coln = 'file_path', @@ -1097,6 +1097,8 @@ fastReseg_internalRef <- function(counts, clust = as.character(clust), s = Matrix::rowSums(as.matrix(counts)), bg = rep(0, nrow(counts))) + }else { + refProfiles <- refProfiles[intersect(rownames(refProfiles), colnames(counts)), ] } # # `get_baselineCT` function gets cluster-specific quantile distribution of transcript number and per cell per molecule transcript score in the provided cell x gene expression matrix based on the reference profiles and cell cluster assignment. @@ -1111,6 +1113,8 @@ fastReseg_internalRef <- function(counts, } else { # reference profiles exists, but no cluster assignment + refProfiles <- refProfiles[intersect(rownames(refProfiles), colnames(counts)), ] + baselineData <- get_baselineCT(refProfiles = refProfiles, counts = counts, clust = NULL, celltype_method = celltype_method) clust = baselineData[['clust_used']] }