app.R

suppressPackageStartupMessages({
  library(shiny)
  library(shinythemes)
  library(shinycssloaders)
  library(dplyr)
  library(glue)
  library(Matrix)
  library(ggplot2)
  library(cowplot)
  library(markdown)
  library(DT)
  library(limma)
  library(edgeR)
  library(tibble)
})

options(spinner.type = 6)
options(spinner.color = "#ababab")
theme_set(theme_cowplot())

# source data
# cluster/sample info per cell - tibble
meta_tbl <- readRDS("data.meta.rds")
# sample info
sample_info <- readRDS("sample.info.rds")
# sctransform-integrated @data expression values - sparse matrix
exp_mat <- readRDS("data.exp.rds")
all.genes <- sort(rownames(exp_mat))
# Differential expression results
de <- readRDS("data.de.rds")
# Normalized expression from limma
v <- readRDS("data.v.rds")
# AUC TF regulon activity scores from SCENIC
auc <- readRDS("data.auc.rds")
# Differential regulon activity among genotypes (Wilcoxon ranked sum tests)
dra <- readRDS("data.dra.rds")
dra.genes <- sort(unique(unlist(lapply(dra, rownames))))

# cell type color palettes
#cell.names <- levels(meta_tbl$Cell_Type)
#colors_cells <- rainbow(length(cell.names))
#colors_cells <- setNames(colors_cells, cell.names)

colors_cells <- readRDS("colors.cells.rds")
colors_groups <- c(WT.F="#F8766D", WT.M="#00BFC4", SOD1G93A.F="red2", SOD1G93A.M="royalblue")
colors_genotypes <- c(SOD1G93A="red2", WT="royalblue")

# check that the cluster order is the same in the meta data table and the colors vector
#clusters_ordered <- meta_tbl %>% pull(Cell_Type) %>% levels()
#if (!identical(clusters_ordered, names(colors_cells))) stop("unexpected color order")

# ui: define UI for dataset viewer app
ui <- fluidPage(

  # layout: header ----
  theme = shinytheme("paper"),

  # layout: title ----
  titlePanel(div(HTML("Single-nucleus RNA-seq of lumbar spinal cord from amyotrophic lateral sclerosis <em>SOD1<sup>G93A</sup></em> model and wild-type mice"))),

  # layout: line break ----
  hr(),

  # layout: main (data) row of content ----
  tabsetPanel(
    tabPanel(
      title = "UMAP plots", fluid=TRUE,

      sidebarLayout(
        sidebarPanel(
          selectInput(
          inputId = "in_gene1",
          label = "Gene:",
          choices = c(Choose = "", all.genes),
          selected = "Gfap",
          selectize = TRUE
          ),

          radioButtons(
          inputId = "in_split_gen1",
          label = "Split plots by genotype?",
          choices = c(
            "Yes" = "Yes",
            "No" = "No"),
          selected = "Yes"
          ),

          radioButtons(
          inputId = "in_split_sex1",
          label = "Split plots by sex?",
          choices = c(
            "Yes" = "Yes",
            "No" = "No"),
          selected = "No"
          ),

          width=2
        ),

        mainPanel(
          fluidRow(
            column(
              width = 5,
              withSpinner(plotOutput("umap_gene_plot"))
            ),

            column(
              width = 5,
              withSpinner(plotOutput("umap_cluster_plot"))
            )
          ),

          fluidRow(
            column(
              width = 10,
              withSpinner(plotOutput("umap_gene_group_plot", height="800px"))
            )
          ),

          width=10
        ),

        position = "left",
        fluid = TRUE
      )
    ),

    tabPanel(
      title = "Violin plot", fluid=TRUE,

      sidebarLayout(
        sidebarPanel(
          selectInput(
            inputId = "in_gene2",
            label = "Gene:",
            choices = c(Choose = "", all.genes),
            selected = "Gfap",
            selectize = TRUE
          ),

          radioButtons(
            inputId = "in_split_gen2",
            label = "Split plots by genotype?",
            choices = c(
              "Yes" = "Yes",
              "No" = "No"),
            selected = "Yes"
          ),

          radioButtons(
            inputId = "in_split_sex2",
            label = "Split plots by sex?",
            choices = c(
              "Yes" = "Yes",
              "No" = "No"),
            selected = "No"
          ),

          width=2
        ),

        mainPanel(
          fluidRow(
            column(
              width = 10,
              withSpinner(plotOutput("vln_plot", height = "700px"))
            )
          ),

          width=10
        ),

        position = "left",
        fluid = TRUE
      )
    ),

    tabPanel(
      title = "Differential expression", fluid=TRUE,

      sidebarLayout(
        sidebarPanel(
          radioButtons(
            inputId = "in_type",
            label = "Query type:",
            choices = c(
              "Search gene in all results" = "gene",
              "Browse by cell type" = "cell"),
            selected = "gene"
          ),

          conditionalPanel(
            condition = "input.in_type == 'gene'",
            selectInput(
              inputId = "in_gene3",
              label = "Gene:",
              choices = c(Choose = "", all.genes),
              selected = "Gfap",
              selectize = TRUE
            )
          ),

            # Only show this panel if cell is selected
          conditionalPanel(
            condition = "input.in_type == 'cell'",
            selectInput(
              inputId = "in_cell",
              label = "Cell type:",
              choices = c(Choose = "", levels(meta_tbl$Cell_Type)),
              selected = "Astrocyte",
              selectize = TRUE
            ),

            selectInput(
              inputId = "in_test",
              label = "Comparison:",
              choices = c(Choose = "", names(de$Astrocyte)),
              selected = "SOD1G93A_v_WT",
              selectize = TRUE
            )
          ),

          width=2
        ),

        mainPanel(
          br(),
          conditionalPanel(
            condition = "input.in_type == 'cell'",
            fluidRow(
              column(
                width = 5,
                withSpinner(plotOutput("ma_plot") #, height = "500px"
                )
              ),

              column(
                width = 5,
                withSpinner(plotOutput("de_plot2") #, height = "500px"
                )
              )
            ),

            br(),

            fluidRow(
              column(
                width = 10,
                withSpinner(DT::dataTableOutput("cell_table"))
              )
            ),
          ),

          conditionalPanel(
            condition = "input.in_type == 'gene'",
            withSpinner(DT::dataTableOutput("gene_table"))
          ),

          width=10
        ),

        position = "left",
        fluid = TRUE
      )
    ),

    tabPanel(
      title = "SCENIC", fluid=TRUE,

      sidebarLayout(
        sidebarPanel(
          radioButtons(
            inputId = "in_type2",
            label = "Query type:",
            choices = c(
              "Search transcription factor in all results" = "tf",
              "Browse by cell type" = "cell2"),
            selected = "tf"
          ),

          conditionalPanel(
            condition = "input.in_type2 == 'tf'",
            selectInput(
              inputId = "in_gene4",
              label = "Gene:",
              choices = c(Choose = "", dra.genes),
              selected = "Fos",
              selectize = TRUE
            )
          ),

          # Only show this panel if cell is selected
          conditionalPanel(
            condition = "input.in_type2 == 'cell2'",
            selectInput(
              inputId = "in_cell2",
              label = "Cell type:",
              choices = c(Choose = "", names(dra)),
              selected = "Astrocyte",
              selectize = TRUE
            )
          ),

          width=2
        ),

        mainPanel(
          br(),
          conditionalPanel(
            condition = "input.in_type2 == 'cell2'",

            fluidRow(
              column(
                width = 10,
                withSpinner(plotOutput("dra_plot")
                )
              )#,
              #
              #column(
              #  width = 5,
              #  withSpinner(plotOutput("density_plot_f")
              #  )
              #)
            ),

            fluidRow(
              column(
                width = 10,
                withSpinner(DT::dataTableOutput("scenic_cell_table"))
              )
            )
          ),

          conditionalPanel(
            condition = "input.in_type2 == 'tf'",
            withSpinner(DT::dataTableOutput("scenic_gene_table"))
          ),

          width=10
        ),

        position = "left",
        fluid = TRUE
      )
    ),

    tabPanel(
      title = "About",
      fluidRow(
        column(width = 2),
        column(width = 8,
          includeMarkdown("about.md")),
        column(width = 2)
      )
    )
  )
)

# server: define server logic to summarize and view selected dataset ----
server <- function(input, output, session) {

  # Sync radio buttons and gene across tabs - must be better way to do this
  observe({
    m <- input$in_gene1
    updateSelectInput(session, "in_gene2", selected = m)
  })

  observe({
    n <- input$in_gene1
    updateSelectInput(session, "in_gene3", selected = n)
  })

  observe({
    o <- input$in_gene2
    updateSelectInput(session, "in_gene1", selected = o)
  })

  observe({
    p <- input$in_gene3
    updateSelectInput(session, "in_gene1", selected = p)
  })

  observe({
    q <- input$in_gene3
    updateSelectInput(session, "in_gene2", selected = q)
  })

  observe({
    r <- input$in_gene2
    updateSelectInput(session, "in_gene3", selected = r)
  })

  observe({
    a <- input$in_split_gen1
    updateRadioButtons(session, "in_split_gen2", selected = a)
  })

  observe({
    b <- input$in_split_gen2
    updateRadioButtons(session, "in_split_gen1", selected = b)
  })

  observe({
    y <- input$in_split_sex1
    updateRadioButtons(session, "in_split_sex2", selected = y)
  })

  observe({
    z <- input$in_split_sex2
    updateRadioButtons(session, "in_split_sex1", selected = z)
  })

  # generate single gene expression values table ----
  exp_tbl <- reactive({
    req(input$in_gene1) #, input$in_split_gen, input$in_split_gen
    set.seed(99)
    # adjust the meta data table based on the requested cells subset
    #meta_tbl <- meta_tbl %>% filter((Genotype == input$in_genotype | input$in_genotype == "Both") & (Sex == input$in_sex | input$in_sex == "Both"))

    tibble(Nucleus = colnames(exp_mat), exp_log = exp_mat[input$in_gene1, ]) %>%
      inner_join(meta_tbl, by = "Nucleus") %>%
      sample_frac()
  })

  # generate UMAP plots ----
  umap_plot_reactive <- reactive({
    req(input$in_gene1, input$in_split_gen1, input$in_split_sex1)

    umap_gene_plot <-
      ggplot(exp_tbl(), aes(x = UMAP_1, y = UMAP_2)) +
      geom_point(aes(color = exp_log), size = 0.08) +
      labs(title = paste(input$in_gene1, "expression in all nuclei")) +
      guides(color = guide_colorbar(title = "Normalized expression")) +
      theme(
        plot.title = element_text(hjust = 0.5),
        aspect.ratio = 1,
        axis.ticks = element_blank(),
        axis.text = element_blank()
      ) +
      scale_color_gradientn(colors = c("gray90", "darkblue"))

    umap_cluster_plot <-
      ggplot(exp_tbl(), aes(x = UMAP_1, y = UMAP_2)) +
      geom_point(aes(color = Cell_Type), size = 0.08, show.legend = TRUE) +
      labs(title = "Cell type") +
      theme(
        plot.title = element_text(hjust = 0.5),
        aspect.ratio = 1,
        axis.ticks = element_blank(),
        axis.text = element_blank()
      ) +
      guides(color = guide_legend(title = "Cell type", override.aes = list(size = 5))) +
      scale_color_manual(values = colors_cells)

    # Split UMAP plots if selected
    if (input$in_split_gen1 == "Yes" & input$in_split_sex1 == "Yes") {
      umap_gene_group_plot <-
        umap_gene_plot +
        labs(title = paste(input$in_gene1, "expression by genotype and sex")) +
        facet_wrap(vars(Group), scales = "fixed", nrow = 2, strip.position = "bottom") +
        theme(axis.text.x = element_blank())
    }  else if (input$in_split_gen1 == "Yes" & input$in_split_sex1 == "No"){
      umap_gene_group_plot <-
        umap_gene_plot +
        labs(title = paste(input$in_gene1, "expression by genotype")) +
        facet_wrap(vars(Genotype), scales = "fixed", nrow = 1, strip.position = "bottom") +
        theme(axis.text.x = element_blank())
    }  else if (input$in_split_gen1 == "No" & input$in_split_sex1 == "Yes"){
      umap_gene_group_plot <-
        umap_gene_plot +
        labs(title = paste(input$in_gene1, "expression by sex")) +
        facet_wrap(vars(Sex), scales = "fixed", nrow = 1, strip.position = "bottom") +
        theme(axis.text.x = element_blank())
    }  else if (input$in_split_gen1 == "No" & input$in_split_sex1 == "No") {
      umap_gene_group_plot <- ggplot() + theme_void()
    }

    align_plots(umap_gene_plot, umap_cluster_plot, umap_gene_group_plot, align = "hv", axis = "tblr")

  })

  # generate violin plot ----
  vln_plot_reactive <- reactive({
    req(input$in_gene2, input$in_split_gen2, input$in_split_sex2) #, input$in_genotype, input$in_sex,

    vln_plot <-
      ggplot(exp_tbl(), aes(x = Cell_Type, y = exp_log)) +
      geom_violin(aes(fill = Cell_Type, color = Cell_Type), scale = "width") +
      labs(title = input$in_gene2, x = "Cell type", y = "Normalized expression") +
      theme(
        plot.title = element_text(hjust = 0.5),
        axis.ticks.x = element_blank(),
        legend.position = "none",
        strip.background = element_blank(),
        strip.placement = "outside",
        strip.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        panel.spacing.x = unit(0.1, "lines")
      ) +
      scale_y_continuous(expand = c(0, 0)) +
      scale_fill_manual(values = colors_cells) +
      scale_color_manual(values = colors_cells)

    # split by genotype and/or sex if requested
    if (input$in_split_gen2 == "Yes" & input$in_split_sex2 == "Yes") {
      vln_plot <-
        vln_plot +
        facet_wrap(vars(Cell_Type, Group), scales = "free_x", nrow = 1, strip.position = "bottom") +
        theme(axis.text.x = element_blank())
    }  else if (input$in_split_gen2 == "Yes" & input$in_split_sex2 == "No"){
      vln_plot <-
        vln_plot +
        facet_wrap(vars(Cell_Type, Genotype), scales = "free_x", nrow = 1, strip.position = "bottom") +
        theme(axis.text.x = element_blank())
    }  else if (input$in_split_gen2 == "No" & input$in_split_sex2 == "Yes"){
      vln_plot <-
        vln_plot +
        facet_wrap(vars(Cell_Type, Sex), scales = "free_x", nrow = 1, strip.position = "bottom") +
        theme(axis.text.x = element_blank())
    }

    # return the plot
    vln_plot
  })

  gene_table_reactive <- reactive({
    req(input$in_gene3)
    temp <- lapply(unlist(de, recursive=F), function(x) x[input$in_gene3, c("logFC", "P.Value", "adj.P.Val")])
    temp.names <- names(temp)
    temp <- Reduce("rbind",temp)
    temp <- cbind(temp.names, temp)
    colnames(temp)[1] <- "CellType.Comparison"
    temp <- temp[complete.cases(temp), ]
    temp  # unsorted, all values
    #temp[temp$P.Value < 0.05]
  })

  cell_table_reactive <- reactive({
    req(input$in_test, input$in_cell)
    res <- de[[input$in_cell]][[input$in_test]]
    res
  })

  ma_plot_reactive = reactive({
    req(input$in_test, input$in_cell)

    res <- de[[input$in_cell]][[input$in_test]]

    col.vec <- ifelse(res$adj.P.Val < 0.1, "red2", "gray20")

    ma_plot <- ggplot(res, aes(x=AveExpr, y=logFC)) +
      geom_hline(yintercept=0, color="gray30") +
      geom_point(shape=16, alpha=0.7, aes(color = col.vec)) +
      labs(title = paste(input$in_test, "in", input$in_cell), x="Mean normalized expression", y="logFC") +
      scale_fill_manual(values = c("gray20", "red2")) +
      scale_color_manual(values = c("gray20", "red2")) +
      theme(plot.title = element_text(hjust = 0.5),
        legend.position = "none",
        strip.background = element_blank(),
        strip.placement = "outside",
        panel.spacing.x = unit(0.1, "lines"))

    if (length(input$cell_table_rows_selected))
      ma_plot <- ma_plot +
        geom_point(aes(x=res[input$cell_table_rows_selected, "AveExpr"], y=res[input$cell_table_rows_selected, "logFC"]), colour="blue", size=2, shape=16)
      #points(res[input$cell_table_rows_selected, c("AveExpr", "logFC"), drop = FALSE], cex = 1, col="blue")

    ma_plot
  })


  de_plot_reactive <- reactive({
    validate(
      need(input$cell_table_rows_selected != "", "Select a gene from table below")
    )

    req(input$in_test, input$in_cell, input$cell_table_rows_selected)

    v_cell <- v[[input$in_cell]]
    v_cell <- v_cell[rownames(de[[input$in_cell]][[input$in_test]]),]
    v_cell_gene <- v_cell[input$cell_table_rows_selected, ]

    gene_selected <- rownames(v_cell[input$cell_table_rows_selected, , drop=F])
    # adjust the meta data table based on the requested cells subset
    meta_sample <- sample_info[ifelse(as.character(unique(meta_tbl$SampleID)) %in% names(v_cell_gene), T, F), ]

    meta_sample_gene <- meta_sample %>% add_column(Expression=v_cell_gene)

    if (input$in_test == "SOD1G93A_v_WT_in_F") {
      meta_sample_gene <- meta_sample_gene %>% filter(Sex == "F")
    }
    else if (input$in_test == "SOD1G93A_v_WT_in_M"){
      meta_sample_gene <- meta_sample_gene %>% filter(Sex == "M")
    }
    else if (input$in_test == "F_v_M_in_SOD1G93A") {
      meta_sample_gene <- meta_sample_gene %>% filter(Genotype == "SOD1G93A")
    }
    else if (input$in_test == "F_v_M_in_WT"){
      meta_sample_gene <- meta_sample_gene %>% filter(Genotype == "WT")
    }


    de_plot <- ggplot(meta_sample_gene, aes(x=Group, y=Expression)) +
      geom_jitter(shape=16, width=0.15, alpha=1, aes(color=Group, size=1.5)) +
      stat_summary(fun=mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.3, color="black") +
      labs(title = paste(gene_selected, "in", input$in_cell), x = "Group", y = "Normalized expression") +
      theme(
        plot.title = element_text(hjust = 0.5),
        axis.ticks.x = element_blank(),
        legend.position = "none",
        strip.background = element_blank(),
        strip.placement = "outside",
        #strip.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        panel.spacing.x = unit(0.1, "lines")
      ) +
      #scale_y_continuous(expand = c(0, 0)) #+
      scale_fill_manual(values = colors_groups) +
      scale_color_manual(values = colors_groups)

    de_plot

  })

  scenic_gene_table_reactive <- reactive({
    req(input$in_gene4)
    temp2 <- lapply(dra, function(x) x[input$in_gene4, ])
    temp.names2 <- names(temp2)
    temp2 <- Reduce("rbind", temp2)
    temp2 <- cbind(temp.names2, temp2)
    colnames(temp2)[1] <- "Cell type"
    #temp2 <- temp2[complete.cases(temp2), ]
    temp2
  })

  scenic_cell_table_reactive <- reactive({
    req(input$in_cell2)
    res2 <- dra[[input$in_cell2]]
    res2
  })

  dra_plot_reactive <- reactive({
    validate(
      need(input$scenic_cell_table_rows_selected != "", "Select a transcription factor from table below")
    )

    req(input$in_cell2, input$scenic_cell_table_rows_selected)

    auc_cell <- auc[[input$in_cell2]]
    auc_cell <- auc_cell[ , rownames(dra[[input$in_cell2]])]
    auc_cell_gene <- auc_cell[ , input$scenic_cell_table_rows_selected, drop=F]

    gene_selected <- colnames(auc_cell_gene)
    # adjust the meta data table based on the requested cells subset
    meta_cell <- meta_tbl[ifelse(meta_tbl$Nucleus %in% rownames(auc_cell_gene), T, F), ]

    meta_cell_gene <- meta_cell %>% add_column(Score=auc_cell_gene[,1])


    dra_plot <- ggplot(data=meta_cell_gene, aes(x=Score, group=Genotype, fill=Genotype)) +
      labs(title = paste("SOD1G93A_v_WT in", input$in_cell2,  ":", gene_selected), x = "Transcription factor regulon activity score (AUC)", y = "Density") +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) +
      geom_density(alpha=0.5) +
      facet_wrap(vars(Sex), scales = "free", nrow = 1, strip.position = "bottom") #+
      #theme(axis.text.x = element_blank())

    dra_plot

  })

  # output: UMAP plot of expression levels ----
  output$umap_gene_plot <- renderPlot({
    # adding ggdraw due to cowplot::align_plots() output format
    ggdraw(umap_plot_reactive()[[1]])
  })

  # output: UMAP plot of clusters ----
  output$umap_cluster_plot <- renderPlot({
    # adding ggdraw due to cowplot::align_plots() output format
    ggdraw(umap_plot_reactive()[[2]])
  })

  output$umap_gene_group_plot <- renderPlot({
    # adding ggdraw due to cowplot::align_plots() output format
    ggdraw(umap_plot_reactive()[[3]])
  })

  # output: violin plot ----
  output$vln_plot <- renderPlot({
    vln_plot_reactive()
  })

  # output: gene de table ----
  output$gene_table <- DT::renderDataTable({
    datatable(gene_table_reactive(), options = list(pageLength = 70), rownames=F) %>%
      formatRound(2, digits = 2) %>% formatSignif(3:4, digits = 3)
  })

  # output: cell de table ----
  output$cell_table <- DT::renderDataTable({
    datatable(cell_table_reactive(), options = list(pageLength = 50), selection = 'single') %>%
      formatRound(c(3:5,8), digits = 2) %>% formatSignif(6:7, digits = 3)
  })

  # output: de jitter plot ----
  output$de_plot2 <- renderPlot({
    de_plot_reactive()
  })

  # output: MA plot ----
  output$ma_plot <- renderPlot({
    ma_plot_reactive()
  })

  # output: gene dra table ----
  output$scenic_gene_table <- DT::renderDataTable({
    datatable(scenic_gene_table_reactive(), options = list(pageLength = 10), rownames=F) %>%
      formatRound(c(2,4), digits = 2) %>% formatSignif(c(3,5), digits = 3)
  })

  # output: cell dra table ----
  output$scenic_cell_table <- DT::renderDataTable({
    datatable(scenic_cell_table_reactive(), options = list(pageLength = 50, order = list(2, 'asc')), selection = 'single') %>%
      formatRound(c(1,3), digits = 2) %>% formatSignif(c(2,4), digits = 3)
  })

  # output: dra jitter plot ----
  output$dra_plot <- renderPlot({
    dra_plot_reactive()
  })

}

# create Shiny app ----
shinyApp(ui = ui, server = server)