Transcripts with beta close to 0 are considered differentially expressed in a likelihood-ratio test

[apcamargo]

I'm comparing the results that I obtain when doing a DE analysis with the Wald test and the likelihood-ratio test. One the thing that I've noticed is that there are many genes with 'beta' close to zero that are considered differentially expressed between the conditions.

I know that the likelihood-ratio test does't use the beta values to calculate the p-values, but I find it strange that transcripts with similar expression values between the conditions are being considered differentially expressed.

Volcano plot (Wald test q-values):

Volcano plot (likelihood-ratio test q-values):

[roryk]

Does your model just have two factors that you are testing, or is it a more complicated model?

[warrenmcg]

Thanks for the feedback! Two things: 1) I would be interested to see an example of the "transcript view" plot from sleuth_live (where we see the point estimates and the bootstrap variability) for one of these close-to-zero-beta transcripts. 2) Also, what is the experimental context here: comparing wild-type to knockout, a response to a drug or stress, comparing diseased cells to normal cells?

[apcamargo]

The model has two factors: tissue and stage. I'm trying to discover transcripts that are differentially expressed between stages, correcting for tissue effects.

This is the code:

model <- sleuth_prep(data, ~ stage + tissue)
model <- sleuth_fit(model)
model <- sleuth_fit(model, ~ tissue, 'reduced')

# Check the models:

models(model)
    [  full  ]
    formula:  ~stage + tissue 
    coefficients:
        (Intercept)
        stagenb
        tissuevno
    [  reduced  ]
    formula:  ~tissue 
    coefficients:
        (Intercept)
        tissuevno

# LTR test (to obtain the q-values):

model <- sleuth_lrt(vno_om_model, 'reduced', 'full')
results_lrt <- model(model, 'reduced:full', test_type = 'lrt')
results_lrt <- results_lrt[order(results_lrt$qval),]

# WT test (to obtain the beta values):

model <- sleuth_wt(model, 'stagenb')
results_wt <- sleuth_results(model, 'stagenb', test_type = 'wt')
results_wt <- results_wt[order(results_wt$qval),]

# Obtain a single data frame

final_results <- merge(results_lrt, results_wt, by="target_id")
final_results <- final_results[,c("target_id", "qval.x", "qval.y", "b")]
colnames(final_results) <- c("target_id", "qval_lrt", "qval_wt", "b")
final_results <- final_results[order(final_results),]

head(final results, 3)

    |        |   target_id    |   qval_lrt    |    qval_wt     |     b      |
    |--------|----------------|---------------|----------------|------------|
    | 168181 | TCONS_00175843 | 5.468570e-14  | 3.299774e-129  |  -7.704499 |
    | 187768 | TCONS_00196546 | 8.096083e-14  | 5.635809e-231  |   5.686079 |
    |  87186 | TCONS_00091156 | 4.401443e-13  | 6.273986e-178  | -3.670974  |

The plots were done using this data frame.

[apcamargo]

Here's one example:

|   target_id    |   qval_lrt    |    qval_wt     |     b      |
|----------------|---------------|----------------|------------|
| TCONS_00028547 | 0.0006324667  | 0.0001756688   |  0.5043041 |

[apcamargo]

It seems to me that sleuth is computing the betas wrong. Maybe it's calculating beta as if the comparison is between tissues.

[apcamargo]

This is what I've got when plotling sleuth LRT's q-values and Log2FoldChange values computed in DESeq2. This makes me think that the problem is not with the betas.

[warrenmcg]

Hi @apcamargo,

The "beta" value is attempting to get exactly what you describe: given your model of ~ stage + tissue (which assumes no interaction between the two factors), what is the effect of stage?

The example you give us is interesting. Looking at stage independent of tissue, there is wide variability in both groups, but there still seems to be ~2-3 fold increase in nb (newborn?) vs adult. Looking at tissue independent of stage (vno vs om), there appears to be no overall difference in the mean, but a huge difference in variability.

However, what is most interesting is that, with this example, there is a clear interaction effect: looking at the tissue plot, and assuming that the first 6 samples are from adult and the last 3 are from nb (based on comparing the tissue plot to the stage plot), there appears to be moderate differential expression between stages in OM, but there is clear and large differential expression in VNO. The fact that stage has a different effect in different tissues indicates an interaction. In this case, it is somewhat misleading to think about an overall effect from tissue without considering an interaction, because the effect is different depending on the context.

It would be interesting to see how your results change overall if you instead use the following model: ~ stage + tissue + stage:tissue (this includes an interaction effect) and compare it to ~ tissue. In this context, you would use LRT to get a sense of which genes behave more closely to the model that stage affects expression. You could also do the comparison of ~stage + tissue to ~stage + tissue + stage:tissue to see which genes are most likely exhibiting an interaction effect. If you would want betas to get a sense of "fold-changes", you would then use stagenb with the full model as your beta if there is a demonstrated interaction effect, or using the model without an interaction term if there is not a clear demonstrated interaction effect. However, remember that these would be fold-changes with an awareness that the effect of stage may be different depending on the tissue (in the case of including the interaction term) or with an awareness that tissue may have an effect independent of stage (in the case of excluding the interaction term).

[warrenmcg]

Addendum: the likelihood ratio test simply compares how the transcripts behave relative to the full model versus relative to the reduced model. In general, if the full model better predicts how the data behaves, the LRT produces a small p-value; in other words, a significant hit means the transcript's expression appears to be associated with both tissue and stage.

Your example with a beta of ~0.5 is roughly equivalent to a 40% increase if it was a straightforward fold-change, which is a reasonably moderate change.

[apcamargo]

Thank you very much, @warrenmcg ! I'll try to use a interaction term right now and will post the results here as soon as I can.

However, I've also done comparisions within each tissue, using '~1' as the reduced model. Even in this case, there seems to be transcripts with beta close to zero and low q-values. Look at the plot in the center:

[warrenmcg]

The scale makes it difficult to think about what's happening close to 0. Recall that I said that your example of ~0.5 corresponds to about a 40% increase if it was a straightforward fold-change, which is respectable and definitely worthy of consideration as biologically meaningful. Your scale goes from +10 to -10, so figuring out what's happening on the scale of 0-2 is difficult to see.

Two things you can try to get a better sense of what's happening are the following: 1) store the result of the plot, and then modify the x-axis like this:

plot_obj <- sleuth::plot_volcano(sleuth_object, [[whatever other options you use]])
plot_obj <- plot_obj + xlim(c(-2, 2))
plot_obj

We can then see in finer detail what's happening close to zero.

2) get the minimum absolute beta from your significant results:

results <- sleuth::sleuth_results(sleuth_object, [[whatever other options you use]])
sig_results <- results[which(results$qval <= 0.05), ]
min(abs(sig_results$b))

What is that value? If it's below ~25% up- or down-regulation (absolute value of b < ~0.32), that's when we can discuss specific examples of where the overall methodology of linear regression is picking up really small changes.

[apcamargo]

There are 9154 betas <=0.3 that belong to significant genes

[warrenmcg]

That's a lot! How many total significant transcripts/genes are there in your experiment?

[apcamargo]

36032 in total. This is the new test, using the interaction in the full model.

[apcamargo]

WT test gives me more credible results:

table(results_wt$qval<=0.05)
     FALSE  TRUE 
     87347  7256 

[warrenmcg]

And what's the minimum absolute beta for the significant wald test results?

sig_results <- results_wt[which(results_wt$qval<=0.05), ]
min(abs(sig_results$b))

[apcamargo]

0.34591856289752304

[warrenmcg]

Ah, that's expected and no issue then. It's a small beta, but above the 25% threshold I mentioned before. The Wald test specifically tests whether a particular beta is non-zero, so I would be surprised if gave many significant hits with betas less than 20-25%.

The betas you're seeing for the 9154 from before are simply from transcripts where the LRT is indicating that the data is better predicted by the full model than the reduced model. I'm not sure which comparison was done ("~ tissue + stage" vs "~ tissue" or "~ tissue + stage + stage:tissue" vs "~ tissue"). If it is with the former, we might have some more exploration to do. However, if it is with the latter, it is within the realm of normal performance for there to be near-zero betas. This can occur when there is a strong interaction effect, but not a strong overall effect (e.g. the expression goes up in one tissue between stages, but goes down in the other tissue between stages, effectively cancelling out on the global scale).

Hopefully this makes sense!

[apcamargo]

That number (9154) is for the "~ tissue + stage + stage:tissue" vs "~ tissue" comparison. Your explanation made perfect sense! Thank you!

As for the comparison "~ tissue + stage" vs "tissue", the minimum beta is 5.23088266008667e-05.

I also can't explain why I get low betas even when I don't have the 'tissue' variate (see plots above)

[warrenmcg]

Well, as I mentioned before, the scale from -10 to +10 makes it difficult to get a good sense of what's happening in the range of -2 to +2.

As I mentioned before, if you do the Wald test on whatever model you choose, I don't anticipate any betas being less than ~0.26-0.32 (20-25% fold change in the simple case).

In both the ~tissue + stage vs ~ tissue LRT and the ~stage vs ~1 LRT comparison, you can often have a case where the full model better predicts the data from a statistical standpoint (i.e. the improvement in the residuals, or degree of error, outweighs the reduction in degrees of freedom) but the resulting beta for the factor of interest is not "biologically meaningful". This can often result with really highly-expressed transcripts, where small changes result in large changes in the residuals when you can account for those small changes with respect to the other transcripts, but those small changes are likely not going to be meaningful when you do follow-up work (or at minimum hard to validate; can you imagine trying to validate a 5% change using qPCR? yikes...).

This distinction between statistically significant and "biologically significant" is an important one that has been discussed a lot in the literature, and it's one to keep in mind. You may get a bunch of results that are statistically significant, but you still need to prioritize those hits based on what may be a biologically meaningful change, which these days usually means what can be validated using qPCR (>20-25% if you have lots of replicates or very low intragroup variation). This is why some of the older studies looking at microarray data used a fold-change cut-off in addition to the statistical cut-off; the major reason not to use a blanket rule for fold-change cut-off is that it is still not clear what constitutes "biologically meaningful change".

[apcamargo]

Oh, I understand! Because of the large number of reads, the test gives these transcripts statistical significance, even if the fold-change is low!

I'll check the counts of those transcripts in sleuth_live.

Thank you so much! This has been truly a marvellous RNA-Seq learning experience!

[warrenmcg]

Hi @apcamargo, as I think we resolved this issue, as it pertains less to a bug and more to an interpretation (biologically meaningful change versus statistically meaningful change), I'm going to close this issue. If there are additional concerns, feel free to reopen the issue to continue the discussion.