Skip to content

Commit

Permalink
Update figures and total sum scaling (to a million)
Browse files Browse the repository at this point in the history
  • Loading branch information
@sdgamboa
sdgamboa committed Nov 19, 2024
1 parent 9d58596 commit d63644a
Show file tree
Hide file tree
Showing 7 changed files with 118 additions and 33 deletions.
2 changes: 1 addition & 1 deletion R/new_DA_methods.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -263,7 +263,7 @@ DA_lefse <- function(
statInfo$adjP <- seq(0.09, 0, length.out = nrow(statInfo))
# statInfo$adj_pval <- stats::p.adjust(statInfo$rawP, method = "fdr")
rownames(statInfo) <- statInfo[["features"]] ## lefser version 1.15
rownames(statInfo) <- statInfo[["Names"]] ## lefser version 1.14
# rownames(statInfo) <- statInfo[["Names"]] ## lefser version 1.14
colnames(statInfo) <- c("Taxa", "LDA_scores", "abs_score", "rawP", "adjP")

pValMat <- statInfo[, c("rawP", "adjP")]
Expand Down
Binary file modified vignettes/articles/Figure1.pdf
View file
Binary file not shown.
Binary file modified vignettes/articles/Figure2.pdf
View file
Binary file not shown.
Binary file modified vignettes/articles/Figure3.pdf
View file
Binary file not shown.
Binary file modified vignettes/articles/Figure4.pdf
View file
Binary file not shown.
48 changes: 44 additions & 4 deletions vignettes/articles/HMP_2012_16S_gingival_V35.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -286,6 +286,17 @@ enrichmentSummary <- purrr::map(enrichment, ~ {
TRUE ~ direction
)
)
# enrichmentSummary <- enrichmentSummary |>
# mutate(
# method = case_when(
# grepl("lefse", method) ~ sub("lefse", "LEfSe", method),
# grepl("wilcox", method) ~ sub("wilcox", "Wilcox", method),
# TRUE ~ method
# )
# )
```

Create enrichment plot
Expand All @@ -298,6 +309,21 @@ enPlot <- enrichmentSummary |>
direction, levels = c("Supragingival", "Subgingival")
)
) |>
mutate(
method = case_when(
grepl("lefse", method) ~ sub("lefse", "LEfSe", method),
grepl("wilcox", method) ~ sub("wilcox", "Wilcox", method),
TRUE ~ method
)
) |>
mutate(
annotation = case_when(
annotation == "aerobic" ~ "Aerobic",
annotation == "anaerobic" ~ "Anaerobic",
annotation == "facultative_anaerobic" ~ "Facultative anaerobic",
TRUE ~ annotation
)
) |>
ggplot(aes(method, n)) +
geom_col(
aes(fill = annotation),
Expand All @@ -315,7 +341,7 @@ enPlot <- enrichmentSummary |>
labs(
x = "DA method", y = "Number of DA taxa"
) +
# theme_minimal() +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "bottom"
Expand Down Expand Up @@ -367,10 +393,24 @@ positives <- map(1:length(DA_output), .f = function(i) {
## Plot TP - FP

```{r, fig.height = 4, fig.width=12}
# names(vec) <- positives$base_method
positives <- positives |>
mutate(diff = jitter(TP - FP, amount = 1.5, factor = 2)) |>
mutate(
base_method = case_when(
grepl("lefse", base_method) ~ sub("lefse", "LEfSe", base_method),
grepl("wilcox", base_method) ~ sub("wilcox", "Wilcox", base_method),
TRUE ~ base_method
),
method = case_when(
grepl("lefse", method) ~ sub("lefse", "LEfSe", method),
grepl("wilcox", method) ~ sub("wilcox", "Wilcox", method),
TRUE ~ method
)
)
vec <- positives$color
names(vec) <- positives$base_method
posPlot <- positives |>
mutate(diff = jitter(TP - FP, amount = 1.5, factor = 2)) |>
ggplot(aes(top, diff)) +
geom_line(
aes(
Expand All @@ -389,7 +429,7 @@ posPlot <- positives |>
scale_shape(name = "Normalization") +
scale_linetype(name = "Normalization") +
scale_color_manual(values = vec, name = "Base method") +
# theme_minimal() +
theme_minimal() +
theme(legend.position = "bottom")
```

Expand All @@ -406,7 +446,7 @@ pp
## Export the figure
ggsave(
filename = "Figure1.pdf", plot = pp, dpi = 300,
height = 9, width = 10,
height = 9, width = 11,
)
```

Expand Down
101 changes: 73 additions & 28 deletions vignettes/articles/Ravel_2011_16S_BV.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -282,6 +282,19 @@ enPlot <- enrichmentSummary |>
)
) |>
filter(annotation != "Unannotated") |>
mutate(
method = case_when(
grepl("lefse", method) ~ sub("lefse", "LEfSe", method),
grepl("wilcox", method) ~ sub("wilcox", "Wilcox", method),
TRUE ~ method
)
) |>
mutate(
annotation = case_when(
annotation == "bv-associated" ~ "BV-associated",
annotation == "hv-associated" ~ "HV-associated"
)
) |>
ggplot(aes(method, n)) +
geom_col(
aes(fill = annotation),
Expand All @@ -299,7 +312,7 @@ enPlot <- enrichmentSummary |>
labs(
x = "DA method", y = "Number of DA taxa"
) +
# theme_minimal() +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "bottom"
Expand Down Expand Up @@ -389,8 +402,22 @@ positives <- map(1:length(DA_output), .f = function(i) {
Create putative positives plot

```{r plot positves, fig.width=15, fig.height=10}
positives <- positives |>
mutate(
base_method = case_when(
grepl("lefse", base_method) ~ sub("lefse", "LEfSe", base_method),
grepl("wilcox", base_method) ~ sub("wilcox", "Wilcox", base_method),
TRUE ~ base_method
),
method = case_when(
grepl("lefse", method) ~ sub("lefse", "LEfSe", method),
grepl("wilcox", method) ~ sub("wilcox", "Wilcox", method),
TRUE ~ method
)
)
vec <- positives$color
names(vec) <- positives$base_method
posPlot <- positives |>
# mutate(diff = jitter(TP - FP, amount = 1.5, factor = 2)) |>
mutate(diff = TP - FP) |>
Expand All @@ -412,7 +439,7 @@ posPlot <- positives |>
scale_shape(name = "Normalization") +
scale_linetype(name = "Normalization") +
scale_color_manual(values = vec, name = "Base method") +
# theme_minimal() +
theme_minimal() +
theme(legend.position = "bottom")
# plots <- plot_positives(positives) |>
# map( ~ {
Expand Down Expand Up @@ -452,15 +479,15 @@ pp
# )
ggsave(
filename = "Figure2.pdf", plot = pp,
dpi = 300, height = 9, width = 10
dpi = 300, height = 9, width = 11
)
```

# Perform DA with lefse, Wilcox, and ZINQ-Cauchy manually

```{r manual transformations}
tssFun <- function(x) {
(x) / sum(x) * 100
(x) / sum(x) * 1e6
}
clrFun <- function(x) {
log(x / exp(mean(log(x))))
Expand Down Expand Up @@ -636,16 +663,27 @@ head(wilcox_raw)
Box plot of incorrect values:

```{r plot abundance values, fig.height=7}
l <- wilcox_raw |>
mutate(taxon_name = sub('genus:', '', taxon_name)) |>
{\(y) split(y, y$transformation)}()
l$`TSS + CLR` <- NULL
l$counts$value <- log(l$counts$value + 1)
l$TSS$value <- log(l$TSS$value + 1)
wilcox_raw <- reduce(l, bind_rows)
wilcox_genus_plot <- wilcox_raw |>
mutate(taxon_name = sub('genus:', '', taxon_name)) |>
mutate(
value = case_when(
transformation %in% c("counts", "TSS") ~ log(value + 1),
TRUE ~ value
)
) |>
# mutate(taxon_name = sub('genus:', '', taxon_name)) |>
# mutate(
# value = case_when(
# transformation %in% c("counts", "TSS") ~ log(value + 1),
# TRUE ~ value
# )
# ) |>
# mutate(value = log(value + 1)) |>
filter(transformation != "TSS + CLR") |>
# filter(transformation != "TSS + CLR") |>
mutate(transformation = factor(
transformation, levels = c('counts', 'TSS', 'CLR'),
labels = c('log(counts + 1)', 'log(TSS + 1)', 'CLR')
Expand All @@ -669,7 +707,8 @@ wilcox_genus_plot <- wilcox_raw |>
theme(
panel.grid.major.x = element_blank(),
legend.title = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1)
axis.text.x = element_text(angle = 45, hjust = 1, face = "italic")
)
wilcox_genus_plot
```
Expand Down Expand Up @@ -1000,13 +1039,13 @@ incorrect_taxa_lefse_clr ## the same as in wilcox.
ANCOM-BC

```{r ancombc}
ancombc <- as.data.frame(DA_output$ancombc.none$statInfo)
ancombc$taxon_name <- rownames(ancombc)
ancombc <- left_join(ancombc, taxa_annotations, by = "taxon_name") |>
relocate(taxon_name, taxon_annotation)
ancombc |>
filter(q_val <= 0.1, lfc < 0, taxon_annotation == 'bv-associated') |>
pull(taxon_name)
# ancombc <- as.data.frame(DA_output$ancombc.none$statInfo)
# ancombc$taxon_name <- rownames(ancombc)
# ancombc <- left_join(ancombc, taxa_annotations, by = "taxon_name") |>
# relocate(taxon_name, taxon_annotation)
# ancombc |>
# filter(q_val <= 0.1, lfc < 0, taxon_annotation == 'bv-associated') |>
# pull(taxon_name)
```

Expand Down Expand Up @@ -1312,11 +1351,12 @@ plot_1b <- data_with_lact |>
) +
labs(
# title = 'Relative abundace vs CLR',
title = 'Actinomyces (BV-associated)',
title = expression(italic('Actinomyces') ~ '(BV-associated)'),
x = 'log(TSS + 1)'
) +
scale_color_discrete(name = 'Condition') +
scale_size(name = 'Lactobacillus Rel. Ab.') +
# scale_size(name = 'Lactobacillus Rel. Ab.') +
scale_size(name = expression(italic('Lactobacillus') ~ 'Rel. Ab.')) +
theme_bw()
plot_2b <- data_with_lact |>
Expand All @@ -1334,11 +1374,12 @@ plot_2b <- data_with_lact |>
) +
labs(
# title = 'Relative abundace vs CLR',
title = 'Corynebacterium (BV-associated)',
title = expression(italic('Corynebacterium') ~ '(BV-associated)'),
x = 'log(TSS + 1)'
) +
scale_color_discrete(name = 'Condition') +
scale_size(name = 'Lactobacillus Rel. Ab.') +
# scale_size(name = 'Lactobacillus Rel. Ab.') +
scale_size(name = expression(italic('Lactobacillus') ~ 'Rel. Ab.')) +
theme_bw()
plot_3b <- data_with_lact |>
Expand All @@ -1356,11 +1397,12 @@ plot_3b <- data_with_lact |>
) +
labs(
# title = 'Relative abundace vs CLR',
title = 'Prevotella (BV-associated)',
title = expression(italic("Prevotella") ~ "(BV-associated)"),
x = 'log(TSS + 1)'
) +
scale_color_discrete(name = 'Condition') +
scale_size(name = 'Lactobacillus Rel. Ab.') +
# scale_size(name = 'Lactobacillus Rel. Ab.') +
scale_size(name = expression(italic('Lactobacillus') ~ 'Rel. Ab.')) +
theme_bw()
plot_4b <- data_with_lact |>
Expand All @@ -1378,11 +1420,13 @@ plot_4b <- data_with_lact |>
) +
labs(
# title = 'Relative abundace vs CLR',
title = 'Lactobacillus (HV)',
# title = 'Lactobacillus (HV-associated)',
title = expression(italic("Lactobacillus") ~ "(HV-associated)"),
x = 'log(TSS + 1)'
) +
scale_color_discrete(name = 'Condition') +
scale_size(name = 'Lactobacillus Rel. Ab.') +
# scale_size(name = 'Lactobacillus Rel. Ab.') +
scale_size(name = expression(italic('Lactobacillus') ~ 'Rel. Ab.')) +
theme_bw()
plotsb <- ggpubr::ggarrange(
Expand All @@ -1392,6 +1436,7 @@ plotsb <- ggpubr::ggarrange(
)
```


```{r, fig.height=6, fig.width=8}
plotsb
Expand Down

0 comments on commit d63644a

Please sign in to comment.