## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, fig.width = 6, fig.height = 6) ## ----install, eval=FALSE------------------------------------------------------ # if (!require("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # # BiocManager::install("mgcca") ## ----load_library------------------------------------------------------------- library(mgcca) ## ----cardio_data-------------------------------------------------------------- data(cardiovascular, package = "mgcca") X <- list(methylation = as.matrix(X1), clinical = as.matrix(X2), other = as.matrix(X3)) X <- lapply(X, function(m) { storage.mode(m) <- "double"; m }) # Different number of individuals per table -> missing individuals sapply(X, dim) length(Reduce(union, lapply(X, rownames))) # size of the individual union ## ----cardio_run--------------------------------------------------------------- h5file <- tempfile(fileext = ".h5") fit <- mgcca(X, filename = h5file, nfac = 2, method = "solve", scores = TRUE, collect = TRUE) fit # print method: a compact overview ## ----cardio_summary----------------------------------------------------------- summary(fit, top = 4) ## ----cardio_clin_load--------------------------------------------------------- round(fit$corsY[["clinical"]], 3) ## ----cardio_group------------------------------------------------------------- ldl <- X$clinical[, "LDL"] # real LDL values (mg/dL) ldl <- ldl[rownames(fit$Y)] # align to the individual union ldl_group <- cut(ldl, breaks = c(-Inf, 100, 160, Inf), labels = c("optimal", "borderline", "high")) table(ldl_group, useNA = "ifany") ## ----cardio_plotInds, fig.width=6, fig.height=4.6----------------------------- plotIndividuals(fit, group = ldl_group) ## ----cardio_plotInds_bmi, fig.width=6, fig.height=4.6------------------------- bmi <- X$clinical[, "BMI"][rownames(fit$Y)] bmi_group <- cut(bmi, breaks = c(-Inf, 25, 30, Inf), labels = c("normal", "overweight", "obese")) plotIndividuals(fit, group = bmi_group, title = "Individuals by BMI category") ## ----cardio_plotVars, fig.width=5.5, fig.height=5.5--------------------------- plotVariables(fit, table = "clinical", top = 6) ## ----cardio_loadings, fig.width=6, fig.height=4------------------------------- plotLoadings(fit, table = "clinical", comp = 1, top = 6) ## ----cardio_ave, fig.width=6, fig.height=4------------------------------------ plotAVE(fit) plotScree(fit) ## ----cardio_signif------------------------------------------------------------ head(getSignif(fit, pval.cut = 1e-3)) topVars(fit, axis = 1, topN = 5) ## ----cardio_cortestp---------------------------------------------------------- r <- fit$corsY[["clinical"]][, 1] p <- corTestP(r, nrow(X$clinical)) # n of the table the correlation came from all.equal(p, fit$pval.cor[["clinical"]][, 1]) signif(head(p, 4), 3) ## ----cardio_restyle, fig.width=6, fig.height=5, fig.retina=1------------------ g <- plotIndividuals(fit, group = ldl_group, ellipse = FALSE) g + ggplot2::labs(subtitle = "LDL category (NCEP ATP III cut-points)") + ggplot2::theme(legend.position = "bottom") ## ----reload------------------------------------------------------------------- # forget the in-memory object; everything needed is on disk rm(fit) fit <- mgccaLoad(h5file) # rebuilt straight from the HDF5 file fit ## ----reload_desc-------------------------------------------------------------- desc <- attr(fit, "desc") desc$method # "solve" desc$route # "auto" desc$route_dual # per-table: covariance vs Gram route desc$eig_values # shared-component eigenvalues ## ----reload_use, fig.width=6, fig.height=4.6---------------------------------- plotIndividuals(fit, group = ldl_group) ## ----results_reader----------------------------------------------------------- res <- mgccaResults(h5file, datasets = names(X), outputs = "Y") dim(res$Y) identical(res$Y, fit$Y) # same numbers, read back from the file ## ----acc_load----------------------------------------------------------------- library(MultiAssayExperiment) library(SummarizedExperiment) data("miniACC", package = "MultiAssayExperiment") miniACC ## ----acc_subset--------------------------------------------------------------- acc_mae <- miniACC[, , c("RNASeq2GeneNorm", "RPPAArray")] ## ----acc_impute--------------------------------------------------------------- acc_mae <- mgccaImpute(acc_mae, method = "knn", k = 10, impute.zero = TRUE) ## ----acc_tables--------------------------------------------------------------- tabs <- getTables(acc_mae) class(tabs) vapply(tabs, dim, integer(2)) # individuals x variables for (a in seq_along(tabs)) rownames(tabs[[a]]) <- substr(rownames(tabs[[a]]), 1, 12) length(Reduce(union, lapply(tabs, rownames))) # individuals in the analysis length(Reduce(intersect, lapply(tabs, rownames))) # measured on both layers ## ----acc_chr2num-------------------------------------------------------------- storage.mode(tabs[[1]]) <- "character" tabs[[1]][1, 1] tabs <- matrixChr2Num(tabs, index = 1) tabs[[1]][1, 1] ## ----acc_lambda, message=TRUE------------------------------------------------- sel <- mgccaSelectLambda(tabs, nfac = 2, K = 5, seed = 1) sel ## ----acc_fit------------------------------------------------------------------ h5acc <- tempfile(fileext = ".h5") acc <- mgcca(tabs, filename = h5acc, nfac = 2, method = "penalized", lambda = sel$lambda, scores = TRUE, route = "auto", collect = TRUE) acc ## ----acc_summary-------------------------------------------------------------- summary(acc, top = 4) ## ----acc_results-------------------------------------------------------------- res <- mgccaResults(h5acc, datasets = names(tabs), outputs = c("Y", "AVE")) dim(res$Y) names(Filter(Negate(is.null), res)) ## ----acc_pheno---------------------------------------------------------------- ids <- rownames(acc$Y) cd <- as.data.frame(colData(miniACC)) stage <- factor(as.character(cd[ids, "pathologic_stage"])) vital <- factor(ifelse(cd[ids, "vital_status"] == 1, "deceased", "alive")) gender <- factor(as.character(cd[ids, "gender"])) table(stage, useNA = "ifany") ## ----acc_scores, fig.width=6, fig.height=4.6---------------------------------- plotScores(acc, table = 1, group = vital) ## ----acc_biplot, fig.width=6, fig.height=5.2---------------------------------- plotBiplot(acc, table = 2, top = 6) ## ----acc_assoc---------------------------------------------------------------- assoc <- mgccaAssociate(acc, list(stage = stage, vital = vital, gender = gender)) assoc ## ----acc_permtest------------------------------------------------------------- set.seed(1) pt <- mgccaPermtest(tabs, nperm = 9, nfac = 2, method = "penalized", lambda = sel$lambda) pt$p.value range(pt$null) ## ----acc_predict-------------------------------------------------------------- new_rows <- tabs$RPPAArray[1:3, , drop = FALSE] predict(acc, newdata = list(RPPAArray = new_rows)) ## ----acc_cleanup-------------------------------------------------------------- BigDataStatMeth::hdf5_close_all() unlink(h5acc) ## ----audit_design------------------------------------------------------------- ids <- sprintf("ind%02d", 1:85) core <- 1:60 # methylation + expression + proteins extra <- 61:85 # proteins + metabolites avail <- cbind( methylation = seq_along(ids) %in% core, expression = seq_along(ids) %in% core, proteins = rep(TRUE, length(ids)), metabolites = seq_along(ids) %in% extra) rownames(avail) <- ids cohort <- rep(c("cohort A", "cohort B"), length.out = length(ids)) avail_w <- ifelse(seq_along(ids) %in% core, 1, 3) # declared design weights colSums(avail) ## ----audit_run---------------------------------------------------------------- audit <- mgccaAudit(avail, L = 2, strata = cohort, weights = avail_w) audit ## ----audit_paths-------------------------------------------------------------- audit$provenance$pairs[, c("block1", "block2", "code", "n_jk", "path", "bottleneck")] ## ----audit_fragility---------------------------------------------------------- audit$graph$bridges audit$graph$articulation ## ----audit_summary------------------------------------------------------------ summary(audit) ## ----audit_plot, fig.wide=TRUE, fig.width=8, fig.height=8.2, out.width="100%"---- plot(audit) ## ----audit_matrix, fig.wide=TRUE, fig.width=8, fig.height=6.6, out.width="100%"---- plot(audit, type = "matrix") ## ----audit_graph, fig.wide=TRUE, fig.width=8, fig.height=6.2, out.width="100%"---- plot(audit, type = "graph") ## ----audit_provenance, fig.wide=TRUE, fig.width=8, fig.height=6.2, out.width="100%"---- plot(audit, type = "graph", pair = c("methylation", "metabolites")) ## ----audit_availability, fig.wide=TRUE, fig.width=8, fig.height=5.0, out.width="100%"---- plot(audit, type = "availability") ## ----audit_representation, fig.wide=TRUE, fig.width=8, fig.height=5.6, out.width="100%"---- plot(audit, type = "representation") ## ----audit_patterns, fig.wide=TRUE, fig.width=8, fig.height=4.2, out.width="100%"---- plot(audit, type = "patterns") ## ----audit_roundtrip---------------------------------------------------------- set.seed(1) blocks <- list( methylation = matrix(rnorm(60 * 4), 60, 4, dimnames = list(ids[core], paste0("cg", 1:4))), expression = matrix(rnorm(60 * 3), 60, 3, dimnames = list(ids[core], paste0("gene", 1:3))), proteins = matrix(rnorm(85 * 3), 85, 3, dimnames = list(ids, paste0("prot", 1:3))), metabolites = matrix(rnorm(25 * 2), 25, 2, dimnames = list(ids[extra], paste0("met", 1:2)))) h5_audit <- tempfile(fileext = ".h5") desc <- mgccaImportHDF5(blocks, filename = h5_audit, overwriteFile = TRUE) stored <- mgccaAudit(desc, L = 2, strata = cohort, weights = avail_w, save_hdf5 = TRUE) # the availability table derived from the blocks is the one built by hand identical(unname(stored$table), unname(avail)) # and the stored audit comes back whole identical(mgccaAuditLoad(h5_audit), stored) ## ----audit_stored------------------------------------------------------------- stored$saved$tables rownames(as.matrix( BigDataStatMeth::hdf5_matrix(h5_audit, "availability_audit/pair_codes"))) BigDataStatMeth::hdf5_close_all() unlink(h5_audit) ## ----tcga_prep, eval=FALSE---------------------------------------------------- # library(curatedTCGAData); library(TCGAutils); library(MultiAssayExperiment) # # mae <- curatedTCGAData("ACC", c("RNASeq2GeneNorm", "Methylation"), # dry.run = FALSE, version = "1.1.38") # mae <- mgccaImpute(mae, method = "knn", k = 10, rowmax = 0.5, colmax = 0.8, # remove.col = TRUE, impute.zero = TRUE) # # # individuals x variables, aligned by the 12-char sample ID; methylation -> numeric # tables_list <- getTables(mae) # for (a in seq_along(tables_list)) # rownames(tables_list[[a]]) <- substr(rownames(tables_list[[a]]), 1, 12) # tables_list <- matrixChr2Num(tables_list, index = 1) ## ----tcga_run, eval=FALSE----------------------------------------------------- # sel <- mgccaSelectLambda(tables_list, nfac = 2, K = 5, seed = 1) # # fit <- mgcca(tables_list, filename = tempfile(fileext = ".h5"), nfac = 2, # method = "penalized", lambda = sel$lambda, # scores = TRUE, route = "auto", collect = TRUE) # # attr(fit, "desc")$route_dual # both TRUE: methylation and RNA-seq both p >> n # summary(fit, top = 4) ## ----tcga_run_plots, eval=FALSE----------------------------------------------- # stage <- factor(mae$pathologic_stage) # vital <- factor(ifelse(mae$vital_status == "1", "deceased", "alive")) # plotIndividuals(fit, group = stage) # figure 1 # plotIndividuals(fit, group = vital) # figure 2 ## ----tcga_fig_inds, echo=FALSE, out.width="49%", fig.show="hold"-------------- knitr::include_graphics(c("figures/tcga_stage.png", "figures/tcga_individuals.png")) ## ----tcga_assoc, eval=FALSE--------------------------------------------------- # mgccaAssociate(fit, list( # stage = factor(mae$pathologic_stage), # vital = factor(mae$vital_status), # gender = factor(mae$gender))) ## ----tcga_run_vars, eval=FALSE------------------------------------------------ # plotVariables(fit, table = "ACC_Methylation-20160128", top = 15) # plotLoadings(fit, table = "ACC_Methylation-20160128", comp = 1, top = 15) # plotAVE(fit) ## ----tcga_fig_vars, echo=FALSE, out.width="49%", fig.show="hold"-------------- knitr::include_graphics(c("figures/tcga_variables.png", "figures/tcga_loadings.png")) ## ----tcga_fig_ave, echo=FALSE, out.width="70%", fig.align="center"------------ knitr::include_graphics("figures/tcga_ave.png") ## ----tcga_ranges, eval=FALSE-------------------------------------------------- # library(GenomicRanges) # library(TxDb.Hsapiens.UCSC.hg19.knownGene) # library(org.Hs.eg.db) # # # MEN1 +- 2kb # chr.coords <- GRanges("11:64550986-64578766") # # # Annotate RNA-seq with ranges; methylation already carries CpG coordinates # mae.annotated <- TCGAutils::symbolsToRanges(mae) # mae.annotated[[3]] <- NULL # drop the unranged part # # mae.subset <- genCoords(multiassayexperiment = mae.annotated, # gen.coords = chr.coords, # meth.index = 1, col.name = "Genomic_Coordinate") # # mae.subset.imputed <- mgccaImpute(mae.subset, method = "knn", k = 10, # rowmax = 0.5, colmax = 0.8, # remove.col = TRUE, impute.zero = TRUE) # # sub_tables <- getTables(mae.subset.imputed) # for (a in seq_along(sub_tables)) # rownames(sub_tables[[a]]) <- substr(rownames(sub_tables[[a]]), 1, 12) # sub_tables <- matrixChr2Num(sub_tables, index = 1) # # sub_sel <- mgccaSelectLambda(sub_tables, nfac = 2, K = 5, seed = 1) # fit.region <- mgcca(sub_tables, filename = tempfile(fileext = ".h5"), # nfac = 2, method = "penalized", lambda = sub_sel$lambda, # scores = TRUE, collect = TRUE) # plotIndividuals(fit.region, group = factor(mae$vital_status)) ## ----sesinfo------------------------------------------------------------------ sessionInfo()