Do the claims hold up? Displaying and analyzing paired data (CC431)
Pat reanalyzes data published in a paper recently published in the journal Nature Microbiology to see if the data support the claims when the data are analyzed using more robust statistics and visualization approaches.
- Critique video: https://youtu.be/7bA4aMENREE
- Related livestream: https://youtube.com/live/ulsljDP4fmw
- The article discussed: https://www.nature.com/articles/s41564-026-02404-5
- Sign up for a complimentary consultation to help me learn your needs: https://calendly.com/pat-riffomonas/30min
- My newsletter: https://shop.riffomonas.org/youtube
- Recorded workshops: https://riffomonas.org/workshops/
-
If you want to cite this video: https://journals.asm.org/doi/10.1128/mra.01310-22
-
packages: tidyverse, readxl, broom, ggplot2, dplyr, purrr
- functions: aes, arrange, as.numeric, download.file, filter, geom_boxplot, geom_col, geom_histogram, geom_jitter, geom_line, geom_point, ggplot, head, if_else, log, map, mutate, nchar, nest, p.adjust, pivot_longer, pivot_wider, pull, read_csv, read_excel, rename, scale_y_continuous, scale_y_log10, select, separate_wider_regex, shapiro.test, str_replace_all, t.test, tidy, unnest, and wilcox.test
This work was supported, in part, by a grant from the National Institutes of Health (R25GM166785)
library(tidyverse)
library(readxl)
library(broom)
download.file(
"https://media.springernature.com/original/springer-static/esm/art%3A10.1038%2Fs41564-026-02404-5/MediaObjects/41564_2026_2404_MOESM3_ESM.xlsx",
"fig1_data.xlsx"
)
download.file(
"https://media.springernature.com/original/springer-static/esm/art%3A10.1038%2Fs41564-026-02404-5/MediaObjects/41564_2026_2404_MOESM8_ESM.xlsx",
"figs1_data.xlsx"
)
# columns...
# * tax_level
# * full_taxonomy
# * patient
# * tissue (PT/T)
# * rel_abund
tax_data <- read_csv("Bacteria_abundance.csv") %>%
rename("full_taxonomy" = "...1") %>%
mutate(
name_length = nchar(full_taxonomy),
no_underscores = str_replace_all(full_taxonomy, "_", "") %>% nchar,
tax_level = name_length - no_underscores + 1,
tax_level = if_else(tax_level > 6, 6, tax_level),
.before = "PT1") %>%
select(-name_length, -no_underscores) %>%
pivot_longer(cols = -c(full_taxonomy, tax_level), values_to = "rel_abund") %>%
separate_wider_regex(
cols = name,
patterns = c(tissue = "P?T", patient = "\\d+"))
tax_data %>%
filter(tax_level == 5) %>%
ggplot(aes(x = paste(tissue, patient), y = rel_abund, fill = full_taxonomy)) +
geom_col(show.legend = FALSE)
bacillota_long <- tax_data %>%
filter(tax_level == 1) %>%
filter(full_taxonomy == "Bacillota")
wilcox.test(rel_abund ~ tissue, data = bacillota_long)
bacillota_wide <- bacillota_long %>%
pivot_wider(names_from = tissue, values_from = rel_abund)
wilcox.test(bacillota_wide$PT, bacillota_wide$T, paired = TRUE)
#unpaired test
tax_data %>%
filter(tax_level == 1) %>%
nest(data = -c(full_taxonomy, tax_level)) %>%
mutate(test = map(
data,
~wilcox.test(rel_abund ~ tissue, data = .x) %>% tidy())) %>%
unnest(test) %>%
mutate(p.adjusted = p.adjust(p.value, method = "BH")) %>%
arrange(p.value)
# filter(full_taxonomy == "Bacillota")
#paired test
tax_data %>%
filter(tax_level == 1) %>%
pivot_wider(names_from = tissue, values_from = rel_abund) %>%
nest(data = -c(full_taxonomy, tax_level)) %>%
mutate(test = map(
.x = data,
~wilcox.test(.x$PT, .x$T, paired = TRUE) %>% tidy())) %>%
unnest(test) %>%
select(-data, -statistic, -method, -alternative) %>%
mutate(p.adjusted = p.adjust(p.value, method = "BH")) %>%
arrange(p.value)
all_tests <- tax_data %>%
pivot_wider(names_from = tissue, values_from = rel_abund) %>%
nest(data = -c(full_taxonomy, tax_level)) %>%
mutate(test = map(
.x = data,
~wilcox.test(.x$PT, .x$T, paired = TRUE) %>% tidy())) %>%
unnest(test) %>%
select(-data, -statistic, -method, -alternative) %>%
mutate(p.adjusted = p.adjust(p.value, method = "BH"), .by = tax_level)
all_tests %>%
filter(tax_level == 6) %>%
arrange(p.value) %>%
head(n = 1) %>%
pull(full_taxonomy)
bfrag <- "Bacteroidota_Bacteroidia_Bacteroidales_Bacteroidaceae_Bacteroides_Bacteroides fragilis"
bfrag_data <- tax_data %>%
filter(full_taxonomy == bfrag)
#wrong!
bfrag_data %>%
ggplot(aes(x = tissue, y = rel_abund, fill = tissue)) +
geom_boxplot(coef = 10, staplewidth = 0.25) +
scale_y_continuous(limits = c(0, 8e-5))
bfrag_data %>%
pivot_wider(names_from = tissue, values_from = rel_abund) %>%
mutate(
first_eight = as.numeric(patient) <= 8,
increasing = T > PT) %>%
pivot_longer(cols = c(PT, T), names_to = "tissue", values_to = "rel_abund") %>%
ggplot(aes(x = tissue, y = rel_abund, group = patient,
color = increasing, linetype = first_eight)) +
geom_line()
d <- bfrag_data %>%
pivot_wider(names_from = tissue, values_from = rel_abund)
wilcox.test(d$PT, d$T, paired = TRUE) #0.0006143
d %>%
mutate(first_eight = as.numeric(patient) <= 8) %>%
nest(data = -first_eight) %>%
mutate(test = map(data, ~wilcox.test(.x$PT, .x$T, paired = TRUE) %>% tidy())) %>%
unnest(test)
long_data <- read_excel("fig1_data.xlsx", sheet = "Figure 1c", skip = 1) %>%
mutate(patient = 1:n()) %>%
select(patient, t_to_pt_ratio = T, PT) %>%
pivot_longer(-patient)
long_data %>%
ggplot(aes(x = name, y = value, group = patient)) +
geom_line() +
geom_point()
wide_data <- long_data %>%
pivot_wider(names_from = name, values_from = value)
# their analysis, wrong test
shapiro.test(wide_data$t_to_pt_ratio)
t.test(wide_data$t_to_pt_ratio, wide_data$PT, paired = TRUE)
t.test(wide_data$t_to_pt_ratio, mu = 1)
# my analysis correct test if data are transformed to be normal
shapiro.test(log(wide_data$t_to_pt_ratio))
t.test(log(wide_data$t_to_pt_ratio), mu = 0)
# correct test for non-normal data
wilcox.test(wide_data$t_to_pt_ratio, wide_data$PT, paired = TRUE)
wilcox.test(wide_data$t_to_pt_ratio, mu = 1)
wide_data %>%
ggplot(aes(x = log(t_to_pt_ratio))) + geom_histogram()
wide_data %>%
ggplot(aes(x = "Ratio", y = t_to_pt_ratio)) +
geom_boxplot(width = 0.6) +
geom_jitter(color = "red", width = 0.3) +
scale_y_log10()