Quick chat: Ask short questions or share comments

[TuomasBorman]

This thread is for quick questions and comments! If you have a short query or just want to share a quick thought without starting a new thread, feel free to post it here.

[mniku]

Stupid question: how do I calculate fold changes for clr transformed values (e.g. counts of microbiota abundances)? Relates to my recent discussion with Leo... I'm trying to compare taxon-wise FCs in our longitudinal data, and Leo recommended using clr transformed data. But I have difficulty understanding the transformed values and how to properly calculate the FCs from them --- doesn't seem intuitive when I compare to simple relative abundace based FCs.

[mniku]

Ah, OK, naturally.. thanks for the super quick reply!

[mniku]

Another stupid question about CLR transformed microbiota count tables. Is it okay to use them for calculating Spearman or Pearson correlations of abundances of a taxon between samples, or is it better to use another transformation or just the relative abundances?

[antagomir]

Do you mean correlations between features (taxa), or correlations between samples?

[mniku]

Sorry, tried to formulate the question clearly but failed :D Between, for example, abundance of Bifidobacterium in sows, and abundance of Bifidobacterium in their piglets. To see if its abundance in sows is associated with its abundance in their own offspring.

[antagomir]

Ok so it seems like correlation between paired samples.

Using correlation between CLR-transformed data should be an acceptable approach for this imo and better than with plain relative abundances.

[mniku]

I'm reading the NetCoMi section of the OMA handbook and I'm a bit confused... is the current example output from a different analysis vs. the accompanying text? The text says almost everything in the output is significant, but the p values don't seem like that...

[TuomasBorman]

Hello!

Sorry for caused confusion. You are probably referring this: https://microbiome.github.io/OMA/docs/devel/pages/network_comparison.html#quantitative-comparison

I tag @stefpeschel as she knows the NetCoMi the best

[mniku]

Yes exactly that part!

[mniku]

@stefpeschel or did I just get it the wrong way... is it actually so that a small p value in the network comparison results mean that the two networks are similar (at least for some of the outputs)? Could you kindly clarify this a bit?

[stefpeschel]

Hi!
Please excuse the confusion. I just realize that there must be an error in the output because all network measures are doubled in the summary, which should not be the case. I will revise the section on network comparison right away.

You are right that small p-values should represent significant differences (the null hypothesis is that the values are equal in both groups). Generally we won't expect significant results for the centrality measures and the associations itself after multiple testing correction due to the low number of permutations. For the global properties, however, we usually see significant differences because there is not correction needed.

[mniku]

Thanks! I’m still struggling a bit in interpreting the output. May I ask if you know any studies which used NetCoMi (or other network analysis tools) to analyze effects of experimental antibiotics on gut microbiota? Or some other experimental, controlled manipulation (comparing inherently similar groups after perturbation)?

[stefpeschel]

You can find explanations on how to interpret the results in the NetCoMi paper itself:
https://academic.oup.com/bib/article/22/4/bbaa290/6017455

A publication I can recommend is this one by Alice Sommer:
https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1010044

Here's a list with papers that cited NetCoMi:
https://scholar.google.de/scholar?oi=bibs&hl=en&cites=16729940613280981842&as_sdt=5
Maybe you can find another publication that is close to your research topic. Unfortunately, I have not looked at most of the papers in detail to be able to make a recommendation.

[mniku]

Thanks Stefanie! After experimenting a bit, I'm slightly worried about the robustness of the results with my data - I got some very interesting observations, but it seems the results vary quite a lot depending on parameters (such as taxon filtration, number of samples included, etc). I have so far been mostly using the SpiecEasi method, now testing with SPRING. Do you think my n is just too small (I have 3 experimental groups which are inherently fairly homogenous, with 22-30 animals per group)? Is there something I can do with NetCoMi parameters to maximize robustness and reproducibility?

[TuomasBorman]

Was this solved? @mniku @stefpeschel

[stefpeschel]

Hi! Please excuse the late reply. I was on vacation for two weeks and didn't find the time to answer your question.

You're absolutely right: microbiome network results can be sensitive to preprocessing and inference parameters – especially with small sample sizes.

In our paper Ullmann & Peschel et al., 2023, we showed that different pipelines can produce quite variable networks from the same data. This is even more pronounced when working with limited sample sizes like yours, where only the strongest associations are reliably detected.

One solution would be correlation shrinkage as discussed in Badri et al. 2020, which is not implemented yet in NetCoMi.

A few suggestions for improving robustness in NetCoMi:

• Filter taxa more aggressively (e.g. keep only those present in ≥30–40% of samples)
• Use SPRING alongside SpiecEasi - it might perform better with small n
• Both, SpiecEasi and SPRING use the StARS stability-based model selection approach. There you can vary the threshold, nlambda, and/or rep_num parameters to get more robust results.
• Use different estimation methods and report a consensus network, where you only keep the edges occuring in all networks.

I hope these ideas help to get more robust results. In any case, your networks will remain sensitive to changes in the parameters for such low sample sizes.

Best,
Stefanie

[mniku]

I'm puzzled about using fixed vs random effects in the differential abundance tools (Maaslin3, ANCOMBC2).

1. It seems the OMA handbook only uses fixed effects although some of the confounders seem more like random effects to me (disclaimer: I'm really ignorant about multivariate statistics!).
2. When I try to use random effects in my Maaslin3 or ANCOMBC2 models (yes same thing in both), I lose most of the significant differences (while if I use the same confounders as fixed effects, the differences get MORE significant).

In this case, I'm analyzing piglet microbiotas from a study with 3 experimental groups, each with ~25 piglets from ~6 dams. I’m trying to use dam parity (6 levels) and room (2 levels) as confounders (as they probably affect the microbiotas, although we aren’t specifically interested in these effects). Could the problem be that I have too small N for too many levels? Is it especially problematic for random effects?

Our statisticians say this kind of confounders MUST be used as random effects. I also asked Maaslin developers and they said it's up to me and either works. So I'm pretty clueless now :D

This is my formula using fixed effects: ~exp_group+parity+room+reads
And random effects: ~exp_group + reads + (1|parity)+(1|room)

[antagomir]

The mixed effect examples in OMA could probably be improved. I would do parity as random effect but if room only has 2 levels then fixed effects might be a better choice (potentially less convergence issues, better interpretability, also the 2 levels should be truly randomly sampled in principle, which I don't think is the case here - usually 2-level cases would go for fixed effects but this can vary).

[mniku]

Thanks. Two more questions about parity:

1. Why should it be used as random rather than fixed? I think it's not strictly a clustering/grouping thing, as we know it affects the sow microbiotas, which then probably affect the piglet microbiotas. Tempted to ask because for some reason both Maaslin and ANCOMBC show much more significant experimental group differences if it's used as fixed (yes I know this is not they way to select models, just trying to understand how it works).

2. Would you use it as continuous or categorical? Alpha diversity seems to decrease pretty linearly by age, but I guess some taxa could very well have nonlinear associations. On the other hand, there are a bit too many categories as such (parity 1-8) and no obvious way to pool categories.

[antagomir]

1. It is common to assume that the effect of subject/parity/such is randomly (normally) distributed. If you can make that assumption then random effects model is justified. Fixed effect is also OK but in general you would need more data to power the analysis because the number of variables is higher (in this case, 1 vs 6). Fixed effects might be more easy to interpret.

2. "Would you use it as continuous or categorical?" -> I am not sure what "it" refers to in here.

In general, this OMA forum is mainly for technical advise on using Bioconductor methods, not for general statistical consultation because that is somewhat out of scope for the technical methods ecosystem and might require more in-depth expert examination of each given application.

[jessieellis]

Hi, I have a quick question on Chapter 10.5: Filtering out zero-variance features.

In the example to remove features with zero variance, features with an sd > 1 are selected to subset, with the explanation that 1 is chosen because log(1)=0. However, I am not quite following this for two reasons:

1. it appears that only the plot is log-transformed, not the calculated standard deviations. So it seems that if the intention is to select features with variances greater than zero, shouldn't sd >0 (not 1) be specified?
2. If I am wrong and the calculated sd variable (and not just the plot) is already log transformed, it still seems wrong to take sd > 1? Since values between 0 and 1 become negative once log-transformed, selecting a log-transformed value of greater than 1 would filter out more than just zero variance features - it would also filter out features that have an sd between 0 and 1.

Please please correct me if I'm wrong! I'm just trying to make sure I understand properly for when I go to conduct this step on my own data. Thanks!

[antagomir]

@raivo-otus could you have a look?

[raivo-otus]

Hey @jessieellis,

Thank you for the question, you are indeed correct. The section discusses removal of zero-variance features, however, the code example essentially filters out low-variance features.

In the chapters example code the sd values are log-transformed for the histogram plot. When selecting; the sd values are not log transformed. The filter of sd > 1 will also filter out features with low, but non-zero variance. If the aim is to remove true zero-variance features the filter should be sd > 0.

Additionally using sd values from raw counts is somewhat suspicious as read depth is not normalized across samples. This should not be an issue when filtering out zero-variance features, but does influence what is considered low-variance.

The motivation to filter out zero-variance features remains sound. I will raise an issue on this section and we will work to correct the section and improve clarity. Thank you for pointing this out!

[jessieellis]

Thanks so much for your quick review and reply!

[jessieellis]

Hi there; dumb question: If I apply a transformation to alternative experiments (such as in Chapter 12, Exercise 8 after agglomerating by ranks), where are the transformed data stored? I can't find them!

[TuomasBorman]

Hello @jessieellis, great question!

When you apply the following command

tse <- agglomerateByRanks(tse)

you will see that altExp slot is now populated with agglomerated data. You can retrieve the agglomerated experiment with altExp(tse, "name_of_the_experiment"). It is TreeSE, but now agglomerated. So it is TreeSE inside TreeSE (this might be confusing).

When you apply

tse <- transformAssay(tse, altexp = "name_of_the_experiment", method = "relabundance")

the transformation is applied to all experiments, including agglomerated data in altExp.

Answer to your question:

The transformed data is here:

altExp(tse, "name_of_the_experiment") |> assay("name_of_your_transformation")

-Tuomas

[jessieellis]

Thank you for the explanation!

[antagomir]

Yes, or just assay(altExp(tse, "name_of_the_experiment"), "name_of_your_transformation").

So it is nested indeed, alternative tse inside the main tse.