heatmapEnrichment(): conversion from wide to long format is incorrect

[toobiwankenobi]

Hi @ncborcherding

I think (not 100% sure though) the conversion from wide to long format performed in the heatmapEnrichment() function is incorrect. I think transposing the data.frame is causing a problem.

escape/R/heatmapEnrichment.R

Lines 97 to 102 in c260ec6

	long <- data.frame(
	variable = rep(gene.set, each = nrow(agg)),
	value = as.vector(t(agg[gene.set])),
	group = rep(agg[[group.by]], times = length(gene.set)),
	stringsAsFactors = FALSE
	)

Example: The line as.vector(t(...)) returns all values of the first row. However, both variable and group are constructed assuming values are returned row-by-row for the first column etc. This leads to misalignment.

I would suggest switching to tidyr::pivot_longer() to avoid this problem (unless you want to keep base R) and make the conversion less error-prone.

If you want to keep base R, I think this chunk would fix the issue

long <- data.frame(
  variable = rep(gene.set, each = nrow(agg)),
  value = unlist(agg[gene.set], use.names = FALSE),
  group = rep(agg[["ident"]], times = length(gene.set)),
  stringsAsFactors = FALSE
)

Could you validate this?

Let me know what you would like to do. I can create a pull request once you confirm this.

Best wishes,
Tobi

[ncborcherding]

Hey @toobiwankenobi

Nice catch in terms of code that could be improved. I just spent a bunch of time trying to reduce the dependencies so base R is the preference right now. Here the difference in defining value() are actually equivalent:

• t(agg[gene.set]) transposes the data, converting it from wide (genes as columns) to long (genes as rows), and then as.vector() flattens it row-wise from the transposed matrix.
• unlist(..., use.names = FALSE) flattens the data frame directly column-wise (which is equivalent to the behavior of as.vector(t(...))).

Using t() is more memory intensive and fragile, so your approach actually makes more sense. If you want to make a pull request, please feel free. Make sure you modify the group back to the original call or it will break group.by, it should look something like this:

long <- data.frame(
  variable = rep(gene.set, each = nrow(agg)),
  value = unlist(agg[gene.set], use.names = FALSE),
  group = rep(agg[[group.by]], times = length(gene.set)), 
  stringsAsFactors = FALSE
)

Thanks man

Nick

[toobiwankenobi]

Hi @ncborcherding

Thanks for your reply. I'm not sure if it is indeed equivalent.

gs <- list(Bcells = c("MS4A1", "CD79B", "CD79A", "IGH1", "IGH2"),
           Tcells = c("CD3E", "CD3D", "CD3G", "CD7","CD8A"))

pbmc <- SeuratObject::pbmc_small |>
  runEscape(gene.sets = gs, min.size = NULL, method = "AUCell")

# assign 3-letter grouping
pbmc$group_new <- rep(LETTERS[1:3], ncol(pbmc)/3, replace = TRUE)

Idents(pbmc) <- "group_new"
heatmapEnrichment(pbmc, assay = "escape", scale = T, palette = "viridis")

Try to run this code several times - for me, the resulting heatmap looks different every time I run the code.

Try 1

Try 2

Admittedly, changing the code to what I suggested also doesn't solve the problem. Unless it's a coloring issue of the values, I think the conversion from wide to long is wrong.

Best,
Tobi

[ncborcherding]

@toobiwankenobi

This is because AUCell is not deterministic and will result in a new set of values each time it runs - its a very frustrating thing and AUCell is so poorly maintained. You should get the same results for the other methods.

Nick

[toobiwankenobi]

Hi @ncborcherding

Good to know, I will stop using AUCell then. I now used fragments of the source code of heatmapEnrich(). Below you can find the reproducible example. The conversion from wide to long does not work properly.

gs <- list(
  Bcells = c("MS4A1", "CD79B", "CD79A", "IGH1", "IGH2"),
  Tcells = c("CD3E", "CD3D", "CD3G", "CD7", "CD8A"),
  CD4Tcells = c("CD3E", "CD4", "IL7R", "CCR7")
)

pbmc <- SeuratObject::pbmc_small |>
  runEscape(gene.sets = gs, min.size = NULL, method = "ssGSEA")

# assign 3-letter grouping
pbmc$group_new <- rep(LETTERS[1:3], ncol(pbmc)/3, replace = TRUE)


df <- escape:::.prepData(pbmc, 
                     assay = "escape", "all",
                     group.by = "group_new",
                     split.by = NULL,
                     facet.by = NULL)
# aggregate data
gene.set <- names(gs)
grp_cols <- "group_new"
agg <- aggregate(df[gene.set],
                 by   = df[grp_cols],
                 FUN  = base::mean,
                 SIMPLIFY = FALSE)

agg
  group_new   Bcells   Tcells CD4Tcells
1         A 64.39774 37.88984  41.93289
2         B 69.40864 36.48005  47.07375
3         C 65.29495 35.41495  39.19314

long <- data.frame(
  variable = rep(gene.set, each = nrow(agg)),
  value    = as.vector(t(agg[gene.set])),
  group    = rep(agg[["group_new"]], times = length(gene.set)),
  stringsAsFactors = FALSE
)

long
   variable    value group
1    Bcells 64.39774     A
2    Bcells 37.88984     B
3    Bcells 41.93289     C
...

The output of the updated code looks correct:

long <- data.frame(
  variable = rep(gene.set, each = nrow(agg)),
  value    = unlist(agg[gene.set], use.names = FALSE),
  group    = rep(agg[["group_new"]], times = length(gene.set)),
  stringsAsFactors = FALSE
)

long
   variable    value group
1    Bcells 64.39774     A
2    Bcells 69.40864     B
3    Bcells 65.29495     C
...

What do you think?

Best, Tobi

[ncborcherding]

@toobiwankenobi

You are absolutely correct!!!! Thank you for the amazing rundown - if you put in a pull request, I will get it fixed ASAP and pushed to bioconductor!

Great catch!

Nick

[toobiwankenobi]

Hi @ncborcherding

Done.

Have a great day, Tobi

[ncborcherding]

@toobiwankenobi Thanks again, got it pushed to the devel branch in bioconductor.

To clarify for any user finding this in the future, I did go an check, the current biocondcutor 3.21 (v2.4) escape does not have this issue. It is something that I introduced with an overhaul for dependencies in >= 2.5. It is now fixed in version 2.5.4.