Questions about choosing analysis method Masslin3

Downstream analysis and statistics MaAsLin
1

Hi,

There are many choices for analysis methods, normalizations, and transformations in Masslin3. I don’t know which will be the optimal methods for my data analysis. My taxa table includes relative abundance values in percentage (ranged 0-100 not 0-1) at species level. Values of each taxa is summed up to 100. By checking the histogram of some species in my data, over 70% of the relative abundance values of a certain species are 0 (highly skewed, and zero inflation). Several bigger values above 5% in some species (outliers?).

Based on the tutorial, I think the optimal analysis method should be CPLM, while I saw the below in the tutorial of Masslin2: Regarding whether to use LM , CPLM , or other models, intuitively, CPLM or a zero-inflated alternative should perform better in the presence of zeroes, but based on our benchmarking, we do not have evidence that CPLM is significantly better than LM in practice. Which method should I try?

Following the analysis method selected from the above, which normalization is suitable? Is the default method TSS good enough?

For the transformation, is AST the best choice? Do I need to scale my original data from 0-100 to 0-1 before conducting AST in Masslin3()?

I would appreciate it if someone could provide any suggestions. Thank you!

2

Hi,

The defaults in MaAsLin 3 (TSS normalization, no filtering, log transforming, split prevalence and abundance analysis) will probably be best for your data as they are for most scenarios. TSS normalization (relative abundances), log (base 2) transformations, and splitting prevalence (presence/absence) from abundance (how much is there, if it’s there) make for very clean interpretations, and we’ve extensively benchmarked these and found them to perform well. To answer your specific question, the TSS normalization will take care of converting your 0-100s to 0-1s.

You might decide to turn on or off the median_comparison_abundance parameter based on reading this section, but I’d recommend leaving median_comparison_prevalence off.

Let me know if you have additional questions.

Will

3

Hi Will,

Thank you for the reply! I tried the models using TSS and LOG in MaAslin3. The features I want to test are time, treatment and their interactions. I used two different formulas. One is ~ Time + Treatment + Time:Treatment + (1|SubjectID). I have six samples from the same subject. Then I created a new interaction column in my metadata named as Trt_time to manually combine the time and treatment. Then run another model with ~ Trt_time + (1|SubjectID). While the two models have different outputs of species (some overlapped), p and q values in significant results for the interaction. Why? I set the reference level for the model 1 as reference = c(‘Treatment,TrtA’ , ‘Time,baseline’) and model 2 as reference = c(‘Trt_time,TrtA_baseline’).

Another question for the model 2: If I set the reference level: reference = c(‘Trt_time,TrtA_baseline’). Any significant associations in the final output were generated by only comparing others with the reference level, right? Then I created a contrast matrix to compare different treatments at the same time points (e.g., trtA_baseline vs trtB_baseline). In my significant results for the model ~ Trt_time + (1|SubjectID), I has one row/species having qval_joint<0.01 for trtB_baseline, I believe this siginificant result is against my reference ‘trtA_baseline’. While in my contrast test output, all the p and q values in the rows of abundance and prevalence results are NA. Is there any possible error? I created the contrast matrix using ‘name’ column in the model result output.

4

Regarding your two models, assuming both time and treatment are categorical variables, both of those models will fit a coefficient to each category. However, what that coefficient means and how significance is tested is very different.

In the first model, [the time groups are tested against the baseline time] for the baseline treatment group, [the treatment groups are tested against the baseline treatment] for the baseline time group, and the interactions then test whether the treatment effect differs from its effect at the reference time. It’s quite complicated and requires being very careful about your interpretations.

In the second model, all groups are simply tested against the single baseline time-treatment group (which may or may not be what you care about).

Regarding your second question for model 2, that is correct: everything is compared against the overall reference. Regarding the NA contrast test, can you provide a reproducible example that I can test on my end? It’s hard to say whether the issue is in defining the contrast matrix, in trying to get coefficients for treatment-time groups that are never present, or in the MaAsLin 3 code itself.

Will

5

For the answer 1, I would like to double check if I understand correctly. I have two levels in treatment : TrtA (reference) and TrtB. I have three levels in time: day0 (reference), day10, day22. If I build a model as below:

DAmodel1= maaslin3(input_data = otu_df,
input_metadata = metadata,
output = “output_damodel1”,
formula = ‘~ Treatment + Time + Treatment:Time + (1|SubjectID)’,
normalization = “TSS”,
transform = “LOG”,
reference = c(“Treatment,TrtA”, ‘Time,day0’),
standardize = FALSE,
augment = TRUE,
min_prevalence = 0.1,
max_significance = 0.1,
median_comparison_abundance = TRUE,
median_comparison_prevalence = FALSE,
max_pngs = 250)

For the main treatment effects: It compared TrtB with the reference TrtA across all the times? You mentioned [the treatment groups are tested against the baseline treatment] for the baseline time group. It makes me a little confused, the baseline time group here is only day0 group or regardless the different time points?
For the main time effects: My understanding is to compare day10 with day0 and compare day22 with day0 regardless different treatments here.
For the interaction effects: Compare TrtA_day10 vs TrtA_day0, TrtA_day22 vs TrtA_day0, TrtB_day10 vs TrtB_day0, and TrtB_day22 vs TrtB_day0, right?

For the model2, my R codes are below: Based on the levels of treatments and times, I have six different levels in Trt_time.
DAmodel2= maaslin3(input_data = otu_df,
input_metadata = metadata,
output = “output_damodel2”,
formula = ‘~ Trt_time + (1|SubjectID)’,
reference = c(‘Trt_time,TrtA_day0’),
normalization = “TSS”,
transform = “LOG”,
standardize = FALSE,
augment = TRUE,
min_prevalence = 0.1,
max_significance = 0.1,
median_comparison_abundance = TRUE,
median_comparison_prevalence = FALSE,
max_pngs = 250)

contrast_mat ← matrix(c(
1, -1, 0, 0, 0, 0, # TrtA_day0 vs. TrtB_day0
0, 0, 1, -1, 0, 0, # TrtA_day10 vs. TrtB_day10
0, 0, 0, 0, 1, -1, # TrtA_day22 vs. TrtB_day22
0, 1, 0, -1, 0, 0, # TrtB_day0 vs. TrtB_day10
0, 1, 0, 0, 0, -1 # TrtB_day0 vs. TrtB_day22
), ncol = 6, byrow = TRUE)

colnames(contrast_mat) ← c(“Trt_timeTrtA_day0”, “Trt_timeTrtB_day0”,
“Trt_timeTrtA_day10”, “Trt_timeTrtB_day10”,
“Trt_timeTrtA_day22”, “Trt_timeTrtB_day22”)

rownames(contrast_mat) ← c(“Baseline_Comparison”,
“10day_Comparison”,
“22day_Comparison”,
“TrtB_day0_vs_4day”,
“TrtB_day0_vs_22day”)

contrast_out ← maaslin_contrast_test(maaslin3_fit = DAmodel2,
contrast_mat = contrast_mat)

significant_abundance_results ← contrast_out$fit_data_abundance$results %>%
filter(qval_individual <= 0.05 | qval_joint <= 0.05) #####???Got 0 rows. All the q and p values are NAs, same for the prevalence results.

6

For model 1, the coefficient that’s reported as TrtB will be the difference between TrtA and TrtB only in the baseline time group (day0). It’s not comparing them across all time points because of the interaction effect you have. If you want to compare across all time points, you should use ~ Treatment + Time + (1|SubjectID), which assumes that the treatment effect is the same at each time but that each time has a difference from day0 that applies equally to both treatment groups.

The time effect as you have it here works analogously, so if you want to compare day10 with day0 and day22 with day0 regardless of treatment, you should not use the interaction term.

The interaction terms you have right now will be:
TrtB:day10, TrtB:day22 which gives the difference in the day10 coefficient for the TrtB group relative to the TrtA group and likewise for day22. The other terms will appear as main effects in your output because they involve either the reference time or the reference group. You can check this in the name column of your results.

Are the matrix names you’re using for the contrast matrix actually all in the names column of your maaslin3 output table? Trt_timeTrtA_day0 for instance should be the baseline (intercept) and therefore not show up in the output table. Does the error column in the contrast fits say anything like Predictors not in the model had non-zero contrast values? Some error should be reported.

If it’s still unclear why the contrast tests are not working, can you provide a small example of the data that I can run so I can evaluate what the issue is?

Will

7

Thank you so much for the detailed explanation, then now I understand why I got different outputs than what I expected. And you are right, I just noticed that I had ‘Predictors not in the model had non-zero contrast values’ in the error column because I set wrong names in the contrast matrix. I also asked a statistician in my department. She suggested not only using an interaction term or a variable representing interactions (like Trt_time in my model2) to build a model. So I think I will stick with the formular like ~ Treatment + Time + Treatment:Time. While based on your explanation, I don’t think how I set up the model and reference now can compare the pairs I want.

To reclarify my study aims: Two treatment including TrtA(Control) and TrtB(Interested treatment). I used Control and Interest instead of TrtA and B for better description in the below. Three time points: day0(Baseline), day10 and day22. The comparisons I want to get in the output are:
Control vs Interest at day0,
Control vs Interest at day10,
Control vs Interest at day22,
Interest_day10 vs day0,
Interest_day22 vs day0,
Interest_day22 vs day10. I want to build a correct model and may need do level wise comparisons by a contrast test based on the correct model. How do I set up the formula and reference in the model? Thank you!

8

Glad you got that figured out with the contrast matrix!

To make those comparisons, you’ll want to fit ~ Treatment + Time + Treatment:Time with Control and day0 as the references and then analyze as follows (the coefficients here won’t literally be what’s in the name column because I’m not sure what your levels are called, but it should be close):

  1. Control vs. interest at day0: TreatmentInterest
  2. Control vs. Interest at day10: TreatmentInterest + TreatmentInterest:Timeday10
  3. Control vs. Interest at day22: TreatmentInterest + TreatmentInterest:Timeday22
  4. Interest_day10 vs day0: Timeday10 + TreatmentInterest:Timeday10
  5. Interest_day22 vs day0: Timeday22 + TreatmentInterest:Timeday22
  6. Interest_day22 vs day10: Timeday22 + TreatmentInterest:Timeday22 - TreatmentInterest:Timeday10 - Timeday10

Put the things that involve sums or differences in the contrast matrix to get proper standard errors and p-values. Linear models with interactions can get thorny :slight_smile:

You can check that these are as expected by doing the addition/subtraction on the coefficients in this toy example:

df <- data.frame('Time' = factor(c('day0', 'day10', 'day22', 'day0', 'day10', 'day22'), 
                                   levels = c('day0', 'day10', 'day22')), 
                 'Treatment' = factor(c('Control', 'Interest', 'Control', 'Interest', 'Control', 'Interest'), 
                                   levels = c('Control', 'Interest')), 
                 'outcome' = c(1, 3, 7, 11, 16, 15))
lm(outcome ~ Time + Treatment + Time:Treatment, df)

Also, it might be easier to do all of this by repeatedly fitting the MaAsLin 3 model with different references and just taking the comparisons you care about. It’ll take longer to run, but the risk of messing something up is lower, and you’ll be more confident in what you did.

9

Thank you for the info, it’s really helpful! Now I got the correct models and outputs for the contrast tests. I would like to make sure if I interpret the significant results correctly. To check the significant differences associated with prevalence or abundance, I look at the coefficients and q values. For example to compare Control vs. interest at day0: TreatmentInterest, if the model is ‘prevalence’, the coefficient is positive and q_individual is very small, it means the prevalence of the correspondering species in Interest group is significantly more than that in Control, right? Then if the model is ‘abundance’, it means the relative abundance of the species in Interest group is significantly higher than that in Control.

Questions 2 about the q_joint, what I understand is it indicates the significant association with either prevalence or abundance. So for the same species with the same value and ‘name’, the q_joint values are the same no matter it’s prevalence or abundance model? Then e.g., if I get a q_individual>0.05 but q_joint<0.05 in a prevalence model, it suggests that the abundance model for the same species and same ‘name’ has a small q individual value <0.05? I am trying to understand the relationship and difference between individual and joint q values.

Question 3 is I got the some big coefficients above 10 (2 absolute values even more than 15, see as below), which are suspicious based on the tutorial. But I didn’t receive any errors, the metadata are not collinear, and 15 samples in each group. I got q values = 0 (both individual and joint) with the big coefficients. How to explain the 0 value here, very different so the q value is extremely small? Or is it still reliable to trust these results?

image
image995×241 6.28 KB

Appreciate for your kind help!