Say we wanted to compare metabolite profiles between our Coculture
and ANG18 groups. We can extract group means from the filtered compounds
using the mpactr function get_group_averages(), isolate the
groups of interest, here Coculture and ANG18 monocoluture, and calculate
compound fold change.
Calculate fold change
get_group_averages(data_filtered) %>%
filter(Biological_Group == "Coculture" |
Biological_Group == "ANG18") %>%
select(Compound, Biological_Group, average) %>%
pivot_wider(names_from = Biological_Group, values_from = average) %>%
mutate(fc = Coculture / ANG18) %>%
head()
#> # A tibble: 6 × 4
#> Compound ANG18 Coculture fc
#> <chr> <dbl> <dbl> <dbl>
#> 1 1000 0 0 NaN
#> 2 1001 0 0 NaN
#> 3 1002 0 0 NaN
#> 4 1003 0 0 NaN
#> 5 1004 0 0 NaN
#> 6 1005 0 0 NaN
Inherently, tandem MS/MS datasets can be filled with many zeros. In
some instances a compound is found in the experimental group but not in
the control group, or vice versa. In these cases, the fold change
calculation yields an infinite number as there is a zero either in the
numerator or denominator (\(fold chage =
experimental / control\)). There is also the chance that the
compound is not found in either group, yielding a fold change on NaN
(0/0).
Compounds that are not in either group are of no interest in this
comparison and can therefore be removed from the analysis.
get_group_averages(data_filtered) %>%
filter(Biological_Group == "Coculture" |
Biological_Group == "ANG18") %>%
select(Compound, Biological_Group, average) %>%
pivot_wider(names_from = Biological_Group, values_from = average) %>%
mutate(nonzero_compound = if_else(Coculture == 0 & ANG18 == 0,
FALSE,
TRUE)) %>%
filter(nonzero_compound == TRUE) %>%
mutate(fc = Coculture / ANG18) %>%
head()
#> # A tibble: 6 × 5
#> Compound ANG18 Coculture nonzero_compound fc
#> <chr> <dbl> <dbl> <lgl> <dbl>
#> 1 1007 0 8929. TRUE Inf
#> 2 1024 0 11317. TRUE Inf
#> 3 103 0 5152. TRUE Inf
#> 4 1039 1145. 0 TRUE 0
#> 5 1040 669. 0 TRUE 0
#> 6 105 0 16501. TRUE Inf
There are multiple ways to handle compounds that are found in the one
group but not the other. Here we will create pseudo-counts by adding a
small number (0.001) to all counts prior to calculating fold change.
This will alleviate division by 0.
First, extract group averages from the mpactr object and
select the two groups of interest. Here we are interested in the columns
Coculture and ANG18:
get_group_averages(data_filtered) %>%
filter(Biological_Group == "Coculture" |
Biological_Group == "ANG18") %>%
select(Compound, Biological_Group, average) %>%
head()
#> Compound Biological_Group average
#> <char> <char> <num>
#> 1: 1000 ANG18 0
#> 2: 1000 Coculture 0
#> 3: 1001 ANG18 0
#> 4: 1001 Coculture 0
#> 5: 1002 ANG18 0
#> 6: 1002 Coculture 0
To calculate fold change, we need to have two columns: one for the
control group, here ANG18, and one for the experimental group,
Coculture:
get_group_averages(data_filtered) %>%
filter(Biological_Group == "Coculture" |
Biological_Group == "ANG18") %>%
select(Compound, Biological_Group, average) %>%
head()
#> Compound Biological_Group average
#> <char> <char> <num>
#> 1: 1000 ANG18 0
#> 2: 1000 Coculture 0
#> 3: 1001 ANG18 0
#> 4: 1001 Coculture 0
#> 5: 1002 ANG18 0
#> 6: 1002 Coculture 0
Now we create pseudo-counts by adding 0.001 to the counts column:
get_group_averages(data_filtered) %>%
filter(Biological_Group == "Coculture" |
Biological_Group == "ANG18") %>%
select(Compound, Biological_Group, average) %>%
mutate(average = average + 0.001) %>%
pivot_wider(names_from = Biological_Group, values_from = average) %>%
head()
#> # A tibble: 6 × 3
#> Compound ANG18 Coculture
#> <chr> <dbl> <dbl>
#> 1 1000 0.001 0.001
#> 2 1001 0.001 0.001
#> 3 1002 0.001 0.001
#> 4 1003 0.001 0.001
#> 5 1004 0.001 0.001
#> 6 1005 0.001 0.001
Now we can remove compounds with an intensity of 0 in both groups.
This equates to a pseuo-count of 0.001:
get_group_averages(data_filtered) %>%
filter(Biological_Group == "Coculture" |
Biological_Group == "ANG18") %>%
select(Compound, Biological_Group, average) %>%
mutate(average = average + 0.001) %>%
pivot_wider(names_from = Biological_Group, values_from = average) %>%
mutate(nonzero_compound = if_else(Coculture == 0.001 & ANG18 == 0.001,
FALSE,
TRUE)) %>%
filter(nonzero_compound == TRUE) %>%
head()
#> # A tibble: 6 × 4
#> Compound ANG18 Coculture nonzero_compound
#> <chr> <dbl> <dbl> <lgl>
#> 1 1007 0.001 8929. TRUE
#> 2 1024 0.001 11317. TRUE
#> 3 103 0.001 5152. TRUE
#> 4 1039 1145. 0.001 TRUE
#> 5 1040 669. 0.001 TRUE
#> 6 105 0.001 16501. TRUE
Next, we calculate fold change for all remaining compounds:
get_group_averages(data_filtered) %>%
filter(Biological_Group == "Coculture" |
Biological_Group == "ANG18") %>%
select(Compound, Biological_Group, average) %>%
mutate(average = average + 0.001) %>%
pivot_wider(names_from = Biological_Group, values_from = average) %>%
mutate(nonzero_compound = if_else(Coculture == 0.001 & ANG18 == 0.001,
FALSE,
TRUE)) %>%
filter(nonzero_compound == TRUE) %>%
mutate(fc = Coculture / ANG18) %>%
head()
#> # A tibble: 6 × 5
#> Compound ANG18 Coculture nonzero_compound fc
#> <chr> <dbl> <dbl> <lgl> <dbl>
#> 1 1007 0.001 8929. TRUE 8.93e+6
#> 2 1024 0.001 11317. TRUE 1.13e+7
#> 3 103 0.001 5152. TRUE 5.15e+6
#> 4 1039 1145. 0.001 TRUE 8.73e-7
#> 5 1040 669. 0.001 TRUE 1.49e-6
#> 6 105 0.001 16501. TRUE 1.65e+7
Finally, transform fold change to log2:
foldchanges <- get_group_averages(data_filtered) %>%
filter(Biological_Group == "Coculture" |
Biological_Group == "ANG18") %>%
select(Compound, Biological_Group, average) %>%
mutate(average = average + 0.001) %>%
pivot_wider(names_from = Biological_Group, values_from = average) %>%
mutate(nonzero_compound = if_else(Coculture == 0.001 & ANG18 == 0.001,
FALSE,
TRUE)) %>%
filter(nonzero_compound == TRUE) %>%
mutate(fc = Coculture / ANG18,
logfc = log2(fc))
head(foldchanges)
#> # A tibble: 6 × 6
#> Compound ANG18 Coculture nonzero_compound fc logfc
#> <chr> <dbl> <dbl> <lgl> <dbl> <dbl>
#> 1 1007 0.001 8929. TRUE 8.93e+6 23.1
#> 2 1024 0.001 11317. TRUE 1.13e+7 23.4
#> 3 103 0.001 5152. TRUE 5.15e+6 22.3
#> 4 1039 1145. 0.001 TRUE 8.73e-7 -20.1
#> 5 1040 669. 0.001 TRUE 1.49e-6 -19.4
#> 6 105 0.001 16501. TRUE 1.65e+7 24.0
Test for compounds with significant fold changes
For the t-tests, we have already calculated variation for compounds
within biological groups as well as within technical replicates in
biological groups. We can calculate the t statistic and degrees of
freedom and use the t distribution (pt()) to calculate the
p-value.
combine biological and technical variation and calculate sample
size
# Satterwaite
calc_samplesize_ws <- function(sd1, n1, sd2, n2) {
s1 <- sd1 / (n1^0.5)
s2 <- sd2 / (n2^0.5)
n <- (s1^2 / n1 + s2^2 / n2)^2
d1 <- s1^4 / ((n1^2) * (n1 - 1))
d2 <- s2^4 / ((n2^2) * (n2 - 1))
d1[which(!is.finite(d1))] <- 0
d2[which(!is.finite(d2))] <- 0
d <- d1 + d2
return(n / d)
}
my_comp <- c("Coculture", "ANG18")
stats <- get_group_averages(data_filtered) %>%
mutate(
combRSD = sqrt(techRSD^2 + BiolRSD^2),
combASD = combRSD * average,
combASD = if_else(is.na(combASD), 0, combASD),
BiolRSD = if_else(is.na(BiolRSD), 0, BiolRSD),
techRSD = if_else(is.na(techRSD), 0, techRSD),
neff = calc_samplesize_ws(techRSD, techn, BiolRSD, Bioln) + 1
) %>%
filter(Biological_Group %in% my_comp)
head(stats)
#> Compound Biological_Group average BiolRSD Bioln techRSD techn combRSD
#> <char> <char> <num> <num> <int> <num> <num> <num>
#> 1: 1000 ANG18 0 0 3 0 3 NA
#> 2: 1000 Coculture 0 0 3 0 3 NA
#> 3: 1001 ANG18 0 0 3 0 3 NA
#> 4: 1001 Coculture 0 0 3 0 3 NA
#> 5: 1002 ANG18 0 0 3 0 3 NA
#> 6: 1002 Coculture 0 0 3 0 3 NA
#> combASD neff
#> <num> <num>
#> 1: 0 NaN
#> 2: 0 NaN
#> 3: 0 NaN
#> 4: 0 NaN
#> 5: 0 NaN
#> 6: 0 NaN
calculate the t statistic
denom <- stats %>%
summarise(den = combASD^2 / (neff),
.by = c("Compound", "Biological_Group")) %>%
mutate(den = if_else(!is.finite(den), 0, den)) %>%
summarise(denom = sqrt(sum(den)), .by = c("Compound"))
t_test <- stats %>%
select(Compound, Biological_Group, average) %>%
pivot_wider(names_from = Biological_Group, values_from = average) %>%
mutate(numerator = (Coculture - ANG18)) %>% # experimental - control
left_join(denom, by = "Compound") %>%
mutate(t = abs(numerator / denom))
head(t_test)
#> # A tibble: 6 × 6
#> Compound ANG18 Coculture numerator denom t
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1000 0 0 0 0 NaN
#> 2 1001 0 0 0 0 NaN
#> 3 1002 0 0 0 0 NaN
#> 4 1003 0 0 0 0 NaN
#> 5 1004 0 0 0 0 NaN
#> 6 1005 0 0 0 0 NaN
calculate degrees of freedom
df <- stats %>%
select(Compound, Biological_Group, neff) %>%
mutate(neff = if_else(!is.finite(neff), 0, neff)) %>%
pivot_wider(names_from = Biological_Group, values_from = neff) %>%
mutate(deg = Coculture + ANG18 - 2) %>%
select(Compound, deg)
head(df)
#> # A tibble: 6 × 2
#> Compound deg
#> <chr> <dbl>
#> 1 1000 -2
#> 2 1001 -2
#> 3 1002 -2
#> 4 1003 -2
#> 5 1004 -2
#> 6 1005 -2
calculate p-value
t <- t_test %>%
left_join(df, by = "Compound") %>%
mutate(
p = (1 - pt(t, deg)) * 2,
logp = log10(p),
neg_logp = -logp
) %>%
select(Compound, t, deg, p, logp, neg_logp)
head(t)
#> # A tibble: 6 × 6
#> Compound t deg p logp neg_logp
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1000 NaN -2 NaN NaN NaN
#> 2 1001 NaN -2 NaN NaN NaN
#> 3 1002 NaN -2 NaN NaN NaN
#> 4 1003 NaN -2 NaN NaN NaN
#> 5 1004 NaN -2 NaN NaN NaN
#> 6 1005 NaN -2 NaN NaN NaN
combine t-test results with fold changes we calculated
num_ions <- t %>%
filter(p <= 1) %>%
count() %>%
pull()
fc <- foldchanges %>%
left_join(t, by = "Compound") %>%
arrange(p) %>%
mutate(
qval = seq_len(length(p)),
qval = p * num_ions / qval
) %>%
arrange(desc(p))
min_qval <- 1
for (i in seq_len(nrow(fc))) {
if (!is.finite(fc$qval[i])) {
next
}
if (fc$qval[i] < min_qval) {
min_qval <- fc$qval[i]
} else {
fc$qval[i] <- min_qval
}
}
fc2 <- fc %>%
mutate(neg_logq = -log10(qval))
Volcano plot
Now we can create a volcano plot to see the magnitude of metabolite
abundance differences between the Cocoulture and ANG18 groups.
We can plot log2 fold change values against log10 p-values
or log10 q-values, both of which we calculated in our t-tests
above.
log2 fold change ~ log10 p-values
The base of a volcano plot is a scatter plot with log2 fold changes
on the x axis and -log10 p-values on the y axis:
fc2 %>%
ggplot() +
aes(x = logfc, y = neg_logp) +
geom_point() +
labs(x = "log2 Fold Change",
y = "-Log~10~ *P*",
color = "Differential Abundance") +
theme_bw()
This is good, but we probably want to know which compounds are above
the p-value threshold and outside the log2 fold change threshold. Fold
change threshold commonly ranges from 1.5 - 2. For this plot, we will
show a p-value threshold of 0.05 and fold change threshold of 1.5.
First add a horizontal line to denote the cutoff of -log10
0.05:
fc2 %>%
ggplot() +
aes(x = logfc, y = neg_logp) +
geom_point() +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
labs(x = "log2 Fold Change",
y = "-Log~10~ *P*",
color = "Differential Abundance") +
theme_bw()
Now we can add vertical lines showing the positive and negative
cutoffs for log2 fold change:
fc2 %>%
ggplot() +
aes(x = logfc, y = neg_logp) +
geom_point() +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
labs(x = "log2 Fold Change",
y = "-Log~10~ *P*",
color = "Differential Abundance") +
geom_vline(xintercept = log2(1.5), linetype = "dashed") +
geom_vline(xintercept = -log2(1.5), linetype = "dashed") +
theme_bw()
So now we can see the thresholds, but it’s difficult to make out
which compounds have significant fold change values. It would be nice to
color the compounds by their significance. Specifically, compounds whose
intensity are significantly higher in the experimental group, compounds
whose intensity are significantly lower in the experimental group,
compounds with significant p-values but the fold change is below the
fold change threshold, and compounds whose fold change are not
significant.
We can do this by adding a new variable to our dataset and using this
variable to color our points.
First, add a new variable sig which has values
Increased (fold change > 1.5 & pvalue < 0.05),
Decreased (fold change < -1.5 & pvalue < 0.05),
Inconclusive (fold change -1.5 - 1.5 & pvalue <
0.05), and not_sig (pvalue > 0.05).
fc2 %>%
mutate(
sig = case_when(
p > 0.05 ~ "not_sig",
p <= 0.05 & logfc > log2(1.5) ~ "Increased",
p <= 0.05 & logfc < -log2(1.5) ~ "Decreased",
TRUE ~ "Inconclusive"
),
sig = factor(sig,
levels = c("Increased", "Decreased", "Inconclusive", "not_sig"),
labels = c("Increased", "Decreased", "Inconclusive", "Not significant")
)
) %>%
select(Compound, ANG18, Coculture, fc, logfc, p, sig) %>%
head()
#> # A tibble: 6 × 7
#> Compound ANG18 Coculture fc logfc p sig
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
#> 1 637 2174. 2022. 9.30e-1 -0.105 0.458 Not significant
#> 2 1284 0.001 312. 3.12e+5 18.3 0.429 Not significant
#> 3 856 130. 0.001 7.69e-6 -17.0 0.429 Not significant
#> 4 1127 5164. 0.001 1.94e-7 -22.3 0.169 Not significant
#> 5 667 0.001 619. 6.19e+5 19.2 0.167 Not significant
#> 6 1302 0.001 2900. 2.90e+6 21.5 0.166 Not significant
Now we can set the color aesthetic to the variable sig
and select custom colors with scale_color_manual():
fc2 %>%
mutate(
sig = case_when(
p > 0.05 ~ "not_sig",
p <= 0.05 & logfc > log2(1.5) ~ "Increased",
p <= 0.05 & logfc < -log2(1.5) ~ "Decreased",
TRUE ~ "Inconclusive"
),
sig = factor(sig,
levels = c("Increased", "Decreased", "Inconclusive", "not_sig"),
labels = c("Increased", "Decreased", "Inconclusive", "Not significant")
)
) %>%
ggplot() +
aes(x = logfc, y = neg_logp, color = sig) +
geom_point() +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
geom_vline(xintercept = log2(1.5), linetype = "dashed") +
geom_vline(xintercept = -log2(1.5), linetype = "dashed") +
scale_color_manual(values = c("red", "blue", "grey", "black")) +
labs(
x = "log2 Fold Change",
y = "-Log~10~ *P*",
color = "Differential Abundance"
) +
theme_bw()
Finally, adjust axis text and titles to your liking:
fc2 %>%
mutate(
sig = case_when(
p > 0.05 ~ "not_sig",
p <= 0.05 & logfc > log2(1.5) ~ "Increased",
p <= 0.05 & logfc < -log2(1.5) ~ "Decreased",
TRUE ~ "Inconclusive"
),
sig = factor(sig,
levels = c("Increased", "Decreased", "Inconclusive", "not_sig"),
labels = c("Increased", "Decreased", "Inconclusive", "Not significant")
)
) %>%
ggplot() +
aes(x = logfc, y = neg_logp, color = sig) +
geom_point() +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
geom_vline(xintercept = log2(1.5), linetype = "dashed") +
geom_vline(xintercept = -log2(1.5), linetype = "dashed") +
scale_color_manual(values = c("red", "blue", "grey", "black")) +
labs(
x = "log2 Fold Change",
y = "-Log~10~ *P*",
color = "Differential Abundance"
) +
theme_bw() +
theme(
axis.title = element_markdown(size = 20),
axis.text = element_text(size = 15)
)
We can also make the plot interactive, showing the compound id for
each point:
volcano <- fc2 %>%
mutate(
sig = case_when(
p > 0.05 ~ "not_sig",
p <= 0.05 & logfc > log2(1.5) ~ "Increased",
p <= 0.05 & logfc < -log2(1.5) ~ "Decreased",
TRUE ~ "Inconclusive"
),
sig = factor(sig,
levels = c("Increased", "Decreased", "Inconclusive", "not_sig"),
labels = c("Increased", "Decreased", "Inconclusive", "Not significant")
)
) %>%
ggplot() +
aes(
x = logfc, y = neg_logp, color = sig,
text = paste0("Compound: ", Compound)
) +
geom_point() +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
geom_vline(xintercept = log2(1.5), linetype = "dashed") +
geom_vline(xintercept = -log2(1.5), linetype = "dashed") +
scale_color_manual(values = c("red", "blue", "grey", "black")) +
labs(
x = "log2 Fold Change",
y = "-Log~10~ *P*",
color = "Differential Abundance"
) +
theme_bw() +
theme(
axis.title = element_markdown(size = 20),
axis.text = element_text(size = 15)
)
ggplotly(volcano)
log2 fold change ~ log10 q-values
You can also look at the volcano plot with the q value to control for
false discovery rate (FDR) using the same approach shown for p
values:
fc2 %>%
mutate(
sig = case_when(
qval > 0.05 ~ "not_sig",
qval <= 0.05 & logfc > log2(1.5) ~ "Increased",
qval <= 0.05 & logfc < -log2(1.5) ~ "Decreased",
TRUE ~ "Inconclusive"
),
sig = factor(sig,
levels = c("Increased", "Decreased", "Inconclusive", "not_sig"),
labels = c("Increased", "Decreased", "Inconclusive", "Not significant")
)
) %>%
ggplot() +
aes(x = logfc, y = neg_logq, color = sig) +
geom_point() +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
geom_vline(xintercept = log2(1.5), linetype = "dashed") +
geom_vline(xintercept = -log2(1.5), linetype = "dashed") +
scale_color_manual(values = c("red", "blue", "grey", "black")) +
labs(
x = "log2 Fold Change",
y = "-Log~10~ q-value",
color = "Differential Abundance"
) +
theme_bw() +
theme(
axis.title = element_markdown(size = 20),
axis.text = element_text(size = 15)
)
