refactor plot heatmap from list

Previous bug: if 1 gene with known name but no expression plotted along two or more genes with expression (e.g. mpst-5,msp-10,ZK512.10), would fail.Largely rewrote the heatmap plotting code in the process for clarity.

server.R

@@ -537,113 +537,182 @@ server <- function(input, output) {
 
   #~ From list ----
   observeEvent(input$PlotHeatmapFromList, {
-    cat("--> heatmap PlotHeatmapFromList\n")
 
-    load_as_needed("med.scaled.long")
-    load_as_needed("L4.TPM.raw.scaled.long")
-
-    ds <- input$dataset_heatmap
-    #ss <- unlist(strsplit(as.character(input$genelist), split = ","))
-    #ss <- gsub(" ", "", as.character(ss))
-    ss <- strsplit(as.character(input$genelist), "\n| |\\,|\t")
-    ss <- as.data.frame(ss)[,1]
-    star <- grep("\\*", ss)
-
-    families <- c()
-    for(i in star){ 
-      families <- c(families, grep( gsub("\\*", "", ss[i]), gene_list$gene_name, value = TRUE) )
-    }
+    cat("--> heatmap PlotHeatmapFromList\n")
 
-    ss <- unique(c(families, ss, filter(gene_list, gene_id %in% ss | seqnames %in% ss)$gene_name))
 
 
-    mis <- ss[ss %in% c(gene_list$gene_id, gene_list$gene_name, gene_list$seqnames) &
-                !ss %in% med.scaled.long$gene_name & ss %in% gene_list$gene_name]
-    mis_all <- ss[ss %in% c(gene_list$gene_id, gene_list$gene_name, gene_list$seqnames) &
-                    !ss %in% L4.TPM.raw.scaled.long$gene_name & ss %in% gene_list$gene_name]
+    #~~ Load dataset ----
+    ds <- input$dataset_heatmap
 
-    if(ds=="Neurons only"){
+    if(ds == "Neurons only"){
 
-      load_as_needed("L4.TPM.medium")
+      load_as_needed("med.scaled.long")
 
-      L4.TPM=L4.TPM.medium
       heatmapdata=med.scaled.long
-      cc = colnames(ths)[-c(1,130,131)]
-      missing = mis
+
+
     } else {
 
-      load_as_needed("L4.all.TPM.raw")
+      load_as_needed("L4.TPM.raw.scaled.long")
 
-      L4.TPM=as(L4.all.TPM.raw,"dgCMatrix")
       heatmapdata=L4.TPM.raw.scaled.long
-      cc=colnames(L4.all.TPM.raw)
-      missing = mis_all
+
     }
 
 
-    flp.neuron.scaled <- heatmapdata[which(heatmapdata$gene_name %in% ss),]
-    flp.ids <- as.character(vlookup(unique(flp.neuron.scaled$gene_name), gene_list, result_column = 1, lookup_column = 2))
-    flp.expr <- L4.TPM[flp.ids,, drop=FALSE]
+    #~~ Process gene input list ----
+    ss <- strsplit(as.character(input$genelist), "\n| |\\,|\t")[[1]] |>
+      setdiff("") |>
+      unique()
 
-    # order
-    if ( nrow(flp.expr) > 1 && input$reorder_rows ) {
-      flp.neuron.order <- pheatmap(flp.expr, scale = "row")
-      flp.neuron.order <- flp.neuron.order[["tree_row"]]$order
-      flp.neuron.order <- flp.ids[flp.neuron.order]
+
+    # expand stars
+    input_has_star <- grepl("\\*", ss)
+
+    expanded_stars <- lapply(which(input_has_star),
+           \(.i){
+             pattern <- paste0("^", gsub("\\*", ".+", ss[[.i]]), "$")
+             grep(pattern, gene_list$gene_name, value = TRUE)
+           }) |>
+      unlist()
+
+    ss <- ss[!input_has_star]
+    ss <- union(ss, expanded_stars)
+
+
+    # convert to gene_names
+    is_a_gene_id <- !(ss %in% gene_list$gene_name) &
+      ss %in% gene_list$gene_id
+
+    is_a_seq_name <- !(ss %in% gene_list$gene_name) &
+      ss %in% gene_list$seqnames
+
+    ss[is_a_gene_id] <- vlookup(ss[is_a_gene_id], gene_list, result_column = 2, lookup_column = 1)
+    ss[is_a_seq_name] <- vlookup(ss[is_a_seq_name], gene_list, result_column = 2, lookup_column = 3)
+
+
+
+    input_name_unknown <- ss[ !(ss %in% gene_list$gene_name) ]
+    input_has_no_expression_data <- ss[!ss %in% heatmapdata$gene_name]
+
+    input_known_but_no_expression <- setdiff(input_has_no_expression_data, input_name_unknown)
+
+    ss_known <- ss |> setdiff(input_has_no_expression_data)
+
+    message("Heatmap from list, ", length(ss_known)," genes: ", paste(ss_known, collapse = ","))
+    if(length(input_has_no_expression_data) > 0){
+      message(length(input_name_unknown)," genes unknown: ",
+              input_name_unknown)
+      message(length(input_known_but_no_expression)," genes known but no data: ",
+              input_known_but_no_expression)
+    }
+
+    #~~ Reorder ----
+    if ( length(ss_known) <= 1 || !input$reorder_rows ){
+      ordered_ss_known <- ss_known
 
-    } else if(nrow(flp.expr) > 1 && ! input$reorder_rows){
-      flp.neuron.order <- as.character(vlookup(unique(ss), gene_list, result_column = 1, lookup_column = 2))
+    } else {
+
+      # get full expression matrix
+      if(ds=="Neurons only"){
+
+        load_as_needed("L4.TPM.medium")
+        expr_matrix <- L4.TPM.medium
+
+      } else {
+
+        load_as_needed("L4.all.TPM.raw")
+        expr_matrix <- as(L4.all.TPM.raw,"dgCMatrix")
+
+      }
+
+      ss_known_gene_id <- vlookup(ss_known, gene_list, result_column = 1, lookup_column = 2)
+
+      order_for_ss_ids <- expr_matrix[ss_known_gene_id, , drop = FALSE] |>
+        t() |>
+        scale() |>
+        t() |>
+        dist() |>
+        hclust() |>
+        purrr::pluck("order")
+
+      ordered_ss_known <- ss_known[order_for_ss_ids]
 
-    } else{
-      flp.neuron.order <- flp.ids[1]
     }
-    flp.neuron.order <- as.character(vlookup(flp.neuron.order, gene_list))
-    flp.neuron.scaled$gene_name <- factor(flp.neuron.scaled$gene_name, levels = c(rev(flp.neuron.order), missing))
-    #flp.neuron.scaled$gene_name <- fct_rev(flp.neuron.scaled$gene_name)
 
-    for( i in missing ){
-      dff <- data.frame(gene_name=i, cell.type= cc, scaled.expr=0, prop=0, Modality="NA")
-      if(ds!="Neurons only"){
-        colnames(dff)[5]<-"tissue"
-        dff$tissue <- filter(L4.TPM.raw.scaled.long, gene_name=="nduo-6")$tissue
-        }
-      flp.neuron.scaled <- rbind(flp.neuron.scaled, dff)
+    #~~ Prepare dataset for plotting ----
+    heatmapdata <- heatmapdata[heatmapdata$gene_name %in% ordered_ss_known, ]
+
+    if(length(input_known_but_no_expression) > 0){
+
+      fake_heatmap_data <- expand.grid(
+        gene_name=input_known_but_no_expression,
+        cell.type= cc,
+        scaled.expr=0,
+        prop=0,
+        Modality="NA"
+      )
+
+      if(ds != "Neurons only"){
+        tissues <- heatmapdata$tissue[heatmapdata$gene_name == "nduo-6"]
+        fake_heatmap_data$tissue <- tissues
+      }
+
+      heatmapdata <- rbind(heatmapdata, fake_heatmap_data)
     }
 
-    if (ds!="Neurons only"){
+    heatmapdata$gene_name <- factor(heatmapdata$gene_name,
+                                    levels = c(
+                                      rev(ordered_ss_known),
+                                      input_known_but_no_expression
+                                    ))
+
+
+
+
+    #~~ Plot ----
+
+    g_base <- ggplot(heatmapdata,
+                aes(y = gene_name, x = cell.type)) +
+      geom_point(aes(color = scaled.expr, size = prop)) +
+      scale_color_gradientn("Scaled TPM", colors = c("orange", "maroon", "navy")) +
+      scale_size_continuous(name = "Proportion", limit = c(0.5, 100), range = c(0,5)) +
+      theme(
+        panel.background = element_blank()
+      )
+
+
+
+    if (ds != "Neurons only"){
 
-      if(nrow(flp.expr) >1){
-        g <- ggplot(flp.neuron.scaled, aes(y = gene_name, x = cell.type)) + 
-          geom_point(aes(color = scaled.expr, size = prop)) +
-          theme(panel.background = element_blank(), axis.text.y.left = element_text(size = 12, color = "black"),
-                legend.key = element_blank(),
-                legend.text = element_text(color = "black", size = 14),
-                legend.title = element_text(color = "black", size = 16),
-                axis.title = element_text(size = 16, color = "black")) + 
-          scale_color_gradientn("Scaled TPM", colors = c("orange", "maroon", "navy")) + 
-          scale_size_continuous(name = "Proportion", limit = c(0.5, 100), range = c(0,5)) + 
-          labs(y = "Gene", x = "Tissue") + theme(panel.grid = element_line(size = 0.2, color = "grey85")) +
-          facet_grid(~tissue, scales = "free_x", space = "free_x", switch = "x") + 
-          theme(strip.placement = "outside", 
-                strip.background.x = element_blank(),
-                axis.text.x.bottom = element_blank(),
-                strip.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 12)) 
-      } else {
-
-        g <- ggplot(flp.neuron.scaled, aes(y = gene_name, x = cell.type)) + 
-          geom_point(aes(color = scaled.expr, size = prop)) +
-          theme(panel.background = element_blank(), axis.text.y.left = element_text(size = 12, color = "black"),
-                legend.key = element_blank(),
-                legend.text = element_text(color = "black", size = 14),
-                legend.title = element_text(color = "black", size = 16),
-                axis.title = element_text(size = 16, color = "black")) + 
-          scale_color_gradientn("Scaled TPM", colors = c("orange", "maroon", "navy")) + 
-          scale_size_continuous(name = "Proportion", limit = c(0.5, 100), range = c(0,5)) + 
-          labs(y = "Gene", x = "Tissue") + theme(panel.grid = element_line(size = 0.2, color = "grey85"))
+      g <- g_base +
+        labs(y = "Gene", x = "Tissue") +
+        theme(
+          axis.text.y.left = element_text(size = 12, color = "black"),
+          legend.key = element_blank(),
+          legend.text = element_text(size = 14, color = "black"),
+          legend.title = element_text(size = 16, color = "black"),
+          axis.title = element_text(size = 16, color = "black"),
+          # grid
+          panel.grid = element_line(linewidth = 0.2, color = "grey85")
+        )
+
+
+      if(length(ordered_ss_known) > 1){
+        g <- g +
+          facet_grid(~tissue, scales = "free_x", space = "free_x", switch = "x") +
+          theme(
+            strip.placement = "outside", 
+            strip.background.x = element_blank(),
+            axis.text.x.bottom = element_blank(),
+            strip.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 12)
+          ) 
       }
 
       pg <- ggplotGrob(g)
 
+
       for(i in which(grepl("strip-b", pg$layout$name))){
         pg$grobs[[i]]$layout$clip <- "off"
       }
@@ -653,16 +722,20 @@ server <- function(input, output) {
           withProgress(message = "Generating heatmap Plot...", value = 0, {
             pg      })
         }
+
+
     } else {
 
-      g<-ggplot(flp.neuron.scaled, aes(y =gene_name, x = cell.type)) +
-        geom_point(aes(color = scaled.expr, size = prop)) +
-        #geom_point(data = NULL, aes(y = vec), pch = NA) +
-        theme(axis.text.x.bottom = element_text(angle = 90, vjust= 0.5, size = 7, color = "black", hjust = 1),
-              panel.background = element_blank(), axis.text.y.left = element_text(size = 7, color = "black")) +
-        scale_color_gradientn("Scaled TPM", colors = c("orange", "maroon", "navy")) +
-        scale_size_continuous(name = "Proportion", limit = c(0.5, 100), range = c(0,5)) +
-        labs(y = "Gene", x= "Neuron") + theme(panel.grid = element_line(size = 0.5, color = "grey85"))
+      g <- g_base +
+        labs(y = "Gene", x= "Neuron") +
+        theme(
+          axis.text.x.bottom = element_text(angle = 90, vjust= 0.5,
+                                            size = 7, color = "black",
+                                            hjust = 1),
+          axis.text.y.left = element_text(size = 7, color = "black"),
+          panel.grid = element_line(linewidth = 0.5, color = "grey85")
+        )
+
 
       fnh <-
         function() {
@@ -718,7 +791,7 @@ server <- function(input, output) {
     output$vals <- renderPrint({
       hover <- input$plot_hover 
       #print(input$plot_hover) # list
-      y <- nearPoints(flp.neuron.scaled, input$plot_hover)
+      y <- nearPoints(heatmapdata, input$plot_hover)
       req(nrow(y) != 0)
       y
     })