anna-neufeld
commented
1 year ago Hi @mihem ,
Thank you so much for this extensive work! It is so exciting to see NBCS in action !
Here are a few comments / answers to your questions.
(1) Thank you so much for writing the faster version of the SSE function! The one included in the very basic tutorial was only intended to be used on the small toy dataset, so it's great to have this faster version on hand. I had been thinking about using an alternate loss function for the real data tutorial, but it looks like SSE worked pretty well for you !
(2) My tutorials have all used cell-by-gene matrices. Having the observations in the rows and the features in the columns is nice because it is compatible with existing R functions such as k-means (used in my tutorials). However, you are 100% correct that this is the opposite of the way that Seurat does it (they do gene-by-cell matrices). So, the fact that you took transposes was correct.
(3) I agree that the resolution parameter suggested by Poisson count splitting is too big. This is probably due to the fact that the data are not Poisson, they are negative binomial. Thus, in the Poisson version, there is dependence between the training and test sets, leading to overfitting.
(4) For using sctransform-- I am so glad that the paper GitHub repository was useful in finding this code, and so sorry that I haven't put this code in a tutorial yet! -- The only step I am not sure about is subsetting the data to only genes with overdispersion estimates. I ran into the issue that sctransform::vst() does not always return an overdispersion estimate for every single gene (e.g. it does not return an estimate for genes that were expressed in <5 cells, by default). I personally decided to retain these genes in the analysis and just set their overdisperson parameter = infinity (i.e. if the estimation fails, assume the genes are Poisson). -- It might not matter here since you are subsetting the genes/cells later anyway!
(5) I think it's really cool that NBCS ends up showing a sharper minimum at a certain (reasonable!) resolution parameter. This is probably because it did a better job making independent train/test sets than Poisson count splitting. The plots are beautiful, and it looks to me like you did the multifold version correctly! It is really nice to see that the minimum around resolution=1 was pretty consistent across folds for multifold splitting.
(5) I have been working on a Seurat negative binomial tutorial for the website using this same dataset and similar ideas, but had not finished it yet. I really like your analysis, as I think that this shows a clear comparison between NBCS and Poisson count splitting on real data, while also showing how to use multifold splitting, etc. If I use this analysis as the basis for a new tutorial, can I credit you on the website and in the package? If so, shoot me an email at aneufeld@fredhutch.org so that you can let me know how best to credit you!
mihem
Hi @anna-neufeld
thanks for this very interesting and very much needed package. I think the approach you take is brilliant and your tutorials are extensive. However, I really miss a tutorial that shows how to use countsplit with real single cell data to find the optimal number of clusters. However, for the end users this is much more important than simulated data.
So I tried to apply your tutorial https://anna-neufeld.github.io/countsplit.tutorials/articles/MSE_tutorial.html to Seurat Data.
First of all, I treid to optimize your
cluster.ssefunction because this has to be run multiple times on pretty large SC datasets.Cluster SSE function
``` # original function cluster.sse <- function(trainDat, testDat, clusters.est, eps.train) { totSS <- 0 for (lab in unique(clusters.est)) { if (sum(clusters.est==lab) > 1) { clustdat.test <- testDat[clusters.est==lab,] clustdat.train <- trainDat[clusters.est==lab,] colmeansTrain <- apply(clustdat.train, 2, mean) pred.means <- 1/eps.train*colmeansTrain ss <- apply(1/(1-eps.train)*clustdat.test, 1, function(u) sum((log(u+1)-log(pred.means+1))^2)) totSS <- totSS+sum(ss) } } return(totSS) } #refactored version: vectorize cluster_sse1 <- function(trainDat, testDat, clusters.est, eps.train) { stopifnot(nrow(trainDat) == length(clusters.est)) totSS <- 0 unique_clusters <- unique(clusters.est) inverse_eps_train <- 1 / (1 - eps.train) for (lab in unique_clusters) { cluster_indices <- which(clusters.est == lab) num_points <- length(cluster_indices) if (num_points > 1) { clustdat.test <- testDat[cluster_indices, ] clustdat.train <- trainDat[cluster_indices, ] colmeansTrain <- colMeans(clustdat.train) pred.means <- colmeansTrain * inverse_eps_train clustdat.test.1 <- 1 / (1 - eps.train) * clustdat.test ss <- rowSums((log1p(t(clustdat.test.1)) - log1p(pred.means))^2) totSS <- totSS + sum(ss) } } return(totSS) } #best version #use the transpose of the trainDat and testDat #sum is the same as sum(rowSums(matrix)) cluster_sse2 <- function(trainDat, testDat, clusters.est, eps.train) { stopifnot(ncol(trainDat) == length(clusters.est)) totSS <- 0 unique_clusters <- unique(clusters.est) inverse_eps_train <- 1 / (1 - eps.train) for (lab in unique_clusters) { cluster_indices <- which(clusters.est == lab) num_points <- length(cluster_indices) if (num_points > 1) { clustdat.test <- testDat[cluster_indices] clustdat.train <- trainDat[,cluster_indices] rowMeansTrain <- rowMeans(clustdat.train) pred.means <- rowMeansTrain * inverse_eps_train clustdat.test.1 <- 1 / (1 - eps.train) * clustdat.test ss <- sum((log1p(clustdat.test.1) - log1p(pred.means))^2) totSS <- totSS + ss } } return(totSS) } ```This led to a nearly 5 times faster code in the
pbmc.countsdataset and the same results.Of note, in the
v2I changed rows and columns because it further improved the performance and is also better compatible with thepbmc.countsmatrix.Here's the code that I used for the poisson version of countsplit with the pbmc.counts dataset using Seurat.
Code Poisson
``` #load libraries library(Seurat) library(countsplit) library(Matrix) library(tidyverse) library(patchwork) data(pbmc.counts, package="countsplit.tutorials") #replace underscores with dashes in rownames rownames(pbmc.counts) <- gsub(x = rownames(pbmc.counts), pattern = "_", replacement = "-") #split data into training and test sets with Poisson set.seed(1) split <- countsplit(pbmc.counts) Xtrain <- split[[1]] Xtest <- split[[2]] ## Set up the training set object and subset pbmcTrain <- CreateSeuratObject(counts = Xtrain, min.cells = 3, min.features = 200) pbmcTrain[["percent.mt"]] <- PercentageFeatureSet(pbmcTrain, pattern = "^MT-") pbmcTrain <- subset(pbmcTrain, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5) # same genes and cells in the train set as in the test set rows <- rownames(pbmcTrain) cols <- colnames(pbmcTrain) XtestSubset <- Xtest[rows,cols] XtrainSubset <- Xtrain[rows,cols] #run Seurat pipeline pbmcTrain <- pbmcTrain |> NormalizeData() |> FindVariableFeatures() |> ScaleData() |> RunPCA() |> FindNeighbors(dims = 1:10) |> RunUMAP(dims = 1:10) # define a range of resolutions resRange <- seq(0.2, 3, by = 0.2) resNames <- paste0("RNA_snn_res.", resRange) for (res in resRange) { pbmcTrain <- FindClusters(pbmcTrain, resolution = res) } # calculcate within-cluster sse for each resolution results <- list() for (k in resNames) { results[k] <- cluster_sse2( trainDat = XtrainSubset, testDat = XtestSubset, clusters.est = pbmcTrain@meta.data[[k]], eps.train = 0.5) } #prepare results for plotting results_df <- data.frame(cluster = names(results), sse = as.numeric(unlist(results))) |> dplyr::mutate(sse_norm = (sse - min(sse)) / (max(sse) - min(sse))) #plot results results_df |> ggplot(aes(x = cluster, y = sse_norm, group = 1)) + geom_point() + geom_line() + theme_bw() + ylab("Within-cluster SSE") + xlab("Number of cluster") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) ```Maybe you can check if the code is correct. I had some troubles with determining what is column and what is row. In your simulated dataset the
clusters.esthave the same length asnrow(trainDat)but in the Seurat Dataset rows are genes and columns are cells, so the cluster annotation have the same length asncol(trainDat)and notnrow(trainDat). So to use your original function I had to uset(Xtrain)andt(Xtest)as input or in my refactored version2Xtrain,Xtest. Correct?The within-cluster SSE suggest as Leiden resolution of 1.2 or 1.4 i would say based on the results.
The Seurat team choose a resolution of 0.5. If you look at the UMAP plots (I used the test dataset without further filtering to produce that plot, is that statistically correct?), I would also say that 1.2 and 1.4 is probably overclustered. But that fit's nicely to your simulated results.
I then wanted to use negative binominal and used sctransform as you suggested. I stole some code from here: https://github.com/anna-neufeld/nbcs_paper_simulations/blob/main/Fig8/reproduce_full_analysis/code/countsplit.R
Code Negative binominal
## countsplit with negative binomial #first get overdispersion estimates sctransform_fit <- sctransform::vst(pbmc.counts, residual_type = "none", min_cells = 10) overdisp_estimates <- sctransform_fit$model_pars_fit[,"theta"] #subset to genes with overdispersion estimates pbmcCountsSubset <- pbmc.counts[names(overdisp_estimates),] split2 <- countsplit(t(pbmcCountsSubset), overdisp = overdisp_estimates) XtrainNB <- t(split2[[1]]) XtestNB <- t(split2[[2]]) ## Set up the training set object and subset pbmcTrainNB <- CreateSeuratObject(counts = XtrainNB, min.cells = 3, min.features = 200) pbmcTrainNB[["percent.mt"]] <- PercentageFeatureSet(pbmcTrainNB, pattern = "^MT-") pbmcTrainNB <- subset(pbmcTrainNB, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5) # same genes and cells in the train set as in the test set rowsNB <- rownames(pbmcTrainNB) colsNB <- colnames(pbmcTrainNB) XtestSubsetNB <- XtestNB[rowsNB,colsNB] XtrainSubsetNB <- XtrainNB[rowsNB,colsNB] #run Seurat pipeline pbmcTrainNB <- pbmcTrainNB |> NormalizeData() |> FindVariableFeatures() |> ScaleData() |> RunPCA() |> FindNeighbors(dims = 1:10) |> RunUMAP(dims = 1:10) # define a range of resolutions resRange <- seq(0.2, 3, by = 0.2) resNames <- paste0("RNA_snn_res.", resRange) for (res in resRange) { pbmcTrainNB <- FindClusters(pbmcTrainNB, resolution = res) } # calculcate within-cluster sse for each resolution resultsNB <- list() for (k in resNames) { resultsNB[k] <- cluster_sse2( trainDat = XtrainSubsetNB, testDat = XtestSubsetNB, clusters.est = pbmcTrainNB@meta.data[[k]], eps.train = 0.5) } #prepare results for plotting results_df_nb <- data.frame(cluster = names(resultsNB), sse = as.numeric(unlist(resultsNB))) |> dplyr::mutate(sse_norm = (sse - min(sse)) / (max(sse) - min(sse))) #plot results results_df_nb |> ggplot(aes(x = cluster, y = sse_norm, group = 1)) + geom_point() + geom_line() + theme_bw() + ylab("Within-cluster SSE") + xlab("Number of cluster") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) ggsave(file.path("countsplit", "cluster_sse_nb.png"), width = 10, height = 5)Maybe you can check if the code is correct, especially concerning the sctransform part. Of note, I struggled again with columns and rows. As you can see in the code i had to 1) subset the matrix based on genes with known overdispersion and 2) tranpose the counts Matrix so that countsplit accepts it.
The resulting plot was similar to your simulation data in the way that it more clearly identified an optimal resolution. And the resolution was also lower (1.0) which is more plausible I think.
These are the results based on the
XestSubsetNBdata.-> I actually think that both are biologically plausible so it's exciting that the algorithm even suggests more clusters than the Seurat team chose.
Finally, I wanted to use the multifold negative binominal approach. Here's the code
Code Multifold Negative binominal
``` #negative binominal with multifold sctransform_fit <- sctransform::vst(pbmc.counts, residual_type = "none", min_cells = 10) overdisp_estimates <- sctransform_fit$model_pars_fit[,"theta"] #subset to genes with overdispersion estimates pbmcCountsSubset <- pbmc.counts[names(overdisp_estimates),] split3 <- countsplit(t(pbmcCountsSubset), folds = 10, overdisp = overdisp_estimates) split3_t <- lapply(split3, t) calcSS <- function(test_data, matrix, resRange, folds) { train_data <- matrix - test_data #create seurat object seu_obj_train <- CreateSeuratObject(counts = train_data, min.cells = 3, min.features = 200) seu_obj_train[["percent.mt"]] <- PercentageFeatureSet(seu_obj_train, pattern = "^MT-") seu_obj_train <- subset(seu_obj_train, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5) # same genes and cells in the train set as in the test set rows <- rownames(seu_obj_train) cols <- colnames(seu_obj_train) train_data <- train_data[rows,cols] test_data <- test_data[rows,cols] # run seurat pipeline seu_obj_train <- seu_obj_train |> NormalizeData() |> FindVariableFeatures() |> ScaleData() |> RunPCA() |> FindNeighbors(dims = 1:10) |> RunUMAP(dims = 1:10) # assign names for each resolution resNames <- paste0("RNA_snn_res.", resRange) for (res in resRange) { seu_obj_train <- FindClusters(seu_obj_train, resolution = res) } # calculcate within-cluster sse for each resolution results <- list() for (k in resNames) { results[k] <- cluster_sse2( trainDat = train_data, testDat = test_data, clusters.est = seu_obj_train@meta.data[[k]], eps.train = 1 - 1/folds) } return(data.frame(results)) } results_multifold <- lapply(1:10, function(x) calcSS(test_data = split3_t[[x]], matrix = pbmcCountsSubset, resRange = seq(0.2, 3, by = 0.2), folds = 10)) results_df <- results_multifold |> bind_rows() |> mutate(fold = 1:10) |> mutate(fold = sprintf("fold_%02d", fold)) |> tidyr::pivot_longer(cols = -fold, names_to = "cluster", values_to = "sse") |> group_by(fold) |> mutate(sse_norm = (sse - min(sse)) / (max(sse) - min(sse))) |> ungroup() |> group_by(cluster) |> mutate(mean_sse = mean(sse_norm)) |> ungroup() results_df |> ggplot(aes(x = cluster, y = sse_norm, group = fold, color = fold)) + geom_point() + geom_line() + geom_line(aes(y = mean_sse), color = "black", size = 1) + theme_bw() + ylab("Within-cluster SSE") + xlab("Number of cluster") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) ggsave(file.path("countsplit", "cluster_sse_multifold.png"), width = 10, height = 5) ```The fat black line is the mean SSE of all 10 folds. This also suggested a resolution of 1.0 (probably 0.8 and 1.2 would also be fine)
@anna-neufeld I am eager to hear you feedback. Especially, if all implementations are correct, like the sctransform, the refactoring of the ss function and the transposition of the matrices.
Thanks!