Is relative abundance an acceptable format for MicrobiomeStat data object? #75
Replies: 2 comments 3 replies
-
|
Hello Max, Thank you for your interest in MicrobiomeStat! I'm happy to address your question about using relative abundance data with CLR transformation in the package. Yes, relative abundance data is an acceptable input format for MicrobiomeStat. Your approach of:
is appropriate for many of our package functions. However, I see that you're encountering an error when running This error occurs because the one-click report function attempts to perform rarefaction on your data, which requires integer count data. Since you're using relative abundance data that has already been CLR-transformed, rarefaction is not appropriate (nor mathematically possible) for your dataset. To resolve this issue, you need to set Your modified code should look like: mStat_generate_report_long(
data.obj = data.obj,
group.var = group.var,
test.adj.vars = test.adj.vars,
vis.adj.vars = vis.adj.vars,
strata.var = strata.var,
subject.var = subject.var,
time.var = time.var,
alpha.obj = NULL,
alpha.name = c("shannon", "observed_species"),
dist.obj = NULL,
dist.name = c("BC", "Jaccard"),
feature.dat.type = "other",
feature.analysis.rarafy = FALSE, # Set this to FALSE for CLR-transformed data
vis.feature.level = c("Phylum", "Family", "Genus"),
test.feature.level = c("Family"),
bar.area.feature.no = 30,
heatmap.feature.no = 40,
feature.mt.method = "fdr",
feature.sig.level = 0.1,
theme.choice = "classic",
output.file = "~/my/directory/of/interest"
)For alpha diversity calculations, you may want to provide pre-calculated alpha diversity indices using the The CLR transformation is excellent for differential abundance and correlation analyses, but some analyses in the package are designed with count data in mind. By setting the appropriate parameters as suggested above, you should be able to use your CLR-transformed relative abundance data effectively with our package. If you have any further questions or encounter other issues, please don't hesitate to ask! Best regards, |
Beta Was this translation helpful? Give feedback.
-
|
Hi Chen, Thank you so much for the prompt response. My only trouble is that I already have feature.analysis.rarafy = FALSE, as well as feature.dat.type = "other" and I still get the same error. Best, Max |
Beta Was this translation helpful? Give feedback.
-
|
Hello Max, Thank you for your follow-up. I apologize for the confusion. After examining the code more thoroughly, I understand why you're still encountering the error despite setting The IssueThe problem is that Specifically, in the rarefy.data.obj <- mStat_normalize_data(data.obj = data.obj, method = 'Rarefy', depth = depth)$data.obj.normThis rarefaction step is performed even when SolutionFor your CLR-transformed relative abundance data, you need to provide pre-calculated alpha and beta diversity objects to bypass the automatic rarefaction step. Here's how to modify your code: # First, calculate alpha diversity on your original (non-CLR) relative abundance data
# If you still have access to the pre-CLR transformed data:
alpha.obj <- mStat_calculate_alpha_diversity(
x = physeq_MetaG_filtered@otu_table, # Use the filtered but not CLR-transformed data
alpha.name = c("shannon", "observed_species")
)
# Calculate beta diversity on your original (non-CLR) relative abundance data
dist.obj <- mStat_calculate_beta_diversity(
data.obj = physeq_MetaG_MbmStat, # Your MicrobiomeStat data object
dist.name = c("BC", "Jaccard")
)
# Now run the report with pre-calculated objects
mStat_generate_report_long(
data.obj = physeq_MetaG_MbmStat,
group.var = group.var,
test.adj.vars = test.adj.vars,
vis.adj.vars = vis.adj.vars,
strata.var = strata.var,
subject.var = subject.var,
time.var = time.var,
alpha.obj = alpha.obj, # Pre-calculated alpha diversity
alpha.name = c("shannon", "observed_species"),
dist.obj = dist.obj, # Pre-calculated beta diversity
dist.name = c("BC", "Jaccard"),
feature.dat.type = "other",
feature.analysis.rarafy = FALSE,
vis.feature.level = c("Phylum", "Family", "Genus"),
test.feature.level = c("Family"),
bar.area.feature.no = 30,
heatmap.feature.no = 40,
feature.mt.method = "fdr",
feature.sig.level = 0.1,
theme.choice = "classic",
output.file = "~/my/directory/of/interest"
)If you no longer have access to the pre-CLR transformed data, you might need to back-transform your CLR data (by exponentiating and then normalizing to sum to 1) before calculating alpha diversity. Alternative ApproachIf the above solution doesn't work for you, an alternative approach would be to use individual analysis functions instead of the one-click report. This way, you can have more control over each step of the analysis and avoid the automatic rarefaction. For example, you could use functions like:
Each of these can be run with your CLR-transformed data directly, without going through the report generation pipeline. I hope this helps! Let me know if you have any further questions. Best regards, |
Beta Was this translation helpful? Give feedback.
-
|
Thanks, Chen. Following up on this issue, I have two questions. I've switched to applying single functions, with mixed success. Yesterday, I ran the generate_taxa_trend_test_long/generate_taxa_trend_volcano_long functions with my physeq_MetaG (non transformed object created above) and had no difficulties. When I try with the CLR transformed version (Mcbm_Stat_physeq_CLR), I get: When I inspect my object, there are no zeros, as I have performed: Are you aware of what might cause this? Can I use the CLR transformed objects, or should I stick with the non-transformed data, and if so is that reliable? Do you know why, with both CLR transformed and non-transformed objects, I cannot perform the longitudinal dotplot nor the area plot functions? I get this error when using the CLR object: and this when using the non-transformed object: |
Beta Was this translation helpful? Give feedback.
-
|
Hello Max, Thank you for your follow-up questions. I've investigated the issues you're encountering with CLR-transformed data in MicrobiomeStat. Let me address each of your concerns: Issue 1: NaNs and "0 (non-NA) cases" Error with CLR-transformed DataThe error you're encountering with your CLR-transformed object ( The Problem:
Solution:When using CLR-transformed data, you need to set Here's a workaround: # For your CLR-transformed object
# First, make a copy of your CLR data
Mcbm_Stat_physeq_CLR_modified <- Mcbm_Stat_physeq_CLR
# Apply exp(x) to convert back from CLR to a positive space
# This counteracts the log2 transformation that will be applied later
otu_table(Mcbm_Stat_physeq_CLR_modified) <- 2^as(otu_table(Mcbm_Stat_physeq_CLR_modified), "matrix")
# Now use this modified object with the functions
test_results <- generate_taxa_trend_test_long(
data.obj = Mcbm_Stat_physeq_CLR_modified,
subject.var = "CaseID",
time.var = "Relative_timepoint",
group.var = "Celiac_Disease",
adj.vars = test.adj.vars,
feature.level = c("Phylum"),
feature.dat.type = "other" # Important to set this
)Issue 2: Problems with Longitudinal Dotplot and Area Plot FunctionsThe error with the non-transformed object in Solution:
DotPlot_res <- generate_taxa_per_time_test_long(
data.obj = physeq_MetaG_MbmStat,
subject.var = "CaseID",
time.var = "Relative_timepoint",
group.var = "Celiac_Disease",
adj.vars = test.adj.vars,
prev.filter = 0.05, # Reduced from 0.1
abund.filter = 0.0005, # Reduced from 0.001
feature.level = c("Phylum"),
feature.dat.type = "proportion"
)
# Check sample counts per time point
table(physeq_MetaG_MbmStat$meta.dat$Relative_timepoint,
physeq_MetaG_MbmStat$meta.dat$Celiac_Disease)
General Recommendations for Working with CLR-transformed Data in MicrobiomeStat
Let me know if you need any clarification or have additional questions! Best regards, |
Beta Was this translation helpful? Give feedback.
-
Hi there!
I am really excited by this package, thank you for putting it together. I have a longitudinal shotgun metagenomics study, though I only have relative abundances (output from HUMAnN/MetaPHLAN pipeline from the Huttenhower lab) and extensive metadata/exposomics. I've transformed these data into a phyloseq object and used the "mStat_convert_phyloseq_to_data_obj" function to create a data obj, but I want to make sure that this data (which has been clr transformed) is an acceptable input.
Phyloseq code:
Finally:
Gives what looks to be an acceptable output. Yet, when I run the one-click analysis:
I get the error:
EDIT: Which, to be clear, appears to confirm my initial concerns about using relative abundances. Despite feature.dat.type = "other",
& feature.analysis.rarafy = FALSE (I also tried "proportional"), the pipeline tries to use vegan::rrarefy().
Beta Was this translation helpful? Give feedback.
All reactions