In andreasmock/MetaboDiff: An R package for differential metabolomic analysis

library("MetaboDiff")
source("http://peterhaschke.com/Code/multiplot.R")

Case study

Priolo, C., Pyne, S., Rose, J., Regan, E. R., Zadra, G., Photopoulos, C., et al. (2014). AKT1 and MYC Induce Distinctive Metabolic Fingerprints in Human Prostate Cancer. Cancer Research, 74(24), 7198–7204. http://doi.org/10.1158/0008-5472.CAN-14-1490

Research question: AKT1 and MYC are two of the most common prostate cancer oncogenes. Priolo and colleagues investigated the metabolic profiles of AKT1- and MYC-driven

Objects

case1_cells
table(case1_cells$group)

case1_mice
table(case1_mice$group)

case1_human <- create_mae(assay,rowData,colData)
case1_human <- case1_human[,colData(case1_human)$group %in% c("Control","MYC-high","AKT1-high")]
case1_human$group <- droplevels(case1_human$group)
table(case1_human$group)

Metabolite annotation

case1_cells <- get_SMPDBanno(case1_cells,3,4,NA)
case1_mice <- get_SMPDBanno(case1_mice,3,4,NA)
case1_human <- get_SMPDBanno(case1_human,6,7,NA)

Create list

case1 = list(case1_cells,case1_mice,case1_human)
names(case1) = c("cells","mice","human")

Add tumor grouping

case1$cells$tumor_groups = case1$cells$group
case1$cells$tumor_groups[case1$cells$tumor_groups=="Control"] = NA
case1$cells$tumor_groups = droplevels(case1$cells$tumor_groups)

case1$mice$tumor_groups = case1$mice$group
case1$mice$tumor_groups[case1$mice$tumor_groups=="Control"] = NA
case1$mice$tumor_groups = droplevels(case1$mice$tumor_groups)

case1$human$tumor_groups = case1$human$group
case1$human$tumor_groups[case1$human$tumor_groups=="Control"] = NA
case1$human$tumor_groups = droplevels(case1$human$tumor_groups)

Imputation of missing values

group_factor = "group"
label_colors = c("orange","dodgerblue","darkseagreen")

#pdf("../../MetaboDiff_paper/re_submission/imputation.pdf",height=4)
sapply(case1,
       na_heatmap,
       group_factor="group",
       label_colors = c("orange","dodgerblue","darkseagreen"))
#dev.off()

case1 <- sapply(case1,knn_impute,cutoff=0.4)

Heatmap

#pdf("../../MetaboDiff_paper/re_submission/hms.pdf",height=4)
sapply(case1,
       outlier_heatmap,
       group_factor="group",
       label_colors = c("orange","dodgerblue","darkseagreen"),
       k=3)
#dev.off()

Normalization

case1 <- sapply(case1, normalize_met)

PCA

#pdf("../../MetaboDiff_paper/re_submission/pca.pdf",width=7,height=5)
multiplot(
pca_plot(case1$cells,group_factor="group",
       label_colors = c("orange","dodgerblue","darkseagreen")) + ggtitle("cells"),
pca_plot(case1$human,group_factor="group",
       label_colors = c("orange","dodgerblue","darkseagreen"))+ ggtitle("human"),
pca_plot(case1$mice,group_factor="group",
       label_colors = c("orange","dodgerblue","darkseagreen"))+ ggtitle("mice"),
cols=2)
#dev.off()

Hypothesis testing

case1 <- sapply(case1,
       diff_test, 
       group_factors=c("group","tumor_groups"))

Volcano plot

#pdf("../../MetaboDiff_paper/re_submission/vp.pdf",width=6,height=7)
par(mfrow=c(3,2))
volcano_plot(case1$cells,
             group_factor="tumor_groups",
             label_colors = c("darkseagreen","orange"),
             main="cells")

volcano_plot(case1$cells,
             group_factor="tumor_groups",
              label_colors = c("darkseagreen","orange"),
             main="cells",
             p_adjust = FALSE)

volcano_plot(case1$mice,
             group_factor="tumor_groups",
             label_colors = c("darkseagreen","orange"),
             main="mice")

volcano_plot(case1$mice,
             group_factor="tumor_groups",
             label_colors = c("darkseagreen","orange"),
             main="mice",
             p_adjust = FALSE)

volcano_plot(case1$human,
             group_factor="tumor_groups",
             label_colors = c("darkseagreen","orange"),
             main="human")
volcano_plot(case1$human,
             group_factor="tumor_groups",
              label_colors = c("darkseagreen","orange"),
             main="human",
             p_adjust = FALSE)
#dev.off()

Figure 3B

#pdf("../../MetaboDiff_paper/re_submission/3b.pdf",width=7,height=3)
ids = case1$human$group %in% c("AKT1-high","MYC-high")
par(mfrow=c(1,3))
plot(assay(case1$human[["norm_imputed"]])[61,ids]~droplevels(case1$human$group[ids]),
     main="Arachidonic acid",xlab="",ylab="Normalized values",frame=FALSE,col=c("orange","darkseagreen"))
text(x=2,y=25.4,"*",cex=2)
t.test(assay(case1$human[["norm_imputed"]])[61,ids]~droplevels(case1$human$group[ids]))[[3]]
plot(assay(case1$human[["norm_imputed"]])[93,ids]~droplevels(case1$human$group[ids]),
     main="Docohexaenoic acid",xlab="",ylab="Normalized values",frame=FALSE,col=c("orange","darkseagreen"))
t.test(assay(case1$human[["norm_imputed"]])[93,ids]~droplevels(case1$human$group[ids]))[[3]]
text(x=2,y=23.7,"*",cex=2)
plot(assay(case1$human[["norm_imputed"]])[179,ids]~droplevels(case1$human$group[ids]),
     main="Oleic acid",xlab="",ylab="Normalized values",frame=FALSE,col=c("orange","darkseagreen"))
t.test(assay(case1$human[["norm_imputed"]])[179,ids]~droplevels(case1$human$group[ids]))[[3]]
text(x=2,y=22.0,"**",cex=2)
#dev.off()

Metabolic correlation networks

case1$cells <- case1$cells %>%
    diss_matrix %>%
    identify_modules(min_module_size=5) %>%
    name_modules(pathway_annotation="SUB_PATHWAY") %>%
    calculate_MS(group_factors=c("group", "tumor_groups"))

case1$mice <- case1$mice %>%
    diss_matrix %>%
    identify_modules(min_module_size=5) %>%
    name_modules(pathway_annotation="SUB_PATHWAY") %>%
    calculate_MS(group_factors=c("group", "tumor_groups"))

case1$human <- case1$human %>%
    diss_matrix %>%
    identify_modules(min_module_size=5) %>%
    name_modules(pathway_annotation="SUB_PATHWAY") %>%
    calculate_MS(group_factors=c("group", "tumor_groups"))

Venn diagram

library(VennDiagram)
A = as.character(rowData(case1$cells[["norm_imputed"]])$BIOCHEMICAL[metadata(case1$cells)$`ttest_tumor_groups_MYC-high_vs_AKT1-high`$pval<0.05])
B = as.character(rowData(case1$mice[["norm_imputed"]])$BIOCHEMICAL[metadata(case1$mice)$`ttest_tumor_groups_MYC-high_vs_AKT1-high`$pval<0.05])
C = as.character(rowData(case1$human[["norm_imputed"]])$BIOCHEMICAL[metadata(case1$human)$`ttest_tumor_groups_MYC-high_vs_AKT1-high`$pval<0.05])
venn <- draw.triple.venn(area1 = length(A),
                         area2 = length(B),
                         area3 = length(C),
                         n12 = length(intersect(A,B)),
                         n13 = length(intersect(A,C)),
                         n23 = length(intersect(B,C)),
                         n123 = length(intersect(intersect(A,B),C)),
                         fill = c("dodgerblue","red3","yellow"),
                         alpha=c(0.1,0.1,0.1),
                         category = c("cells","mice","human"),
                         lwd = c(0.5,0.5,0.5),
                         cex = 1.3,
                         fontfamily = "sans",
                         cat.cex = 1.3
                         )
#pdf("../../MetaboDiff_paper/re_submission/venn.pdf",width=4,height=4)
grid.draw(venn)
grid.newpage()
#dev.off()

Properties of correlation networks

table(metadata(case1$cells)$modules)

MS plot

#pdf("../../MetaboDiff_paper/re_submission/ms.pdf",width=8,height=4)
sapply(case1,
       MS_plot,
       group_factor="tumor_groups",
       p_value_cutoff=0.1,
       p_adjust=FALSE
)
#dev.off()

MOI plots

#pdf("../../MetaboDiff_paper/re_submission/MOI.pdf",width=8,height=13)
multiplot(
    MOI_plot(case1$cells,
         group_factor="group",
         MOI = 2,
         label_colors=c("darkseagreen","orange"),
         p_adjust = TRUE) + xlim(c(-1,7.5)) + ggtitle("cells") +
    ylim(c(0.5,1)),
MOI_plot(case1$mice,
         group_factor="group",
         MOI = 1,
         label_colors=c("darkseagreen","orange"),
         p_adjust = TRUE) + xlim(c(-1,7.5)) + ggtitle("mice") +
    ylim(c(0.2,0.5)),
MOI_plot(case1$human,
         group_factor="group",
         MOI = 3,
         label_colors=c("darkseagreen","orange"),
         p_adjust = TRUE) + xlim(c(-1,7.5)) + ggtitle("human")+
    ylim(c(0.6,0.9)),
cols=1)
#dev.off()

Session information {.unnumbered}

sessionInfo()

[^1]: Zhang, B., & Horvath, S. (2005). A general framework for weighted gene co-expression network analysis. Statistical Applications in Genetics and Molecular Biology, 4(1), Article17. http://doi.org/10.2202/1544-6115.1128

[^2]: Langfelder, P., & Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics, 9, 559–559. http://doi.org/10.1186/1471-2105-9-559

[^3]: Priolo, C., Pyne, S., Rose, J., Regan, E. R., Zadra, G., Photopoulos, C., et al. (2014). AKT1 and MYC Induce Distinctive Metabolic Fingerprints in Human Prostate Cancer. Cancer Research, 74(24), 7198–7204. http://doi.org/10.1158/0008-5472.CAN-14-1490

[^4]: Zheng, C.-H., Yuan, L., Sha, W., & Sun, Z.-L. (2014). Gene differential coexpression analysis based on biweight correlation and maximum clique. BMC Bioinformatics, 15 Suppl 15(Suppl 15), S3. http://doi.org/10.1186/1471-2105-15-S15-S3

[^5]: Langfelder, P., Zhang, B., & Horvath, S. (2008). Defining clusters from a hierarchical cluster tree: the Dynamic Tree Cut package for R. Bioinformatics, 24(5), 719–720. http://doi.org/10.1093/bioinformatics/btm563.

[^6]: Horvath, S., & Dong, J. (2008). Geometric Interpretation of Gene Coexpression Network Analysis. PLoS Computational Biology (PLOSCB) 4(8), 4(8), e1000117–e1000117. http://doi.org/10.1371/journal.pcbi.1000117

andreasmock/MetaboDiff documentation built on Oct. 31, 2020, 8:18 a.m.