This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Making more complex design in limma removes significant genes

I have a proteomics data set from 35 patients divided into three treatment groups (placebo, drug 1, drug 2) and recorded at 5 visits (timepoints) with the first visit making a baseline. Here is an illustration of the experimental design:

> library(tidyverse)
> metadata <- expand_grid(treatment = c("placebo", "drug1", "drug2"), visit = as.character(1:5)) |>
  mutate(subject = list(sprintf("%02d", 1:12))) |>
  unnest(subject) |>
  unite(participant, c(treatment, subject), remove = FALSE) |>
  select(-subject) |>
  arrange(participant)

> metadata
# A tibble: 180 × 3
   participant treatment visit
   <chr>       <chr>     <chr>
 1 drug1_01    drug1     1    
 2 drug1_01    drug1     2    
 3 drug1_01    drug1     3    
 4 drug1_01    drug1     4    
 5 drug1_01    drug1     5    
 6 drug1_02    drug1     1    
 7 drug1_02    drug1     2    
 8 drug1_02    drug1     3    
 9 drug1_02    drug1     4    
10 drug1_02    drug1     5    
# 170 more rows

First, I tried a simple model ~ treatment + visit, which returned about 400 proteins significantly changing at FDR < 0.05 level, mostly from one of the drugs. My approach was

> design_mat <- model.matrix(~ treatment + visit, data = metadata)
> tab |>
  limma::lmFit(design_mat) |>
  limma::eBayes()

where tab is the matrix with median-normalised log-intensities. Then, I looked a model with interactions: ~ treatment * visit. And now all the significant proteins, but one, disappeared.

As these approaches ignore longitudinal design of the experiment with the same patient giving 5 samples in 5 visits, I tried a block model:

> block <- metadata$participant
> corfit <- limma::duplicateCorrelation(tab, design_mat, block = block)
> fit <- tab |>
  limma::lmFit(design_mat, block = block, correlation = corfit$consensus.correlation) |>
  limma::eBayes()

This approach resulted in no statistically significant proteins.

My question is as follows. Are the results from the first, simplest model, suspicious? There is no interaction or patient blocking in this design. In the two more complicated model, do I lose statistical power by including more model components?

Here is one of the proteins picked up as significant by the simple model (Benjamini-Hochberg FDR = 6e-13), but not by the more complicated models (FDR = 0.2 and 0.96, for interactions and block design, respectively). To my eye, the change between placebo and drug1 looks real.

enter image description here

limma
• 2,565 views • link
updated by Gordon Smyth Gordon Smyth • written by marek.gierlinski • marek.gierlinski

When you refer to "significant proteins", significant for what comparison? You don't show any code that tests any hypotheses.

The interaction model is a totally different model from the additive model, and the coefficients have different interpretations. There are no coefficients in the interaction model that correspond to the treatment term in the additive model, so I suspect that you have not just fitted a more complex model, but tested a different hypothesis as well. As you may already know, the limma authors strongly recommend that you fit interaction models using a oneway layout in order to simplify and clarify the process of forming constrasts. Most analysts are not skilled at forming meaningful contrasts from factorial models. In our opinion, the interaction model concept of a "main effect" has little scientific relevance to omics analyses.

If you're using limma for mass spectrometry data, it would make sense to use the limpa package (https://doi.org/doi:10.18129/B9.bioc.limpa), which gives an enhanced limma pipeline specifically for this type of data.

• 0 views • link
written by Gordon Smyth Gordon Smyth

Thank you for your answer. To find "significantly changing" proteins I tabulate the p-values and fold changes from eBayes() for all coefficients except the intercept. See the code snippet below.

Would you recommend to use specific contrasts driven by the biological questions instead? I tried this, with no significant results. For example one of the questions is to compare proteins between drug and placebo at each visit. I used the following approach. This is an example for just one contrast, but I tested multiple contrasts - at each visit and comparing drug1 and drug2 vs placebo:

metadata <- metadata |>
  unite(group, c(treatment, visit), remove = FALSE)
groups <- unique(metadata$group)
design_mat <- model.matrix(~ 0 + group, data = meta)
colnames(design_mat) <- groups

ctr <- "drug1_5 - placebo_5" # other visits tested as well
contrast_mat <- limma::makeContrasts(contrasts = ctr, levels = design_mat)
fit <- tab |>
  limma::lmFit(design_mat) |>
  limma::contrasts.fit(contrasts = contrast_mat) |> 
  limma::eBayes()

In any approach I would proceed with the following to tabulate the results for all coefficients, except for the intercept:

coefs <- colnames(fit$coefficients) |> str_subset("Intercept", negate = TRUE)
map(coefs, function(ctr) {
  limma::topTable(fit, coef = ctr, number = 1e6, sort.by = "none") |>
    as_tibble(rownames = "id") |>
    add_column(contrast = ctr)
}) |>
  list_rbind()
• 0 views • link
written by marek.gierlinski • marek.gierlinski

Your code doesn't correspond to either of the design matrices defined in your question. I will assume that the interaction formula shown in your question was just a red herring, and you actually formed visit-specific contrasts as shown above.

• 0 views • link
written by Gordon Smyth Gordon Smyth

Apologies for not being clear. What I wanted to say was that I tried different approaches: first with additive, then with interaction formula (the original question) and then I also tried a fews models with selected contrasts (the followup).

The second part of the followup post (the one with topTable()) is my attempt to answer how I found "significant" proteins. You said "You don't show any code that tests any hypotheses". Correct me if I'm wrong: eBayes() performs statistical testing and returns t-statistics and the corresponding p-values for all coefficients under the null hypothesis that the mean(coefficient level) - mean(coefficient baseline) = 0. This is how I understand it. I simply extract the p-values using topTable(). I don't perform any additional tests.

• 0 views • link
written by marek.gierlinski • marek.gierlinski

2 answers

Yes, the results from the simplest model are suspect because that model ignores the repeated measures on the patients, so it is over-stating (by a factor of 5) the number of independent observations in your dataset.

On the other hand, I am not that convinced about the need to adjust for visit, because the visit times are likely to be participant-specific and because the plot shown in your question does not indicate any discernible visit effect. I also don't see much merit in comparing drugs by visit time, because that approach throws away lots of statistical power and (I suspect) is not of primary biological interest or relevance. No one is going to be interested in a drug that changes or reverses its effect from one visit to another. You could be interested in a drug effect that develops over time, but that can be more efficiently tested by a different model -- see below.

Meanwhile, there are also lots of ways to increase power in limma, such as sample weights and variance trends. There is also the limpa package, that I mentioned above, which is tuned to mass spec data and which is more powerful than limma alone.

In limpa, I would be trying an additive model with blocking:

treatment <- factor(treatment, levels=c("placebo","drug1","drug2"))
visit <- factor(visit)
y <- dpcQuant(peptidedata)
design <- model.matrix(~treatment+visit)
fit <- dpcDE(y, design, block=participant, sample.weights=TRUE)
fit <- eBayes(fit, robust=TRUE)
topTable(fit, coef=2) # drug1
topTable(fit, coef=3) # drug2
topTable(fit, coef=4:8) # visit effect

You have said that an additive model is difficult to interpret, but I don't see why. It simply tests for a consistent drug effect across the visits.

If the baseline visit effect gives no significant results, then I would consider removing it from the model and just using

design <- model.matrix(~treatment)

If you want to specifically test for a drug effect that increases over repeated visits, I would do it like this:

visit <- numeric(visit)
design <- model.matrix(~treatment+visit+treatment:visit)

In this model, coefficients 2 and 3 correspond to the baseline drug effects at the first visit, while coefficients 5 and 6 correspond to drug vs placebo trends over visits.

• 0 views • link
written by Gordon Smyth Gordon Smyth

I am aware of limpa, but at the moment I have an entire pipeline developed for proteomics based on limma, so it was easier for me to do it this way. But I will definitely give limpa a go later.

I have tried all the approaches you outline in this post. The one with participant as a block variable (see my original post) removes a lot of significance (in the original post I stated mistakenly that there were no significant proteins, which was due to a bug in my code). Here is what I get. Below I quote the numbers of significant (FDR < 0.05) proteins.

  • For the model ~ treatment + visit: 257 (drug1), 247 (drug2) and 8 across all visits.
  • For the model ~ treatment + visit with participant as a block variable: 3 (drug1), 0 (drug2) and 33 across all visits.

From duplicateCorrelation() I get the consensus correlation of 0.38.

Essentially, including participant as a block variable destroys nearly all drug effect. Does it mean the results from the non-block model are mostly false positives? I think I understand now what you mean by saying "[non-block] model ignores the repeated measures on the patients, so it is over-stating (by a factor of 5) the number of independent observations in your dataset."

Does it have to do anything with the fact that participant and treatment are confounded (each participant receives only one treatment)?

• 0 views • link
written by marek.gierlinski • marek.gierlinski

I have tried all the approaches you outline in this post.

You don't mention having tried sample weights or robust empirical Bayes.

Does it mean the results from the non-block model are mostly false positives?

The non-block model is not a correct representation of the true nature of the data, and it over-states significance, so the true FDR may be much higher than the nominal 5% rate, perhaps much higher. The adjusted p-values from the blocked model for the same proteins will give a better idea of the likely true FDR rate.

Does it have to do anything with the fact that participant and treatment are confounded (each participant receives only one treatment)?

Yes. The comparisons of interest (drug vs placebo) are between participants, so only the between-subject variability is relevant, not within-subject variability. Yet the non-block statistical model is estimating variability primarily from the within-subject variability. The majority of the residual degrees of freedom arise from within-subject (between visit) comparisons. The blocking model is statistically somewhat subtle, but it builds in the fact that within-subject comparisons are likely to be more consistent than between subjects, hence putting the between and within subject comparisons on the same footing and reducing the bias in the variance estimation.

I will definitely give limpa a go later.

limpa will give most benefit when run on peptide precursor level data, but it can also work ok at the protein level. If you have a limma EList of protein LFQ values, including missing values, you can convert it to a limpa object by:

yimp <- dpcImpute(y)

Then everything will work exactly as for limma except that you use dpcDE() instead of lmFit(), and you don't need to run duplicateCorrelation() or arrayWeights().

• 0 views • link
written by Gordon Smyth Gordon Smyth

Thank you for this clarification. Your clear explanations and generous help are greatly appreciated by everyone here.

• 0 views • link
written by marek.gierlinski • marek.gierlinski

Adding more terms to your model can decrease statistical power, because the same amount of data is being used to estimate more parameters. In general, you will get more statistical power if you include terms that are explainitory, and lose statsitical power if you add terms that are not explaintory.

Its not clear to me from your description what contrasts it is you are interested in, and what your expected behavior for the data is.

I assume when you test the ~treatment + visit model you performing an annova-like test on the "treament" coefficint, with the reduced model being "~visit".

When you test the ~treatment * visit model are you testing the "treatment" terms, or the "treatment:visit" terms? In general, I'd not add an interaction term to a model unless you were going to test it.

What is the intresting biological question here? Are you interested in identifying proteins that change over the course of several visits, proteins that are different between the drugs, or proteins where the difference between drugs changes over a series of visits?

• 0 views • link
written by i.sudbery i.sudbery

To answer your question about the biological questions: yes! The biologists I work with are interested in the comparison between drugs and placebo at each visit, in how the proteome changes over the course of visits for each branch (placebo, drug1, drug2) and how the change in each branch compare between branches. To answer these specific questions I tried various contrasts (see my reply to Gordon for an example). For the last biological question I tried something like (drug1_5 - drug1_1) - (placebo_5 - placebo_1).

Alas, none of the contrast approaches returned any significantly changing proteins. I suppose this due to the loss of statistical power. For example, when testing contrast drug1_5 - placebo_5 I have 24 samples instead of 180, as in the full additive model. The only approach returning any significantly changing proteins is the full additive model ~ treatment + visit. But this one is difficult to interpret.

• 0 views • link
written by marek.gierlinski • marek.gierlinski

Log in to answer this question.

More posts like this