07-marker-B-sets.Rmd

```{r parameters-and-defaults, include = FALSE}
module <- "scRNAseq"
section <- "marker_genes"
```

```{r parameter-merge, include = FALSE}
local_params <- module %>%
  options() %>%
  magrittr::extract2(module) %>%
  magrittr::extract2(section) %>%
  ReporteR.base::validate_params(parameters_and_defaults)
```

```{r scRNAseq-marker-B-sets-check, include = FALSE, echo = FALSE}
if (length(local_params$set_facets) > 0) {
  assertive.sets::is_subset(local_params$set_facets, colnames(SummarizedExperiment::colData(object_filtered)))
}
```

### Marker gene sets

In contrast to individual genes, sets of genes are sometimes more appropriate to assign cellular identities to (sub-) populations. Figure \@ref(fig:scRNAseq-marker-B-sets-figure) shows aggregated expression of sets of marker genes, across `r ifelse(length(local_params$set_facets) > 0, ReporteR.base::itemize(local_params$set_facets), "all")` cells.

```{r scRNAseq-marker-B-sets-processing, include = FALSE, echo = FALSE}
assertive.files::assert_all_are_readable_files(local_params$feature_sets)

feature_sets <- yaml::read_yaml(local_params$feature_sets)

collapse <- TRUE
if (length(local_params$set_facets) > 0) {
  collapse <- FALSE
}

feature_set_expression <- lapply(names(feature_sets), function(name) {
  s <- feature_sets[[name]]
  m <- match(s, rownames(object_filtered))
  
  # Number of matched features
  n <- sum(!is.na(m))
  
  name <- glue::glue("{name} (n={n})")
  
  if (collapse) {
    exprs <- Matrix::colMeans(SummarizedExperiment::assay(object_filtered, local_params$assay)[na.omit(m), ])
    set_df <- data.frame(feature_set = rep(name, length(exprs)),
                         facet = rep(name, length(exprs)),
                         level = rep(name, length(exprs)),
                         avg_exprs = exprs)
  } else {
    facet_dfs <- lapply(local_params$set_facets, function(f) {
      levels <- levels(SummarizedExperiment::colData(object_filtered)[, f])
      facet_expression <- lapply(levels, function(l) {
        i <- which(SummarizedExperiment::colData(object_filtered)[, f] == l)
        Matrix::colMeans(SummarizedExperiment::assay(object_filtered, local_params$assay)[na.omit(m), i])
      })
      
      unlisted_facet_expression <- unlist(facet_expression)
      data.frame(feature_set = rep(name, times = length(unlisted_facet_expression)),
                 facet = rep(f, times = length(unlisted_facet_expression)),
                 level = rep(levels, times = sapply(facet_expression, length)),
                 avg_exprs = unlisted_facet_expression)
    })
    
    set_df <- do.call("rbind", facet_dfs)
  }
  return(set_df)
})

marker_sets_data <- do.call("rbind", feature_set_expression)
```

```{r scRNAseq-marker-B-sets-figure-params, include = FALSE, echo = FALSE}
theme_marker_sets <- ggplot2::theme(plot.background = ggplot2::element_blank(),
                                 panel.grid.major = ggplot2::element_blank(),
                                 panel.grid.minor = ggplot2::element_blank(),
                                 panel.border = ggplot2::element_rect(color = "black", size = 0.8, fill=NA),
                                 panel.background = ggplot2::element_blank(),
                                 axis.line.x = ggplot2::element_blank(),
                                 axis.line.y = ggplot2::element_blank(),
                                 axis.text = ggplot2::element_text(size = 5),
                                 axis.title.y = ggplot2::element_text(size = 7, margin = ggplot2::margin(0, -2, 0, 0)),
                                 axis.title.x = ggplot2::element_blank(),
                                 legend.text = ggplot2::element_text(size = 5),
                                 legend.title = ggplot2::element_text(size = 5),
                                 plot.title = ggplot2::element_text(face="bold", color="black", size=5),
                                 legend.key.size =  grid::unit(2, "mm"),
                                 legend.margin = ggplot2::margin(b = 2),#grid::unit(-50, "mm"),
                                 #legend.position = c(0,1),
                                 legend.justification = c(0,0),
                                 legend.background = ggplot2::element_rect(fill = "white"),
                                 #legend.direction = "horizontal"),
                                 strip.background = ggplot2::element_blank(),
                                 strip.text = ggplot2::element_text(size = 7))


fig_height <- ReporteR.base::estimate_figure_height(
  height_in_panels = ifelse(collapse, 1, length(feature_sets)),
  panel_height_in_in = params$formatting_defaults$figures$panel_height_in,
  axis_space_in_in = params$formatting_defaults$figures$axis_space_in,
  mpf_row_space = as.numeric(grid::convertUnit(grid::unit(5, 'mm'), 'in')),
  max_height_in_in = params$formatting_defaults$figures$max_height_in)

caption <- glue::glue("Expression across marker gene sets")
```

```{r scRNAseq-marker-B-sets-figure, echo = FALSE, message=FALSE, warning=FALSE, fig.height = fig_height$global, fig.cap = caption}
figure_marker_sets <- ggplot2::ggplot(marker_sets_data, ggplot2::aes(x = level, y = avg_exprs, fill = level)) +
  ggplot2::geom_boxplot(outlier.alpha = 0.6) +
  ggplot2::ylab("Average marker gene expression") +
  ggplot2::guides(fill = FALSE) +
  theme_marker_sets

if (!collapse) {
  figure_marker_sets <- figure_marker_sets +
    ggplot2::facet_grid(ggplot2::vars(feature_set), ggplot2::vars(facet), scales = "free")
}

figure_marker_sets
```