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/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/get_baseline.R b/R/get_baseline.R index 1379a11..aebaa1d 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,56 @@ get_baselineCT <- function(refProfiles, counts <- as.matrix(counts)[, common_genes] # get score matrix based on refProfiles for each gene and cell ---- - # replace zero in mean profiles with 1E-5 - refProfiles <- pmax(refProfiles, 1e-5) - # tLL score - transcript_loglik <- scoreGenesInRef(genes = common_genes, ref_profiles = refProfiles) - # tLLR score, re-center on maximum per row/transcript - tmp_max <- apply(transcript_loglik, 1, max) - tLLRv2_geneMatrix <- sweep(transcript_loglik, 1, tmp_max, '-') - rm(tmp_max, transcript_loglik) - - # get cell x cell-cluster score matrix = counts (cell x gene) %*% tLLR_score (gene x cell-cluster) - tLLRv2_cellMatrix <- counts %*% tLLRv2_geneMatrix + # replace zero in mean profiles with 1E-8 + refProfiles <- pmax(refProfiles, 1e-8) - # assign cell type for each cell if not provided ---- - if(is.null(clust)){ - message('Perform cluster assignment based on maximum transcript score given the provided `refProfiles`.') + if(celltype_method == 'LogLikeRatio'){ + # tLL score + transcript_loglik <- scoreGenesInRef(genes = common_genes, ref_profiles = refProfiles) + # tLLR score, re-center on maximum per row/transcript + tmp_max <- apply(transcript_loglik, 1, max) + tLLRv2_geneMatrix <- sweep(transcript_loglik, 1, tmp_max, '-') + rm(tmp_max, transcript_loglik) + + # get cell x cell-cluster score matrix = counts (cell x gene) %*% tLLR_score (gene x cell-cluster) + tLLRv2_cellMatrix <- counts %*% tLLRv2_geneMatrix + + # assign cell type for each cell if not provided + if(is.null(clust)){ + message('Perform cluster assignment based on maximum transcript score given the provided `refProfiles`.') + + # assign cell type based on max values + max_idx_1st <- max.col(tLLRv2_cellMatrix, ties.method="first") + clust <- colnames(tLLRv2_cellMatrix)[max_idx_1st] - # assign cell type based on max values - max_idx_1st <- max.col(tLLRv2_cellMatrix, ties.method="first") - clust <- colnames(tLLRv2_cellMatrix)[max_idx_1st] + rm(max_idx_1st) + } + + + } else if (celltype_method == 'NegBinomial'){ + # get logliks for cell under all cell types + nb_res <- quick_celltype(counts, bg = 0.01, + reference_profiles = refProfiles, + align_genes = FALSE) + + # per cell logliks for all cells, cell x cell-cluster score matrix, exclude cells of zero count + tLLRv2_cellMatrix <- nb_res[['logliks']][1: (length(nb_res[['clust']]) - length(nb_res[['zeroCells']])), ] + + # assign cell type for each cell if not provided + if(is.null(clust)){ + message('Perform cluster assignment based on negative binomial model given the provided `refProfiles`.') + clust <- nb_res[['clust']][1: (length(nb_res[['clust']]) - length(nb_res[['zeroCells']]))] + + } + rm(nb_res) - common_celltypes <- unique(clust) - rm(max_idx_1st) + + } else { + stop(sprintf('The provided `celltype_method` = `%s` is not supported.', celltype_method)) } + common_celltypes <- unique(clust) + # get transcript number quantile profile --- all_transNum <- rowSums(counts) span_transNum_CellType <- tapply(all_transNum, @@ -274,6 +305,7 @@ get_baselineCT <- function(refProfiles, rowidx <- which(clust == each_celltype) all_tLLRv2[rowidx] <- tLLRv2_cellMatrix[rowidx, each_celltype] } + # normalized by transcript number to get per molecule transcript score for each cell all_tLLRv2 <- all_tLLRv2/all_transNum span_tLLRv2_CellType <- tapply(all_tLLRv2, diff --git a/R/neighborhood_spatstat.R b/R/neighborhood_spatstat.R index 9939753..2c1e428 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,33 @@ 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) + # 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) + + # 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) + } + ## 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,20 +422,26 @@ 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))) + 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] } - 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), - 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 @@ -382,8 +450,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..3f8aad5 --- /dev/null +++ b/R/quick_celltyping.R @@ -0,0 +1,183 @@ +#' @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 +#' @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)) { + 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) & nrow(mat) >1) { + bgsub <- apply(mat, 2, function(xx) xx - bg) + bgsub <- pmax(bgsub, 0) + } else { + # bg is a matrix or mat only has 1 row + 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: +#' \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)) { + 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: + 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), + 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 67c6a8d..62fd71a 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,46 @@ 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]], ] + if(nrow(subTransDF)<1){ + newCellTypes <- NULL + } else 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) + # 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) + + # 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']] - 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 + rm(exprMat, missingGenes, nb_res) + + } + + # 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..73e4622 100644 --- a/R/wrapper_FastReseg.R +++ b/R/wrapper_FastReseg.R @@ -9,10 +9,10 @@ #' @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') +#' @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} @@ -65,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, @@ -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,47 @@ 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']] - colnames(select_cellmeta) <- c(cellID_coln,'tLLRv2_maxCellType') - rm(tmp_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) + # 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']]) + + # 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) + } + + + colnames(select_cellmeta) <- c(cellID_coln,'tSum_maxCellType') + 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))) @@ -296,7 +329,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) @@ -307,7 +340,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) @@ -322,7 +355,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) @@ -331,10 +364,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)) } @@ -360,12 +393,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 @@ -431,7 +464,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, @@ -530,15 +563,47 @@ 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) - 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) + 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) + 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) + # 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 <- 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) <- 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) @@ -559,7 +624,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 ---- @@ -571,17 +636,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 +697,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 @@ -714,7 +811,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') @@ -732,6 +829,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} @@ -865,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', @@ -880,7 +978,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, @@ -901,10 +999,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) @@ -995,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. @@ -1005,11 +1109,13 @@ 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) + refProfiles <- refProfiles[intersect(rownames(refProfiles), colnames(counts)), ] + + baselineData <- get_baselineCT(refProfiles = refProfiles, counts = counts, clust = NULL, celltype_method = celltype_method) clust = baselineData[['clust_used']] } @@ -1237,7 +1343,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, @@ -1246,7 +1352,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' @@ -1394,12 +1501,13 @@ 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. #' @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} @@ -1495,14 +1603,18 @@ 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, 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 +1715,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']] } @@ -1617,6 +1729,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() @@ -1633,7 +1747,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 +1814,47 @@ 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) + # 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']]) + + # 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) + + } - select_cellmeta <- tmp_df[['cellType_DF']] - colnames(select_cellmeta) <- c('UMI_cellID','tLLRv2_maxCellType') - rm(tmp_df) + 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') @@ -1724,19 +1867,19 @@ 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) - ## (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']], 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) @@ -1747,16 +1890,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 @@ -1781,7 +1924,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/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..c196c3b --- /dev/null +++ b/man/quick_celltype.Rd @@ -0,0 +1,36 @@ +% 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: +\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{ +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 diff --git a/vignettes/0_flagErrorOnly_on_SMIobject.R b/vignettes/0_flagErrorOnly_on_SMIobject.R index be940e7..eec76d3 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')) @@ -78,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 @@ -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 ) @@ -230,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 ea78d54..c818b8b 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')) @@ -77,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 @@ -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 ) 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){