1
$\begingroup$

I would like to know how to interpret the output of the formula although i been thorough loads of literature but I'm not confident yet.

So this is my full model I'm running

dds <- DESeqDataSetFromMatrix(countData=rpkm_ordered, colData=coldata, design= ~ FAB + Sex + Age + TMB + WBC + BM_percentage)

So if I have to interpret the output of this how it would be ?

What I have read is this

Here in the above formula we are trying or controlling the FAB variable while we look for variation that is coming from other variables?

Now my confusion is when I made this dds object the FAB was the variable which I fixed while I look for variation coming from other variable.

Next thing is I would like to run a reduced model.

I tried this

test.Age <- as_tibble(results(DESeq(dds, test="LRT", reduced= ~Sex + Age )))

Now my question is since I'm using the same dds object here where the FAB is the variable i fixed and with new formula which is

reduced= ~Sex + Age this so how would i interpret this model?

Any suggestion or help would be really appreciated

Improve this question
$\endgroup$
3
  • $\begingroup$ You know that RPKM is inappropriate for DESeq2? $\endgroup$
    – swbarnes2
    Commented May 30, 2022 at 17:03
  • $\begingroup$ yes..are you rpkm_ordered telling about this i suppose its just dataframe i named in my code which i didn't change.. $\endgroup$
    – kcm
    Commented May 30, 2022 at 18:15
  • $\begingroup$ @swbarnes2 sorry for this rpkm_ordered confusion since once i had to reorder my dataframe...which I did quite long time back so its just in my code.. $\endgroup$
    – kcm
    Commented May 30, 2022 at 18:16

1 Answer 1

Reset to default
2
$\begingroup$

By fitting the full model ~FAB + Sex + Age + TMB + WBC + BM_percentage you can estimate the effect of FAB (or any of the other variables) on RNA expression after adjusting for differences caused by Sex, Age, TMB, WBC, and BM percent. See this post for exactly what it means to adjust for other variables

When you fit a reduced model (in your case ~Sex + Age) and then conduct a likelihood ratio test (LRT), you are comparing the fit of the full model (~FAB + Sex + Age + TMB + WBC + BM_percentage) to the fit of the reduced model (~Sex + Age). If the fit is substantially better, the likelihood of the full model will be larger than the reduced model and the p-value will be significant. This means that adding the extra variables that are in the full model (FAB, TMB, WBC, BM_percentage) improved the fit of the model to the data because together these variables have some type of relationship with RNA expression even after adjusting for the variables in the reduced model.

The challenge with using a LRT is that it doesn't tell you which of the variables added in the full model are responsible for the improved fit - it just indicates that the full model is better. This type of test is useful when you are trying to make broad statements about several variables having a relationship with RNA expression. From the DESeq2 vignette:

The LRT is therefore useful for testing multiple terms at once, for example testing 3 or more levels of a factor at once, or all interactions between two variables. The LRT for count data is conceptually similar to an analysis of variance (ANOVA) calculation in linear regression, except that in the case of the Negative Binomial GLM, we use an analysis of deviance (ANODEV), where the deviance captures the difference in likelihood between a full and a reduced model.

Definitely read the DESeq2 vignette for more details. In your case if the LR test is significant, it means that FAB, TMB, WBC, and BM percent combined are associated with RNA expression.

Most people are interested in specific comparisons when performing differential expression analyses with DESeq2. These comparisons are usually like do males have higher gene expression than females (comparing sex differences), do patients with higher TMB have greater expression than patients with lower TMB (comparing effect of TMB on RNA expression). If you want to do a specific comparison, you will likely want to use the Wald test, which calculates p-values for individual coefficients (male vs female coefficient, TMB coefficient etc). The Wald test determines whether a coefficient is equal to 0 or some other null hypothesis (see lfcThreshold parameter).

The Wald statistic and accompanying p-values are reported by default when using results() and not specifying test="LRT"

There are some additional differences between Wald tests vs LRT, including accuracy differences when operating on small sample sizes and computational efficiency issues when operating with large sample sizes. Discussion on Wald vs LRT

TL;DR Use LRT to test if several variables or factor levels are associated with RNA expression. Use the Wald test to determine if specific comparisons are significant.

Improve this answer
$\endgroup$
8
  • $\begingroup$ thank you for the detailed write up...this is more than enough for my data reporting which would be now be more statistical rather than using loose terms thanks a ton $\endgroup$
    – kcm
    Commented Jun 2, 2022 at 21:42
  • $\begingroup$ Little doubt "The design formula should have all of the factors in your metadata that account for major sources of variation in your data. The last factor entered in the formula should be the condition of interest." are these saying the same thing you described? hbctraining.github.io/DGE_workshop_salmon_online/lessons/… $\endgroup$
    – kcm
    Commented Jun 15, 2022 at 6:43
  • 1
    $\begingroup$ The last factor in the formula is the variable whose comparison (if sex is the last variable then male vs female) will show up by default when you call results(). You can change which comparison you get with the contrast argument. $\endgroup$
    – Tomas Bencomo
    Commented Jun 16, 2022 at 18:10
  • $\begingroup$ if I have a single factor in my design for example age and i ran deseq2 what would be my interpretation would I be saying here Im testing the effect of age on the differences in gene expression?is that statistically correct statement? $\endgroup$
    – kcm
    Commented Jun 24, 2022 at 1:09
  • 1
    $\begingroup$ Yes that is correct. If the age variable is continuous, the age coefficient for a gene $g$ would represent the log2FoldChange of $g$ when age increases by 1 unit. If age has no effect on gene expression for $g$, the coefficient will be very small and not significant (assuming adequate power) $\endgroup$
    – Tomas Bencomo
    Commented Jun 25, 2022 at 6:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.