In insightsengineering/hermes: Preprocessing, analyzing, and reporting of RNA-seq data

This vignette explores the calculation of correlation between samples and the resulting heatmap plot.

Idea

• One function cor() should calculate the correlations between samples.
• The other function draw_cor_heatmap() takes those results and creates the plot.
• The advantage of this modular design is that we can better test the correlation calculation and the user can also modify things / have their function before putting into the drawing function.

Annotation

Let's explore the annotation needed for the heatmap.

object <- hermes_data
anno <- data.frame(colData(object)[c("TechnicalFailureFlag", "LowDepthFlag")])
head(anno)
anno.char <- as.data.frame(lapply(anno, function(x) ifelse(is.na(x), "NA", as.character(x))))
rownames(anno.char) <- rownames(anno)
head(anno.char)

So effectively this just gets the two flag variables, converts them to character with explicit missings. We could maybe use the tern function for that:

tern::explicit_na(as.character(c(FALSE, NA)))

Correlation Calculation

Many users will know already the cor function, therefore the idea is to set that as a generic function and define a method for our hermes objects. We not only return the correlation matrix, but as an attribute flag_data also the two flags for the samples. This will be later used for the heatmap annotation.

setGeneric("cor")

.HermesDataCor <- setClass(
  Class = "HermesDataCor",
  contains = "matrix",
  slots = c(flag_data = "DataFrame")
)

setMethod(
  f = "cor",
  signature = signature(x = "AnyHermesData"),
  definition = function(x, y = "counts") {
    assert_that(is.string(y))

    chosen_assay <- assay(x, y)
    sample_cor_matrix <- stats::cor(chosen_assay, method = "pearson")

    .HermesDataCor(
      sample_cor_matrix,
      flag_data = colData(x)[, c("TechnicalFailureFlag", "LowDepthFlag")]
    )
  }
)

Note that we cannot use the argument name assay instead of y here since the original function does not take ... arguments.

Let's try it out:

result <- cor(x = object)
result[1:3, 1:3]
result@flag_data
class(result)

Note that we give the result an S4 class, just that we can check in the downstream plotting function that it comes from here. Actually... then we can just define a plot method to make it even more easy. We could also use S3 classes, but then it gets confusing in the package where we otherwise use S4 classes.

Heatmap Plotting

Fortunately we already have ComplexHeatmap in the docker. Also fortunately there is a translation from the previously used pheatmap arguments to this, see here.

library(ComplexHeatmap)

left_annotation <- rowAnnotation(
  LowDepthFlag = factor(anno$LowDepthFlag),
  col = list(LowDepthFlag = c("FALSE" = "green", "TRUE" = "red"))
)

top_annotation <- HeatmapAnnotation(
  TechnicalFailureFlag = factor(anno$TechnicalFailureFlag),
  col = list(TechnicalFailureFlag = c("FALSE" = "green", "TRUE" = "red"))
)

Heatmap(
  matrix = result,
  col = circlize::colorRamp2(c(0, 0.5, 1), c("red", "yellow", "green")),
  name = "Correlation",
  left_annotation = left_annotation,
  top_annotation = top_annotation
)

Now we can define the plot method for the correlation object:

setMethod(
  f = "plot",
  signature = "HermesDataCor",
  definition = function(x, 
                        flag_colors = c("FALSE" = "green", "TRUE" = "red"),
                        cor_colors = circlize::colorRamp2(c(0, 0.5, 1), c("red", "yellow", "green")),
                        ...) {
    df <- x@flag_data
    left_annotation <- rowAnnotation(
      LowDepthFlag = factor(df$LowDepthFlag),
      col = list(LowDepthFlag = flag_colors)
    )
    top_annotation <- HeatmapAnnotation(
      TechnicalFailureFlag = factor(df$TechnicalFailureFlag),
      col = list(TechnicalFailureFlag = flag_colors)
    )
    Heatmap(
      matrix = x,
      col = cor_colors,
      name = "Correlation",
      left_annotation = left_annotation,
      top_annotation = top_annotation,
      ...
    )
  }
)

Let's try it out:

plot(result)

# We can also pass options to `Heatmap()`:
plot(result, show_column_names = FALSE, show_row_names = FALSE)

insightsengineering/hermes documentation built on March 11, 2024, 11:04 p.m.