@@ -24,7 +24,7 @@ Imports:

     utils

 Suggests: 

     knitr,

-    SingleCellClassR.data,

+    scRNAseq,

     testthat

 VignetteBuilder: knitr

 Depends: R (>= 4.0), Seurat, SingleCellExperiment, SummarizedExperiment

@@ -170,7 +170,8 @@ setMethod("train_classifier", c("train_obj" = "Seurat"),

   clf$trainingData <- clf$resample <- clf$resampledCM <- NULL 

   p_thres <- 0.5

-  object <- SingleCellClassR(cell_type, clf, features, p_thres, NA_character_)

+  object <- SingleCellClassR(cell_type, clf, labels(clf$terms), p_thres, 

+                             NA_character_)

   # only assign parent if pretrained model for parent cell type is avai

   if ((!is.null(parent_process$parent.clf) 

@@ -266,7 +267,7 @@ setMethod("train_classifier", c("train_obj" = "SingleCellExperiment"),

   clf <- train_func(balance_ds$mat, balance_ds$tag)

   p_thres <- 0.5

-  object <- SingleCellClassR(cell_type, clf, features, p_thres, 

+  object <- SingleCellClassR(cell_type, clf, labels(clf$terms), p_thres, 

                              NA_character_)

   # only assign parent if pretrained model for parent cell type is avai

@@ -665,7 +666,7 @@ setMethod("classify_cells", c("classify_obj" = "Seurat"),

   if (any(pred_cells != "")) {

     pred_cells <- gsub("/$", "", pred_cells)

-  

+    

     # ignore ambiguous results

     if (ignore_ambiguous_result == TRUE) {

       pred_cells <- unlist(lapply(pred_cells, 

@@ -728,12 +729,16 @@ setMethod("classify_cells", c("classify_obj" = "SingleCellExperiment"),

     prediction <- make_prediction(filtered_mat, classifier, pred_cells, ignore_ambiguous_result)

     pred <- prediction$pred

     pred_cells <- prediction$pred_cells

-    

     # add prediction to meta data: 2 cols: p, class 

     for (colname in colnames(pred)) {

-      classify_obj[[colname]] <- pred[, colname, drop = FALSE]

+      classify_obj[[colname]] <- unlist(lapply(colnames(classify_obj), 

+                                               function(x) 

+                                                 if (x %in% rownames(pred)) {pred[x, colname]} 

+                                               else {NA}))

+      #pred[, colname, drop = FALSE, ]

     }

   }

+  

   if (any(pred_cells != "")) {

     pred_cells <- gsub("/$", "", pred_cells)

@@ -26,7 +26,7 @@

 #' @description Pretrained classifier obtained by training and testing on 

 #' The Tabula Muris Consortium.

 #' @docType data

-#' @usage default_models_n

+#' @usage default_models_m

 #' @format a list of \code{\link{SingleCellClassR}} objects

 #' @author Vy Nguyen, November 2020

 #' @keywords datasets

@@ -615,17 +615,16 @@ make_prediction <- function(mat, classifier, pred_cells, ignore_ambiguous_result

       if(x[1] >= classifier@p_thres) {"yes"} else {"no"}))

   rownames(pred) <- rownames(mat)

-  # append a sumary to whole predicted cell type

+  # append a summary to whole predicted cell type

   pred_cells <- unlist(lapply(cells,

   function(i)

     if (i %in% rownames(pred) && pred[i, "class"] == "yes") {

       if (ignore_ambiguous_result == TRUE && !is.na(classifier@parent) &&

           gsub("/", "", pred_cells[i]) == classifier@parent)

-      { classifier@cell_type }

+      { paste0(classifier@cell_type, "/") }

       else { paste0(pred_cells[i], classifier@cell_type, "/") }

     }

     else { pred_cells[i] }))

-  

   names(pred_cells) <- cells

   # remove no column and rename yes column to p

@@ -706,12 +705,13 @@ verify_parent <- function(mat, classifier, meta.data) {

   # parent clf, if avai, always has to be applied before children clf.

   parent_slot <- paste0(c(unlist(strsplit(classifier@parent, split = " ")), "class"), collapse = "_")

   if (parent_slot %in% colnames(meta.data)) {

-    parent_pred <- meta.data[, parent_slot, drop = FALSE]

-    pos_parent <- rownames(parent_pred[parent_pred == 'yes', , drop=FALSE]) 

+    parent_pred <- meta.data[, parent_slot]

+    pos_parent <- colnames(mat)[parent_pred == 'yes'] 

   } else {

     warning(paste0('Parent classifier of ', classifier@cell_type, 'cannot be applied.\n 

                    Please list/save parent classifier before child(ren) classifier.\n

-                   Skip applying classification models for ', classifier@cell_type, ' and its parent cell type.\n'), call. = FALSE, immediate. = TRUE)

+                   Skip applying classification models for ', classifier@cell_type, 

+                   ' and its parent cell type.\n'), call. = FALSE, immediate. = TRUE)

   }

   if (!is.null(pos_parent)) {

@@ -1 +1 @@

-c03fea4600edc3605663fbcfda019841a4dc9df3

\ No newline at end of file

+17113b58d0b1fd5d2f5297d30e533276223aca7a

\ No newline at end of file

@@ -8,7 +8,7 @@

 a list of \code{\link{SingleCellClassR}} objects

 }

 \usage{

-default_models_n

+default_models_m

 }

 \description{

 Pretrained classifier obtained by training and testing on 

@@ -29,41 +29,67 @@ any other cell type.

 ## Preparing train object and test object

-The workflow starts from a couple of Seurat objects where cells have been 

+The workflow starts from a couple of objects where cells have been 

 assigned to be different cell types. To do this, users may have annotated 

-scRNA-seq data (by a FACS-sorting process, for example), create a Seurat 

-object based on the sequencing data and assign the predetermined cell types 

-as Seurat meta data. If the scRNA-seq data has not been annotated yet, 

-another possible approach is to follow the basic Seurat workflow until 

-assigning cell type identity to clusters.

+scRNA-seq data (by a FACS-sorting process, for example), create a Seurat/

+SingleCellExperiment (SCE) object based on the sequencing data and assign the 

+predetermined cell types as cell meta data. If the scRNA-seq data has not 

+been annotated yet, another possible approach is to follow the basic 

+workflow (Seurat, for example) until assigning cell type identity to clusters.

+

+In this vignette, we use the human lung dataset from Zilionis et al., 2019,

+which is available in the scRNAseq (2.4.0) library. The dataset is stored as a

+SCE object.

 To start the training workflow, we first load the neccessary libraries. 

 ```{r}

 library(SingleCellClassR)

-library(SingleCellClassR.data)

+library(scRNAseq)

 ```

-One Seurat object will be used as train object, while the other is the test 

-object. In this example, we used Sade-Feldman dataset to create the train 

-object.

+Load the dataset:

 ```{r}

-data("feldman_seurat")

-feldman_seurat

+zilionis <- ZilionisLungData()

 ```

-We load Jerby-Arnon dataset for the testing object.

+We cut this dataset into two parts, one for the training and the other for the testing.

 ```{r}

-data("jerby_seurat")

-jerby_seurat

+pivot = ncol(zilionis)%/%2

+train_set <- zilionis[, 1:pivot]

+test_set <- zilionis[, (1+pivot):ncol(zilionis)]

 ```

-In our example, the cell type meta data is indicated as the active 

-identification of the Seurat object (in both train object and test 

-object). If cell type is stored in another slot of object meta data, 

+

+In this dataset, the cell type meta data is stored in the *Most likely LM22 cell type*

+slot of the SingleCellExperiment object (in both train object and test object). 

+If cell type is stored in another slot of object meta data, 

 the slot/tag slot name must be then provided as a parameter in the 

 train and test method. 

 ```{r}

-head(Idents(feldman_seurat))

+unique(train_set$`Most likely LM22 cell type`)

+```

+```{r}

+unique(test_set$`Most likely LM22 cell type`)

+```

+

+We want to train a classifier for B cells and their phenotypes. Considering memory B cells,

+naive B cells and plasma cells as B cell phenotypes, we convert all those cells to a uniform 

+cell label, ie. B cells. All non B cells are converted into 'others'.

+

+```{r}

+# change cell label

+train_set$B_cell <- unlist(lapply(train_set$`Most likely LM22 cell type`,

+                                  function(x) if (is.na(x)) {'ambiguous'} else if (x %in% c('Plasma cells', 'B cells memory', 'B cells naive')) {'B cells'} else {'others'}))

+

+test_set$B_cell <- unlist(lapply(test_set$`Most likely LM22 cell type`,

+                                 function(x) if (is.na(x)) {'ambiguous'} else if (x %in% c('Plasma cells', 'B cells memory', 'B cells naive')) {'B cells'} else {'others'}))

+```

+

+We observe that there are cells marked NAs. Those can be understood as 1/different from all indicated cell types or 2/any unknown cell types. Here we consider the second case, ie. we don't know whether they are positive or negative to B cells. To avoid the affect of those NAs cells, we can assign them as 'ambiguous'. All cells tagged 'ambiguous' will be ignored by SingleCellClassR from training and testing. 

+

+We may want to check the number of cells in each category:

+```{r}

+table(train_set$B_cell)

 ```

 ## Defining set of features

@@ -83,13 +109,13 @@ selected_features_B <- c("CD19", "MS4A1", "SDC1", "CD79A", "CD79B",

 ## Train model

 When the model is being trained, three most important information must be 

-provided are: the Seurat object used for training, the set of applied features

+provided are: the Seurat/SCE object used for training, the set of applied features

 and the cell type defining the trained model. 

 Cell type corresponding to the trained model must exist among identities 

 assigned to cells in the trained Seurat object. Remember if cell types 

 are not indicated as active identification of the trained object, name 

-of the tag slot in object meta data must be provided to the tag_slot parameter. 

+of the tag slot in object meta data must be provided to the sce_tag_slot parameter. 

 When training on a imbalanced dataset, the trained model may bias toward the 

 majority group and ignore the presence of the minority group. To avoid this, 

@@ -99,8 +125,8 @@ from the majority group. To use the same set of cells while training multiple

 times for one model, users can use set.seed. 

 ```{r}

 set.seed(123)

-clf_B <- train_classifier(train_obj = feldman_seurat, 

-features = selected_features_B, cell_type = "B cells")

+clf_B <- train_classifier(train_obj = train_set, cell_type = "B cells", features = selected_features_B,

+                          sce_assay = 'counts', sce_tag_slot = 'B_cell')

 ```

 ```{r}

 clf_B

@@ -115,7 +141,8 @@ clf(clf_B)

 ## Test model

 ```{r}

-clf_B_test <- test_classifier(test_obj = jerby_seurat, classifier = clf_B)

+clf_B_test <- test_classifier(test_obj = test_set, classifier = clf_B, 

+                              sce_assay = 'counts', sce_tag_slot = 'B_cell')

 ```

 ### Interpreting test model result

@@ -138,7 +165,7 @@ With the same classification model, the sensitivity and the specification of cla

 clf_B_test$overall_roc

 ```

-In this example of B cell classifier, the current threshold is at 0.5. The sensitivity is 0.9932203, and the specificity is 0.9890493 (FPR = 0.010950643). The higher sensitivity (0.9943503) can be reached if we set the p_thres at 0.4. However, we will have lower specificity (FPR = 0.013966037), which means that we misclassify more stranger cells as B cells. In contradiction, we may not retrieve all actual B cells with higher p_thres (0.6, for example).

+In this example of B cell classifier, the current threshold is at 0.5. The higher sensitivity can be reached if we set the p_thres at 0.4. However, we will then have lower specificity, which means that we misclassify more stranger cells as B cells. In contradiction, we may not retrieve all actual B cells with higher p_thres (0.6, for example).

 There is of course a certain trade-off between the sensitivity and the specificity of the model. Depending on the need of the project or the user-own preference, a probability threshold giving higher sensitivity or higher specificity can be chosen. In our perspective, p_thres at 0.5 is a good choice for the current B cell model.

@@ -50,33 +50,74 @@ clf_B <- default_models[['B cells']]

 clf_B

 ```

 ## Preparing train object and test object

+

 Same as training for basic models, training for a child model also requires 

-a train (Seurat) object and a test (Seurat) object. All Seurat objects must 

+a train (Seurat/SCE) object and a test (Seurat/SCE) object. All objects must 

 have a slot in meta data indicating the type of cells. Tag slot indicating 

 parent cell type can also be provided. In this case, parent cell type tag 

 will further be tested for coherence with the provided parent classifier.

+

+Cell tagged as child cell type but incoherent to parent cell type will be 

+removed from training and testing for the child cell type classifier. 

+

+In this vignette, we use the human lung dataset from Zilionis et al., 2019,

+which is available in the scRNAseq (2.4.0) library. The dataset is stored as a

+SCE object.

+

+To start the training workflow, we first load the neccessary libraries. 

 ```{r}

-library(SingleCellClassR.data)

-data("feldman_seurat")

+library(scRNAseq)

 ```

+

+Load the dataset:

+

 ```{r}

-# view train data

-feldman_seurat

+zilionis <- ZilionisLungData()

 ```

+

+We cut this dataset into two parts, one for the training and the other for the testing.

 ```{r}

-# tag slot indicating parent cell type, B cells

-head(Idents(feldman_seurat))

+pivot = ncol(zilionis)%/%2

+train_set <- zilionis[, 1:pivot]

+test_set <- zilionis[, (1+pivot):ncol(zilionis)]

 ```

+

+In this dataset, the cell type meta data is stored in the *Most likely LM22 cell type*

+slot of the SCE object (in both train object and test object). 

+If cell type is stored in another slot of object meta data, 

+the slot/tag slot name must be then provided as a parameter in the 

+train and test method. 

+

 ```{r}

-# tag slot indicating cell subset, plasma cells

-head(feldman_seurat[['Plasma_cells']])

+table(train_set$`Most likely LM22 cell type`)

 ```

-A second object for testing:

 ```{r}

-data("jerby_seurat")

+table(test_set$`Most likely LM22 cell type`)

+```

+

+Unlike the example of the training basic model, we will remove all NAs cells in order to reduce computationalc complexity.

-# view test data

-jerby_seurat

+```{r}

+# remove NAs cells

+train_set <- train_set[, !is.na(train_set$`Most likely LM22 cell type`)]

+test_set <- test_set[, !is.na(test_set$`Most likely LM22 cell type`)]

+```

+

+```{r}

+# convert cell label: 

+# 1 - positive to plasma cells, 

+# 0 - negative to plasma cells

+train_set$plasma <- unlist(lapply(train_set$`Most likely LM22 cell type`,

+                                  function(x) if (x == 'Plasma cells') {1} else {0}))

+

+test_set$plasma <- unlist(lapply(test_set$`Most likely LM22 cell type`,

+                                 function(x) if (x == 'Plasma cells') {1} else {0}))

+```

+

+We may want to check the number of cells in each category:

+```{r}

+table(train_set$plasma)

+# 1: plasma cells, 0: not plasma cells

 ```

 ## Defining set of features

@@ -109,17 +150,17 @@ path.to.models = '.'*

 Train the child classifier:

 ```{r}

 set.seed(123)

-clf_plasma <- train_classifier(train_obj = feldman_seurat, 

+clf_plasma <- train_classifier(train_obj = train_set, 

 features = selected_features_plasma, cell_type = "Plasma cells", 

-seurat_tag_slot = 'Plasma_cells', parent_clf = clf_B)

+sce_assay = 'counts', sce_tag_slot = 'plasma', parent_clf = clf_B)

 ```

 If B cells classifier has not been loaded to current working space, 

 an equivalent training process should be:

 ```{r}

 set.seed(123)

-clf_plasma <- train_classifier(train_obj = feldman_seurat, 

+clf_plasma <- train_classifier(train_obj = train_set, 

 features = selected_features_plasma, cell_type = "Plasma cells", 

-seurat_tag_slot = 'Plasma_cells', parent_cell = 'B cells')

+sce_assay = 'counts', sce_tag_slot = 'plasma', parent_cell = 'B cells')

 ```

 ```{r}

 clf_plasma

@@ -132,35 +173,17 @@ clf(clf_plasma)

 Parent classifier must be also indicated in test method.

 ```{r}

-clf_plasma_test <- test_classifier(test_obj = jerby_seurat, 

-classifier = clf_plasma, seurat_tag_slot = 'Plasma_cells', parent_clf = clf_B)

+clf_plasma_test <- test_classifier(test_obj = test_set, 

+classifier = clf_plasma, sce_assay = 'counts', sce_tag_slot = 'plasma', 

+parent_clf = clf_B)

 ```

 ### Interpreting test model result

+The test result obtained from a child model can be interpreted in the same way

+as we do with the model for basic cell types. We can change the prediction 

+probability threshold according to the research project or personal preference 

+and plot a roc curve.

-Same as testing a basic classification models, testing a child classifier also returns an object, which is a list of: *test_tag*, *pred*, *acc*, *auc*, *overall_roc*

-  

-The *overall_roc* is a summary of True Positive Rate (sensitivity) and False Positive Rate (1 - specificity) obtained by the trained model according to different thresholds:

-

-```{r}

-clf_plasma_test$overall_roc

-```

-

-We see that with the same TPR (*= 1.0*), the FPR increases if *p_thres* increases to *0.7*. Therefore, this is the better result as compare to the result produced by *p_thres = 0.5*. We change the prediction probability threshold as:

-

-```{r}

-p_thres(clf_plasma) <- 0.7

-```

-

-Hence the prediction result was changed as:

-```{r}

-clf_plasma_test <- test_classifier(test_obj = jerby_seurat, 

-classifier = clf_plasma, seurat_tag_slot = 'Plasma_cells', parent_clf = clf_B)

-```

-

-Now the specificity increased from *0.963759909399774* to *0.994337485843715* and with the same number of actual plasma cells (*65*), the number of predicted cells decreased from *97* to *70*.

-

-The ROC curve remains the same with the same AUC score:

 ```{r}

 print(clf_plasma_test$auc)

 roc_curve <- plot_roc_curve(test_result = clf_plasma_test)

@@ -172,11 +195,13 @@ plot(roc_curve)

 In order to save child classifier, parent classifier must have existed in the 

 classifier database, either in the package default database or in user-defined 

 database.

+

 ```{r}

 # see list of available model in package

 data("default_models")

 names(default_models)

 ```

+

 In our package, default models include already models classifying plasma cells.

 Therefore, we will save this model to a new local database specified by the 

 *path.to.models* parameter. If you start with a fresh new local database, 

@@ -195,61 +220,39 @@ save_new_model(new_model = clf_plasma, path.to.models = getwd(),

 When we save the B cells' classifier and the plasma cells' classifier, a local database is newly created. We can use this new database to classify cells in a Seurat or SingleCellExperiment object.

-Let's try to classify Jerby-Arnon dataset:

+Let's try to classify cells in the test set:

 ```{r}

-classified_jerby <- classify_cells(classify_obj = jerby_seurat, 

-                                   cell_types = 'all', path_to_models = getwd())

+classified <- classify_cells(classify_obj = test_set, sce_assay = 'counts', 

+                             cell_types = 'all', path_to_models = getwd())

 ```

 Using the *classify_cells()* function, we have to indicate exactly the repository containing the database that the models has recently been saved to. In the previous section, we saved our new models to the current working directory. 

-In the *classified_jerby* object, the classification process added new columns to the cell meta data, including the *predicted_cell_type* and *most_probable_cell_type* columns. Let's take a look at the original plasma cells tag:

-

-```{r}

-# get a summary of the plasma cells tag

-table(classified_jerby[['Plasma_cells']][,1])

-# 1 corresponds to cells positive to B cells and 0 corresponds to cells negative to B cells

-```

-65 cells (labeled 1) are plasma cells and 7121 cells are not plasma cells.

+In the *classified* object, the classification process added new columns to the cell meta data, including the *predicted_cell_type* and *most_probable_cell_type* columns. 

 If we use the full prediction to compare with actual plasma tag, we obtain this result: 

 ```{r}

 # compare the prediction with actual cell tag

-table(classified_jerby[['predicted_cell_type']][,1], classified_jerby[['Plasma_cells']][,1])

+table(classified$predicted_cell_type, classified$plasma)

+# plasma cell is child cell type of B cell

+# so of course, all predicted plasma cells are predicted B cells 

 ```

-We find that 65 actual plasma cells were assigned as B cells and as plasma cells at the same time. This is a reasonable result because in this example, we consider all plasma cells are B cells.

-

-However, comparing the actual tag with the most probable prediction, we obtain: 

+When comparing the actual tag with the most probable prediction, we obtain: 

 ```{r}

 # compare the prediction with actual cell tag

-table(classified_jerby[['most_probable_cell_type']][,1], classified_jerby[['Plasma_cells']][,1])

+table(classified$most_probable_cell_type, classified$plasma)

 ```

-This is a quite surprise result. Among the 65 actual plasma cells, only 1 was identified but 64 others were identified as B cells. In contradiction, the testing process of *clf_plasma* proved that this model could detect all plasma cells in the same dataset. Why this can happen? This contradictory can be explained.

-

-The *predicted_cell_type* takes all predictions having the probabilities satisfying the corresponding probability thresholds. In this case, SingleCellClassR takes all predicted B cells having probability from 0.5 and all predicted plasma cells having probability from 0.7. Meanwhile, the *most_probable_cell_type* takes only the cell type which gives highest prediction probability. Here, 64 predicted plasma cells were also predicted as B cells. However, they have a higher probability to be B cells than to be plasma cells. Therefore, they were classified as B cells in the *most_probable_cell_type*.

+The number of identified plasma cells is different in the *predicted_cell_type* slot and in the *most_probable_cell_type*. This is because the *predicted_cell_type* takes all predictions having the probabilities satisfying the corresponding probability thresholds. Meanwhile, the *most_probable_cell_type* takes only the cell type which gives highest prediction probability. 

 To have all plasma cells specified as plasma cells, we can set the *ignore_ambiguous_result* to TRUE. This option will actually hide all ambiguous prediction in case we have more distinct cell types. In the parent-chid(ren) relationship of cell types, the more specified cell types/phenotypes will be reported. Of course, we don't obtain the *most_probable_cell_type* in cell meta data.

 ```{r}

-classified_jerby <- classify_cells(classify_obj = jerby_seurat, 

-                                   cell_types = 'all', path_to_models = getwd(),

-                                   ignore_ambiguous_result = TRUE)

-table(classified_jerby[['predicted_cell_type']][,1], classified_jerby[['Plasma_cells']][,1])

-```

-

-To check the B cell classification with actual B cells, we can do as follows:

-

-```{r}

-classified_jerby <- classify_cells(classify_obj = jerby_seurat, 

-                                   cell_types = 'all', path_to_models = getwd())

-# convert B cells label to binaries to ignore details of other cell types

-label_B <- Idents(classified_jerby)

-label_B <- unlist(lapply(label_B, function(x) if (x == 'B cells') {1} else 0))

-

-# compare prediction with cell label

-table(classified_jerby[['most_probable_cell_type']][,1], label_B)

+classified <- classify_cells(classify_obj = test_set, sce_assay = 'counts',

+                             cell_types = 'all', path_to_models = getwd(),

+                             ignore_ambiguous_result = TRUE)

+table(classified$predicted_cell_type, classified$plasma)

 ```

 ## Session Info