-
Notifications
You must be signed in to change notification settings - Fork 3
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
LinDA for longitudinal microbiome analysis #7
Comments
For longitudinal data, one may use the mixed-effect model to incorporate subject as a random effect, since one subject always has multiple observations causing correlations among data. The use of mixed-effect models in LinDA is the same as in linear regression (lm), such as 'formula = '~Smoke+Sex+(1|SubjectID)' in the example in the readme page. |
Dear @zhouhj1994 I am using LinDA for longitudinal microbiome analysis. I have questions as a few things are not clear to me. This is also related to another issue #5 I have 3 groups (Control, Placebo, Treatment) and 3 timepoints (0, 22, 52). Thus in total there are 18 pairwise comparisons for which I want to know what taxa are differentially abundant. I have used the following code:
This code should allow me to see if there is a difference between Control and Treatment at Day 22. linda.object$covariance is NULL so I cant complete the calculation. However, even with this method its not possible to get all pairwise comparisons. For example - Treatment at Day 0 vs Treatment at Day 22 or to check if there is a difference between Control and Placebo at any time. How should I proceed to get all pairwise comparisons? Your help is greatly appreciated. Many thanks Kind regards, P |
Thanks for the use of our package. I am not sure whether your codes will work. If the "Treatment" variable has three levels: Control, Placebo, Treatment, then "Control" will be treated as the reference level (alphabetizing rule), so there shouldn't be such results as "linda.obj$output$TreatmentControl" but "linda.obj$output$TreatmentPlacebo" and "linda.obj$output$TreatmentTreatment" instead. To get all pairwise comparisons, for example,
alp1 <- linda.obj$output$Time22$log2FoldChange
se1 <- linda.obj$output$Time22$lfcSE
df1 <- linda.obj$output$Time22$df
alp2 <- linda.obj$output$TreatmentTreatment:Time22$log2FoldChange
se2 <- linda.obj$output$TreatmentTreatment:Time22$lfcSE
df2 <- linda.obj$output$TreatmentTreatment:Time22$df
cov12 <- linda.obj$covariance$Time22$TreatmentTreatment:Time22
stat <- (alp1+alp2) / (sqrt(se1^2+se2^2+2*cov12))
pvalue <- 2 * pt(-abs(stat), (df1+df2)/2)
padj <- p.adjust(pvalue, method = 'BH')
rownames(linda.obj$output[[1]])[padj <= 0.1]
# Time 0
alp1 <- linda.obj$output$TreatmentPlacebo$log2FoldChange
se1 <- linda.obj$output$TreatmentPlacebo$lfcSE
df1 <- linda.obj$output$TreatmentPlacebo$df
stat <- alp1/se1
pvalue <- 2 * pt(-abs(stat), df1)
padj <- p.adjust(pvalue, method = 'BH')
rownames(linda.obj$output[[1]])[padj <= 0.1]
# Time 22
alp2 <- linda.obj$output$TreatmentPlacebo:Time22$log2FoldChange
se2 <- linda.obj$output$TreatmentPlacebo:Time22$lfcSE
df2 <- linda.obj$output$TreatmentPlacebo:Time22$df
cov12 <- linda.obj$covariance$TreatmentPlacebo$TreatmentPlacebo:Time22
stat <- (alp1+alp2) / (sqrt(se1^2+se2^2+2*cov12))
pvalue <- 2 * pt(-abs(stat), (df1+df2)/2)
padj <- p.adjust(pvalue, method = 'BH')
rownames(linda.obj$output[[1]])[padj <= 0.1]
# Time 52: similar to Time 22
|
Thanks for the response. This is very very helpful. However, I cannot find the covariances in the output.
There is no output called covariance in linda.obj. How are you calculating this or what am I missing? Thanks very much. P |
Sorry for the late reply. Are you using the latest version of the package? The old version didn't include the covariance in the output. Have you tried this example? https://github.com/zhouhj1994/LinDA Can you see the covariances in the output? Thanks. |
No problem. Yes I was using wrong version. I was using the linda function implemented in the microbiomeStat package. I have now installed the LinDA package correctly via devtools and its al working perfectly. Thanks again for all your help. P |
Great. |
Thank you for the insightful discussion. I have a question related to this. It would be highly appreciated if you could help me with the explanation of the LinDA results. @zhouhj1994 |
If you want to understanding the effect of time, in other words, how is the abundance changed by one unit change of time, you may treat it as continuous. |
Hi @zhouhj1994, thank you for the answer. I've executed the following command: fit <- linda(feature.dat=feature.dat, meta.dat=meta.dat, formula = "~ treatment * time"), where "time" includes the values 0, 22, and 52, and "treatment" includes two levels, "yes" and "no." Could you help me check whether my use of linda output is correct? Much appreciated! If I consider "time" as a continuous variable: If I consider "time" as a categorical variable: |
Your interpretation is not bad. But you may want to say "...while keeping other variables fixed" or something like that and be careful of the interaction term. Your model is "treatment + time + treatment * time"
The effect of time: for treatment "no", the effect is "fit$output$time$log2FoldChange", while for treatment "yes", the effect is "fit$output$time$log2FoldChange" + "fit$output$treatmentyes:time$log2FoldChange" The effect of treatment: while fixing time, the treatment "yes" yields "fit$output$treatmentyes$log2FoldChange" + "fit$output$treatmentyes:time$log2FoldChange" change in comparison to the treatment "no".
The effect of time22: for treatment "no", the effect is "fit$output$time22$log2FoldChange", while for treatment "yes", the effect is "fit$output$time22$log2FoldChange" + "fit$output$treatmentyes:time22$log2FoldChange" The effect of treatment: for fixed time22, the treatment "yes" yields "fit$output$treatmentyes$log2FoldChange" + "fit$output$treatmentyes:time22$log2FoldChange" change in comparison to the treatment "no". |
Thank you for the correction and the details! @zhouhj1994 |
Hi, I have a naive question. What does |
When time is a continuous variable, fit$output$time$log2FoldChange means the effect of time on the response variable (if there are no interactions considered). We might say, a one unit increase in time corresponds to change of fit$output$time$log2FoldChange in response. |
Here, the response variable is the log2-transformed microbiome abundance. |
Hi,
I am currently using MassLin2 to analyze my longitudinal microbiome dataset; however, I saw your recent paper describing LinDA and its strengths in comparison to MassLin2. Therefore, I would like to compare my current results to those I would get from LinDA. However, after browsing through your page, I do not see an example of how to use LinDA to analyze longitudinal microbiome data. Could you point me in the direction of a tutorial or webpage that demonstrates how to adapt LinDA for longitudinal analysis?
Thanks in advance.
The text was updated successfully, but these errors were encountered: