--- a
+++ b/code_final/visutils.R
@@ -0,0 +1,100 @@
+loadC2L = function(d,sd2zero=0){
+  r = read.csv(paste0(d,'/q05_cell_abundance_w_sf.csv'),row.names = 1,check.names = FALSE)
+  sd = read.csv(paste0(d,'/stds_cell_abundance_w_sf.csv'),row.names = 1,check.names = FALSE)
+  colnames(r) = sub('q05cell_abundance_w_sf_','',colnames(r))
+  colnames(sd) = sub('stdscell_abundance_w_sf_','',colnames(sd))
+  sd = sd [rownames(r),colnames(r)]
+  r[r<sd*sd2zero] = 0
+  mm = do.call(rbind,strsplit(rownames(r),'|',T))
+  r$barcode = mm[,2]
+  r = split(r,mm[,1])
+
+  for(i in 1:length(r)){
+    rownames(r[[i]]) = r[[i]]$barcode
+    r[[i]]$barcode = NULL
+    r[[i]] = as.matrix(r[[i]])
+  }
+  r
+}
+plotNMFConsSummary <- function(
+      coefs, cons,
+      clcols = NULL,
+      max.cex = 4.5/14*8,
+      colfun = function(x) {
+        visutils::num2col(
+          x, c('blue','gray','orange','violet','black'), minx = 0, maxx = 1
+        )
+      }, ylab.cex = 1, xlab.cex = 1
+){
+  cons = cons[colnames(coefs),colnames(coefs)]
+  cls = apply(coefs,2,which.max)
+  if(is.null(clcols))
+    clcols = RColorBrewer::brewer.pal(nrow(coefs),'Set3')
+  o = names(cls)[order(cls)]
+  parl = par(no.readonly=TRUE)
+  par(mar=c(3,9,2.5,6),bty='n',oma=c(0,0,1,0),cex=1,tcl=-0.2,mgp=c(1.2,0.3,0),las=1)
+  dotPlot(
+    t(coefs[, o]),
+    rowColours = cbind(clcols[cls[o]]),
+    colColours = clcols,
+    max.cex = max.cex, ylab.cex = ylab.cex,
+    xlab.cex = xlab.cex,
+    scaleWM = F,
+    colfun = colfun,
+    plot.legend = TRUE
+  )
+  par(parl)
+}
+plotFeatureProfiles_fun <- function (
+  m, features, features_facet = NULL, cols = NULL, sd.mult = 2, legend. = TRUE,
+  ylim = NULL, scaleY = TRUE, area.opacity = 0.2, lwd = 2,
+  xlab = "Distance (spots)", ylab = "Relative abundance", main = "",
+  xlim = NULL, ...
+) {
+  if (is.null(features_facet))
+    features_facet = sapply(features, function(x) dimnames(m)[[3]], simplify = FALSE)
+  if (is.null(cols))
+      cols = char2col(features)
+  if (is.null(names(cols)))
+      names(cols) = features
+  x = as.numeric(dimnames(m)[[1]])
+  areas = lapply(features, function(ct) {
+      do.call(rbind, apply(
+        m[, gsub("Healthy ", "", ct), unlist(features_facet[ct])], 1,
+        function(x) {
+          x = x[!is.na(x)]
+          data.frame(mean = mean(x), sd = sd(x)/sqrt(length(x)),
+              n = length(x))
+      }))
+  })
+  names(areas) = features
+  if (scaleY)
+      areas = lapply(areas, function(a) {
+          a[, 1:2] = a[, 1:2]/max(a$mean, na.rm = TRUE)
+          a
+      })
+  if (is.null(ylim))
+      ylim = range(sapply(areas, function(a) c(min(a$mean -
+          a$sd * sd.mult, na.rm = TRUE), max(a$mean + a$sd *
+          sd.mult, na.rm = TRUE))))
+  if (is.null(xlim))
+      xlim = range(x)
+  plot(1, t = "n", xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab,
+      main = main, ...)
+  for (n in names(areas)) plotArea(x, areas[[n]][, 1:2], sd.mult = sd.mult,
+      col = cols[n], new = FALSE, lwd = lwd, area.transp = area.opacity)
+  if (!is.null(legend.) && (!is.logical(legend.) || legend.)) {
+      if (!is.list(legend.))
+          legend. = list()
+      if (is.null(legend.$x)) {
+          legend.$x = grconvertX(1, "npc", "user")
+          legend.$y = grconvertY(1, "npc", "user")
+      }
+      legend.$fill = cols
+      legend.$legend = names(cols)
+      legend.$bty = par("bty")
+      legend.$border = NA
+      legend.$xpd = NA
+      do.call(legend, legend.)
+  }
+}
\ No newline at end of file