From 96eaf100d21af43adf2a7ba5dc37902e0a7939b4 Mon Sep 17 00:00:00 2001
From: Jacob Wagner <jpwagner@fredhutch.org>
Date: Wed, 22 Jul 2020 13:25:44 -0700
Subject: [PATCH] #34: Handle some edge cases in warpSet

---
 R/landmarkMatrix.R | 11 +++++++++++
 R/warpSet.R        | 14 ++++++++++++--
 2 files changed, 23 insertions(+), 2 deletions(-)

diff --git a/R/landmarkMatrix.R b/R/landmarkMatrix.R
index 6ed4cd2..3ade956 100644
--- a/R/landmarkMatrix.R
+++ b/R/landmarkMatrix.R
@@ -104,6 +104,17 @@ landmarkMatrix <- function(data, fres, parm, border=0.05, peakNr=NULL, densities
         matD[i,cl$cluster] <- cl$dist/diff(ranges[[1]][,parm])
     }
     resRegions <- resRegions[rownames(mat)]
+    
+    # TODO: Somehow it is possible for a landmark column to be NA for all samples
+    # and this is seed-dependent. This is a deeper bug, but for now remove these
+    # problematic landmarks so they don't cause obscure errors later
+    valid_landmarks <- apply(mat, 2, function(col) !all(is.na(col)))
+    mat <- mat[,valid_landmarks]
+    regRegions <- lapply(resRegions, function(sample_regions){
+      sample_regions[valid_landmarks,]
+    })
+    
+    matD <- matD[,valid_landmarks]
     attr(mat, "regions") <- resRegions
     attr(mat, "cdists") <- matD
     return(mat)
diff --git a/R/warpSet.R b/R/warpSet.R
index 84d8e18..9d9a4c2 100644
--- a/R/warpSet.R
+++ b/R/warpSet.R
@@ -704,8 +704,18 @@ warpSet.flowSet <- function(x, stains, grouping=NULL, monwrd=TRUE, subsample=NUL
 				stop("'", grouping, "' is not a phenoData variable.")
 			grps <- as.factor(pData(x)[,grouping])
 			anv <- numeric(ncol(landmarks))
-			for(i in seq_len(ncol(landmarks)))
-				anv[i] <- anova(lm(landmarks[,i] ~ grps))$Pr[1]
+			
+			for(i in seq_len(ncol(landmarks))){
+			  grps_present <- unique(grps[!is.na(landmarks[,i])])
+			  if(length(grps_present == 1)){
+			    warning(paste0("The following landmark is only present in a single group --",
+			                   "\nstain: ", p, 
+			                   "\nmean value: ", mean(landmarks[!is.na(landmarks[,i]),i])))
+			    anv[i] = 1.0 # pass the variance check below to avoid double warning
+			  }else{
+			    anv[i] <- anova(lm(landmarks[,i] ~ grps))$Pr[1] 
+			  }
+			}
 			if(any(anv < sig))
 				warning("Within-group variances are smaller than ",
 						"across-group variances for stain ", p,