fixed bugs due to empty samples, closes #13

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: premessa
 Type: Package
 Title: R package for pre-processing of flow and mass cytometry data
-Version: 0.1.6
+Version: 0.1.7
 Author: "Pier Federico Gherardini <pfgherardini@parkerici.org> [aut, cre]"
 Description: This package includes panel editing/renaming for FCS files, bead-based normalization and debarcoding.
 Imports: shiny (>= 0.14), flowCore, reshape, ggplot2, hexbin, gridExtra, rhandsontable, jsonlite

R/debarcode_cytof.R

@@ -147,7 +147,7 @@ get_sample_idx <- function(label, bc.results, sep.threshold, mahal.threshold = N
 #'  and the centroid of the debarcoded sample, as assigned using this threshold.
 #'
 #' @return Returns a vector of squared Mahalanobis distances, as calculated by \code{stats::mahalanobis}. The distance
-#'  values are capped at 30
+#'  values are capped at 30. If the distance calculation fails the corresponding vector values will be \code{NA}
 #'
 #' @export
 get_mahalanobis_distance <- function(m, bc.res, sep.threshold) {
@@ -158,9 +158,19 @@ get_mahalanobis_distance <- function(m, bc.res, sep.threshold) {
     for(lab in unique(assignments)) {
         sel.rows <- assignments == lab
         temp <- m[sel.rows, bc.channels]
-        cov.m <- cov(temp)
+
         #Here we are not taking the sqrt (different from normalizer)
-        mahal <- mahalanobis(temp, colMeans(temp), cov.m)
+        mahal <- tryCatch({
+                cov.m <- cov(temp)
+                mahalanobis(temp, colMeans(temp), cov.m, tol = 1e-30)
+            },
+            error = function(e) {
+                message(sprintf("Mahalanobis distance calculation failed for sample %s", lab))
+                message("Mahalanobis distance filtering will not be performed for this sample")
+                flush.console()
+                return(NA)
+            })
+
         ret[sel.rows] <- mahal
     }
     ret[ret > 30] <- 30
@@ -271,10 +281,12 @@ debarcode_fcs <- function(fcs, bc.key, output.dir, output.basename, sep.threshol
     assignments <- get_assignments_at_threshold(bc.res, sep.threshold, mahal.dist.threshold, mahal.dist)
 
     for(lab in unique(assignments)) {
-        temp <- m[assignments == lab,]
+        temp <- m[assignments == lab,, drop = F]
+
         out.fcs <- as_flowFrame(temp, fcs)
         out.fname <- file.path(output.dir, sprintf("%s_%s.fcs", output.basename, lab))
         write_flowFrame(out.fcs, out.fname)
+
     }
 }
 

README.md

@@ -226,7 +226,7 @@ Upon launching the GUI you will have access to the following controls:
 - *Current barcode key*: the CSV file containing the barcode key. By default the standard Fluidigm Cell-ID 20-Plex Barcoding key is used. If you want to select a different one, press the *Select key* button. 
 - *Current FCS file*: the FCS that is currently being debarcoded. Press the *Select FCS* button to select the FCS file you want to debarcode. Upon selecting both the FCS file and the key, the preliminary debarcoding process will start immediately. After a few seconds a number of diagnostic plots will appear in the right portion of the window (see [below](#plot-types))
 - *Minimum separation*: the minimum seperation between the positive and negative barcode channels that an event needs to have in order to be assigned to a sample. Events where the separation is less than this threshold are left unassigned. This filtering is done after rescaling the intensity of the barcode channels, and therefore the threshold should be a number between 0 and 1
-- *Maximum Mahalanobis distance*: the maximum distance between a single cell event, and the centroid of the sample the event has been assigned to. Events with distance greather than the threshold are left unassigned. The distance is capped at 30, so the default value of this option does not apply a filter based on Mahalanobis distance.
+- *Maximum Mahalanobis distance*: the maximum distance between a single cell event, and the centroid of the sample the event has been assigned to. Events with distance greather than the threshold are left unassigned. The distance is capped at 30, so the default value of this option does not apply a filter based on Mahalanobis distance. Note that in some cases (for instance when there are very few events in a sample), the Mahalanobis distance calculation may fail. In this case a diagnostic message is printed in the R console, and filtering based on Mahalanobis distance is not performed for the corresponding sample
 - *Plot type*: selects the type of plot to be displayed. Please see [below](#plot-types) for a description of the plots. Depending on the plot type, a few additional controls may be displayed:
   - *Select sample*: select a specific sample for plotting. Sample names are taken from the barcode key
   - *Select x axis*: select the channel to be displayed on the x axis