agglomerateByPrevalence introduces wrong relative abundances

[artur-sannikov]

Describe the bug
After running agglomerateByPrevalence on a relative abundance assay, the total sum of relative abundances per sample is no longer 1.

To Reproduce

library(mia)
data("GlobalPatterns", package = "mia")
tse <- GlobalPatterns

tse <- transformAssay(tse, method = "relabundance")
print(sum(assay(tse, "relabundance")[, "CL3"])) # [1] 1

tse <- agglomerateByPrevalence(
  tse,
  assay.type = "relabundance",
  rank = "Genus"
)
print(sum(assay(tse, "relabundance")[, "CL3"])) # [1] 0.3980999

Expected behavior
It's supposed to return 1 after the agglomeration.

Additional context
Independently verified with another dataset: https://ebi-metagenomics.github.io/HoloFoodR/articles/case_study.html#analyzing-microbiota-time-and-treatment-effects

[vsuksi]

Problem seems to be that there are lower- level taxa (e.g. strains) that don't have higher-level taxonomy. Use empty.rm = TRUE to keep them in agglomeration.

# Relabundance of non-empty fields in rowData(tse)$Genus
data(GlobalPatterns)
tse <- GlobalPatterns
tse <- transformAssay(tse, method = "relabundance") 
empty_genera <- is.na(rowData(tse)$Genus)
sum(assay(tse, "relabundance")[!empty_genera, "CL3"]) [1] 0.3980999

# Working example
tse <- agglomerateByPrevalence(
  tse,
  assay.type = "relabundance",
  rank = "Genus",
  empty.rm = FALSE
)
sum(assay(tse, "relabundance")[, "CL3"]) # [1] 1

[antagomir]

One question is what should be the default.

I agree that this can be confusing, especially because also the rows with missing genus name are real hits, real data. Filtering them out by default in the sum calculation can be misleading.

Can we also show both in the examples?

[artur-sannikov]

@vsuksi, thanks, it does solve the issue. However, it breaks the plotAbundance functionality now.

library(mia)
data("GlobalPatterns", package = "mia")
tse <- GlobalPatterns
tse <- transformAssay(tse, method = "relabundance")
tse <- agglomerateByPrevalence(
  tse,
  assay.type = "relabundance",
  rank = "Genus",
  empty.rm = FALSE,
  # Prevalence and detection are  to reduce the number of taxa
  prevalence = 20 / 100,
  detection = 0.2 / 100
)
sum(assay(tse, "relabundance")[, "CL3"]) # [1] 1

plotAbundance(tse, assay.type = "relabundance", order.row.by = "abund")

The code above just outputs a bunch of taxa squares with different colors, while the example of HoloFoodR merges all taxa in one line, so they look like "NA_Fusobacteria_Fusobacteriia_Fusobacteriales_Fusobacteriaceae_Cetobacterium".

[TuomasBorman]

One question is what should be the default.

I agree that this can be confusing, especially because also the rows with missing genus name are real hits, real data. Filtering them out by default in the sum calculation can be misleading.

Can we also show both in the examples?

empty.rank (onRankOnly) was discussed here, for instance: #419 (comment)

I think that it is more intuitive to have empty.rm = TRUE by default; user gets only taxa from requested rank.

@vsuksi, thanks, it does solve the issue. However, it breaks the plotAbundance functionality now....

The reason is that there are taxa which do not have information on the Kingdom level even though they have lower level information. This kind of taxonomy table is not common.

Some of the datasets from MGnify includes "Kingdom"-level information on "Superkingdom" column.

[TuomasBorman]

The most straight-forward solution is to do first the agglomeration and then data transformations. I think that order is most justified and avoids problems in agglomeration of transformed data.

[artur-sannikov]

Some of the datasets from MGnify includes "Kingdom"-level information on "Superkingdom" column.

So should we add the Kingdom level manually to avoid this issue?

[TuomasBorman]

Some of the datasets from MGnify includes "Kingdom"-level information on "Superkingdom" column.

So should we add the Kingdom level manually to avoid this issue?

The problem in your case was that the transformation was not applied after the agglomeration to prev_genus data:

Modify the according lines (add altexp = altExpNames(mae[[2]]))

#| label: transformation_metagenomics

# Transform microbiome with centered log-ratio method
mae[[2]] <- transformAssay(
    mae[[2]],
    assay.type = "counts",
    method = "relabundance",
    MARGIN = "cols",
    altexp = altExpNames(mae[[2]])
)
mae[[2]] <- transformAssay(
    mae[[2]],
    assay.type = "counts",
    method = "clr",
    pseudocount = TRUE,
    MARGIN = "cols",
    altexp = altExpNames(mae[[2]])
)

That should align with the goal of the analysis