MetaboDiff case study
8 July 2019
Abstract
This markdown documents presents the code to generate the plots presented in the supplementary data of the MetaboDiff publication.
Package
MetaboDiff 0.9.3
library("MetaboDiff")
source("http://peterhaschke.com/Code/multiplot.R")1 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
1.1 Objects
case1_cells## A MultiAssayExperiment object of 1 listed
## experiment with a user-defined name and respective class.
## Containing an ExperimentList class object of length 1:
## [1] raw: SummarizedExperiment with 133 rows and 15 columns
## Features:
## experiments() - obtain the ExperimentList instance
## colData() - the primary/phenotype DataFrame
## sampleMap() - the sample availability DataFrame
## `$`, `[`, `[[` - extract colData columns, subset, or experiment
## *Format() - convert into a long or wide DataFrame
## assays() - convert ExperimentList to a SimpleList of matricestable(case1_cells$group)##
## AKT1-high Control MYC-high
## 5 5 5case1_mice## A MultiAssayExperiment object of 1 listed
## experiment with a user-defined name and respective class.
## Containing an ExperimentList class object of length 1:
## [1] raw: SummarizedExperiment with 170 rows and 18 columns
## Features:
## experiments() - obtain the ExperimentList instance
## colData() - the primary/phenotype DataFrame
## sampleMap() - the sample availability DataFrame
## `$`, `[`, `[[` - extract colData columns, subset, or experiment
## *Format() - convert into a long or wide DataFrame
## assays() - convert ExperimentList to a SimpleList of matricestable(case1_mice$group)##
## AKT1-high Control MYC-high
## 6 6 6case1_human <- create_mae(assay,rowData,colData)## Warning in MultiAssayExperiment(experiments = experiment_list, colData =
## colData, : sampleMap[['primary']] coerced to character()## Warning in MultiAssayExperiment(experiments = experiment_list, colData =
## colData, : sampleMap[['colname']] coerced to character()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)##
## AKT1-high Control MYC-high
## 21 25 91.2 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)1.3 Create list
case1 = list(case1_cells,case1_mice,case1_human)
names(case1) = c("cells","mice","human")1.4 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)1.5 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"))## $cells
##
## $mice
##
## $human
#dev.off()case1 <- sapply(case1,knn_impute,cutoff=0.4)1.6 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)## $cells
##
## $mice
##
## $human
#dev.off()1.7 Normalization
case1 <- sapply(case1, normalize_met)## Warning in vsnSample(v): 6 rows were removed since they contained only NA
## elements.## Warning in vsnSample(v): 4 rows were removed since they contained only NA
## elements.1.8 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()1.9 Hypothesis testing
case1 <- sapply(case1,
diff_test,
group_factors=c("group","tumor_groups"))1.10 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()1.11 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]]## [1] 0.03871038plot(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]]## [1] 0.01298886text(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]]## [1] 0.00636392text(x=2,y=22.0,"**",cex=2)
#dev.off()1.12 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"))## ## alpha: 1.000000
## ..cutHeight not given, setting it to 0.982 ===> 99% of the (truncated) height range in dendro.
## ..done.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"))## alpha: 1.000000
## ..cutHeight not given, setting it to 0.984 ===> 99% of the (truncated) height range in dendro.
## ..done.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"))## alpha: 1.000000
## ..cutHeight not given, setting it to 0.989 ===> 99% of the (truncated) height range in dendro.
## ..done.1.13 Venn diagram
library(VennDiagram)## Loading required package: futile.loggerA = 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()1.14 Properties of correlation networks
table(metadata(case1$cells)$modules)##
## 0 1 2 3 4 5 6 7 8 9 10 11
## 1 22 19 13 10 9 9 8 8 7 7 51.15 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
)
## cells mice human
## 0 0 0#dev.off()1.16 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)## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 2 rows containing missing values (geom_point).## Warning: Removed 2 rows containing missing values (geom_text).## Warning: Removed 7 rows containing missing values (geom_point).## Warning: Removed 7 rows containing missing values (geom_text).## Warning: Removed 7 rows containing missing values (geom_point).## Warning: Removed 7 rows containing missing values (geom_text).
#dev.off()Session information
sessionInfo()## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] VennDiagram_1.6.20 futile.logger_1.4.3
## [3] MetaboDiff_0.9.3 forcats_0.3.0
## [5] stringr_1.4.0 dplyr_0.8.0.1
## [7] purrr_0.3.2 readr_1.3.1
## [9] tidyr_0.8.2 tibble_2.1.1
## [11] ggplot2_3.1.0 tidyverse_1.2.1
## [13] ComplexHeatmap_1.99.8 MultiAssayExperiment_1.8.1
## [15] SummarizedExperiment_1.12.0 DelayedArray_0.8.0
## [17] BiocParallel_1.16.5 matrixStats_0.54.0
## [19] Biobase_2.42.0 GenomicRanges_1.34.0
## [21] GenomeInfoDb_1.18.1 IRanges_2.16.0
## [23] S4Vectors_0.20.1 BiocGenerics_0.28.0
## [25] BiocStyle_2.10.0
##
## loaded via a namespace (and not attached):
## [1] readxl_1.2.0 backports_1.1.3 circlize_0.4.6
## [4] fastmatch_1.1-0 Hmisc_4.2-0 igraph_1.2.2
## [7] plyr_1.8.4 lazyeval_0.2.2 splines_3.5.2
## [10] robust_0.4-18 digest_0.6.19 foreach_1.4.4
## [13] BiocInstaller_1.30.0 htmltools_0.3.6 GO.db_3.6.0
## [16] phytools_0.6-60 magrittr_1.5 checkmate_1.9.1
## [19] memoise_1.1.0 fit.models_0.5-14 cluster_2.0.7-1
## [22] doParallel_1.0.14 limma_3.36.5 fastcluster_1.1.25
## [25] annotate_1.58.0 modelr_0.1.2 colorspace_1.4-1
## [28] blob_1.1.1 rvest_0.3.2 rrcov_1.4-7
## [31] haven_2.0.0 xfun_0.7 crayon_1.3.4
## [34] RCurl_1.95-4.11 jsonlite_1.6 genefilter_1.62.0
## [37] impute_1.54.0 phangorn_2.4.0 ape_5.2
## [40] survival_2.43-3 iterators_1.0.10 glue_1.3.1
## [43] gtable_0.3.0 zlibbioc_1.28.0 XVector_0.22.0
## [46] GetoptLong_0.1.7 shape_1.4.4 maps_3.3.0
## [49] DEoptimR_1.0-8 scales_1.0.0 futile.options_1.0.1
## [52] vsn_3.48.1 mvtnorm_1.0-10 DBI_1.0.0
## [55] Rcpp_1.0.1 plotrix_3.7-4 xtable_1.8-3
## [58] htmlTable_1.13.1 clue_0.3-57 foreign_0.8-71
## [61] bit_1.1-14 preprocessCore_1.42.0 Formula_1.2-3
## [64] animation_2.6 htmlwidgets_1.3 httr_1.4.0
## [67] RColorBrewer_1.1-2 acepack_1.4.1 pkgconfig_2.0.2
## [70] XML_3.98-1.16 nnet_7.3-12 dynamicTreeCut_1.63-1
## [73] tidyselect_0.2.5 labeling_0.3 rlang_0.3.4
## [76] AnnotationDbi_1.44.0 munsell_0.5.0 cellranger_1.1.0
## [79] tools_3.5.2 cli_1.1.0 generics_0.0.2
## [82] RSQLite_2.1.1 broom_0.5.1 evaluate_0.14
## [85] yaml_2.2.0 knitr_1.23 bit64_0.9-7
## [88] robustbase_0.93-4 nlme_3.1-137 formatR_1.5
## [91] xml2_1.2.0 compiler_3.5.2 rstudioapi_0.10
## [94] png_0.1-7 affyio_1.50.0 clusterGeneration_1.3.4
## [97] pcaPP_1.9-73 stringi_1.4.3 lattice_0.20-38
## [100] Matrix_1.2-15 pillar_1.3.1 BiocManager_1.30.4
## [103] combinat_0.0-8 GlobalOptions_0.1.0 data.table_1.12.0
## [106] bitops_1.0-6 R6_2.4.0 latticeExtra_0.6-28
## [109] affy_1.58.0 bookdown_0.11.1 gridExtra_2.3
## [112] codetools_0.2-15 lambda.r_1.2.3 MASS_7.3-51.1
## [115] assertthat_0.2.1 rjson_0.2.20 withr_2.1.2
## [118] mnormt_1.5-5 GenomeInfoDbData_1.2.0 expm_0.999-3
## [121] hms_0.4.2 quadprog_1.5-5 rpart_4.1-13
## [124] coda_0.19-2 rmarkdown_1.13 scatterplot3d_0.3-41
## [127] numDeriv_2016.8-1 WGCNA_1.68 lubridate_1.7.4
## [130] base64enc_0.1-3