You have your abundance table formatted as what I think are supposed to be percentages. You can convert them to proportions in [0,1] by dividing out the sum of each column like so:
Taxa = Taxa %>% apply(2, function(x) x / sum(x))
You also seem to have duplicate study identifiers. Below I deduplicate the identifiers in the metadata the same way read_tsv does it when it reads in the abundances, but you should make sure that this is appropriate.
META$ID[duplicated(META$ID)] = paste0(META$ID[duplicated(META$ID)], '_1')
META = META %>% as.data.frame() %>% column_to_rownames("ID")
fit_adjust_batch runs after these two steps for me.
Dear Andrew, thank you very much for your response. I used the correction and it worked!
A follow up question, I used this code since I want to see differences between response, considering the study and controlling for age & sex.
fit_adjust_batch <- adjust_batch(feature_abd = TaXa,
batch = "study",
covariates = c("age", "sex"),
data = MeTa,
control = list(verbose = FALSE))
TaXa <- TaXa %>% column_to_rownames("ID") %>% as.matrix()
MeTa <- MeTa %>% column_to_rownames("ID")
TaXa = TaXa %>% apply(2, function(x) x / sum(x))
CRC_abd_adj <- fit_adjust_batch$feature_abd_adj # I just recycled this
fit_lm_meta <- lm_meta(feature_abd = TaXa,
batch = "study",
exposure = "response",
data = MeTa,
control = list(verbose = FALSE))
meta_fits <- fit_lm_meta$meta_fits
Could you expand on why the results in the response_R_forest.pdf (23 species) differ from the bar plot (with just 2 features).
Inspecting the meta_fits object I found 40 features/species with p-value <0.05 (and just 2 with FDR<0.05). Please find attached these results.
For a bug to get an entry in the forest plot it needs to show up in at least two studies. So you have 40 bugs with a nominal p-value below 0.05, 23 of which show up in at least two studies, and two of which have a q-value below 0.05.