In ProMetIS/ProMetIA: Proteomics and Metabolomics Integrated Analysis

knitr::opts_chunk$set(fig.width = 8,
                      fig.height = 8,
                      fig.path = 'figures/temp/5_aggregated_results/',
                      message = FALSE,
                      warning = FALSE)

Loading datasets

latmx2.mset <- phenomis::reading(LATMX2::statistics_intraomics_dir.c())

Building the dataset and .tsv table for each gene and set

for (gene.c in LATMX2::genes.vc()) {

  gene.mset <- ProMetIA::sub_mset(latmx2.mset,
                                  genes.vc = c("WT", gene.c))

  for (set.c in names(gene.mset)) {

    gene.eset <- gene.mset[[set.c]]
    fdata.df <- Biobase::fData(gene.eset)

    gene_discard.vc <- grep(setdiff(LATMX2::genes.vc(), gene.c),
                            colnames(fdata.df), value = TRUE)

    fdata.df <- fdata.df[, setdiff(colnames(fdata.df), gene_discard.vc)]

    if (set.c == "proteomics_liver" && gene.c == "LAT") {

      m_gene_dif.vn <- fdata.df[, paste0("limmaM_WT.", gene.c, "_diff")]
      m_gene_fdr.vn <- fdata.df[, paste0("limmaM_WT.", gene.c, "_BH")]
      m_gene_sig.vi <- fdata.df[, paste0("limmaM_WT.", gene.c, "_signif")]

      f_gene_dif.vn <- fdata.df[, paste0("limmaF_WT.", gene.c, "_diff")]
      f_gene_fdr.vn <- fdata.df[, paste0("limmaF_WT.", gene.c, "_BH")]
      f_gene_sig.vi <- fdata.df[, paste0("limmaF_WT.", gene.c, "_signif")]

      w_sex_dif.vn <- fdata.df[, paste0("limmaWT_M.F_diff_", gene.c)]
      w_sex_fdr.vn <- fdata.df[, paste0("limmaWT_M.F_BH_", gene.c)]
      w_sex_sig.vi <- fdata.df[, paste0("limmaWT_M.F_signif_", gene.c)]

      k_sex_dif.vn <- fdata.df[, paste0("limma", gene.c, "_M.F_diff")]
      k_sex_fdr.vn <- fdata.df[, paste0("limma", gene.c, "_M.F_BH")]
      k_sex_sig.vi <- fdata.df[, paste0("limma", gene.c, "_M.F_signif")]

      gene_fdr_min.vl <- m_gene_fdr.vn < f_gene_fdr.vn

      gene_dif.vn <- as.integer(gene_fdr_min.vl) * m_gene_dif.vn +
        as.integer(!gene_fdr_min.vl) * f_gene_dif.vn
      gene_fdr.vn <- as.integer(gene_fdr_min.vl) * m_gene_fdr.vn +
        as.integer(!gene_fdr_min.vl) * f_gene_fdr.vn
      gene_sig.vi <- as.integer(gene_fdr_min.vl) * m_gene_sig.vi +
        as.integer(!gene_fdr_min.vl) * f_gene_sig.vi      

      sex_fdr_min.vl <- w_sex_fdr.vn < k_sex_fdr.vn

      sex_dif.vn <- as.integer(sex_fdr_min.vl) * w_sex_dif.vn +
        as.integer(!sex_fdr_min.vl) * k_sex_dif.vn
      sex_fdr.vn <- as.integer(sex_fdr_min.vl) * w_sex_fdr.vn +
        as.integer(!sex_fdr_min.vl) * k_sex_fdr.vn
      sex_sig.vi <- as.integer(sex_fdr_min.vl) * w_sex_sig.vi +
        as.integer(!sex_fdr_min.vl) * k_sex_sig.vi

    } else {

      gene_dif.vn <- fdata.df[, paste0("limma2ways_WT.", gene.c, "_diff")]
      gene_fdr.vn <- fdata.df[, paste0("limma2ways_WT.", gene.c, "_BH")]
      gene_sig.vi <- fdata.df[, paste0("limma2ways_WT.", gene.c, "_signif")]

      sex_dif.vn <- fdata.df[, paste0("limma2ways_M.F_diff_", gene.c)]
      sex_fdr.vn <- fdata.df[, paste0("limma2ways_M.F_BH_", gene.c)]
      sex_sig.vi <- fdata.df[, paste0("limma2ways_M.F_signif_", gene.c)]

    }

    hypotest.df <- data.frame(gene_dif = gene_dif.vn,
                              gene_fdr = gene_fdr.vn,
                              gene_sig = gene_sig.vi,
                              sex_dif = sex_dif.vn,
                              sex_fdr = sex_fdr.vn,
                              sex_sig = sex_sig.vi)
    colnames(hypotest.df) <- paste0(rep(c("WT.KO", "M.F"), each = 3),
                                    "_",
                                    rep(c("fold", "BH", "signif"), 2))

    fdata.df <- cbind.data.frame(name = rownames(fdata.df),
                                 set = rep(set.c, nrow(fdata.df)),
                                 gene = rep(gene.c, nrow(fdata.df)),
                                 hypotest.df,
                                 fdata.df)

    fdata.df <- fdata.df[order(gene_fdr.vn,
                               -abs(gene_dif.vn)), ]


    readr::write_tsv(fdata.df,
                     file.path("../../LATMX2/inst/extdata/5_supplementary/tables",
                               paste0(gene.c, "_", set.c, ".tsv")))


    stopifnot(identical(sort(rownames(fdata.df)),
                        sort(Biobase::featureNames(gene.eset))))

    gene.eset <- gene.eset[rownames(fdata.df), ]
    Biobase::fData(gene.eset) <- fdata.df

    gene.mset <- MultiDataSet::add_eset(gene.mset,
                                        gene.eset,
                                        dataset.type = set.c,
                                        GRanges = NA,
                                        overwrite = TRUE,
                                        warnings = FALSE)

  }

  phenomis::writing(gene.mset,
                    file.path("../../LATMX2/inst/extdata/5_supplementary/datasets",
                              gene.c),
                    overwrite.l = TRUE)

}

Graphics

latmx2_mset.ls <- lapply(LATMX2::genes.vc(),
                         function(gene.c)
                           phenomis::reading(file.path(LATMX2::supplementary_datasets_dir.c(), gene.c)))
names(latmx2_mset.ls) <- LATMX2::genes.vc()

feat_maxview.i <- 30

for (gene.c in LATMX2::genes.vc()) {

  message(gene.c)

  gene.mset <- latmx2_mset.ls[[gene.c]]

  for (set.c in names(gene.mset)) {

    message(set.c)

    wtgene.vc <- c("WT", gene.c)

    gene.eset <- gene.mset[[set.c]]

  # discarding the features without statistics for this specific gene

  stat_col.vc <- grep("^(WT.KO_|M.F_)", Biobase::fvarLabels(gene.eset), value = TRUE)

  feat_allna.vl <- apply(Biobase::fData(gene.eset)[, stat_col.vc], 1, function(x) all(is.na(x)))

  gene.eset <- gene.eset[!feat_allna.vl, ]

  # building the data frame for display

  fdr.vn <- Biobase::fData(gene.eset)[, "WT.KO_BH"]
  fc.vn <- Biobase::fData(gene.eset)[, "WT.KO_fold"]
  sex_fdr.vn <- Biobase::fData(gene.eset)[, "M.F_BH"]

  sex_fdr.vc <- sapply(sex_fdr.vn,
                       function(sex_fdr.n) {
                         if (is.na(sex_fdr.n)) {
                           return(", sexNA")
                         } else if (sex_fdr.n <= 0.001) {
                           return(", sex***")
                         } else if (sex_fdr.n <= 0.01) {
                           return(", sex**")
                         } else if (sex_fdr.n <= 0.05) {
                           return(", sex*")
                         } else
                           return("")
                       })

  labels.vc <- paste0(Biobase::featureNames(gene.eset),
                    "\nFDR = ", signif(fdr.vn, 2),
                    ", FC = ", signif(fc.vn, 2),
                    sex_fdr.vc)

  fig_pfx.c <- file.path("../../LATMX2/inst/extdata/5_supplementary/tables",
                         paste0(gene.c, "_", set.c))

  # volcano plot

  for (fig_type.c in c("plotly", "html", "pdf")) {

    figure.c <- switch(fig_type.c,
                       plotly = "interactive_plotly",
                       html = paste0(fig_pfx.c, "_volcano.html"),
                       pdf = paste0(fig_pfx.c, "_volcano.pdf"))

    if (fig_type.c == "pdf") {
      feat_show.vl <- fdr.vn <= 0.05
      fig_labels.vc <- Biobase::featureNames(gene.eset)
      if (sum(feat_show.vl, na.rm = TRUE) <= feat_maxview.i) {
        fig_labels.vc[!feat_show.vl] <- ""
      } else {
        feat_fdr.vn <- fdr.vn
        names(feat_fdr.vn) <- fig_labels.vc
        feat_ord.vn <- feat_fdr.vn[order(feat_fdr.vn)]
        fig_labels.vc[!(fig_labels.vc %in% names(feat_ord.vn)[1:feat_maxview.i])] <- ""
      }
    } else
      fig_labels.vc <- labels.vc

    phenomis::gg_volcanoplot(fold_change.vn = fc.vn,
                             adjusted_pvalue.vn = fdr.vn,
                             adjust_method.c = "BH",
                             adjust_thresh.n = 0.05,
                             label.vc = fig_labels.vc,
                             title.c = set.c,
                             xlab.c = "Fold Change",
                             class_name.vc = c("WT", gene.c),
                             class_color.vc = LATMX2::palette.vc()[c("WT", gene.c)],
                             size.ls = list(class.i = 10,
                                            lab.i = 16,
                                            point.i = 3,
                                            tick.i = 14,
                                            title.i = 20),
                             figure.c = figure.c)
  }

  # individual boxplots

  maxview.i <- max(sum(fdr.vn <= 0.05, na.rm = TRUE), feat_maxview.i)
  view.vi <- 1L:maxview.i
  data.mn <- t(Biobase::exprs(gene.eset))[, view.vi]
  colnames(data.mn) <- paste0("v", view.vi)

  data.df <- data.frame(x = factor(Biobase::pData(gene.eset)[, "gene"],
                                   levels = c("WT", gene.c)),
                        data.mn)

  grDevices::pdf(paste0(fig_pfx.c, "_boxplot.pdf"))

   for (k in view.vi) {

      var.c <- paste0("v", k)

      if (k <= sum(fdr.vn <= 0.05, na.rm = TRUE)) {
        var_palette.vc <- as.character(LATMX2::palette.vc()[c("WT", gene.c)])
      } else
        var_palette.vc <- as.character(LATMX2::palette.vc()[c("WT", "WT")])

      phenomis::gg_boxplot(data.df,
                           x.c = "x",
                           y.c = var.c,
                           color.c = "x",
                           title.c = labels.vc[k],
                           xlab.c = NA,
                           ylab.c = "",
                           palette.vc = var_palette.vc,
                           size.ls = list(dot.n = 0.7,
                                          lab.i = 20,
                                          tick.i = 20,
                                          title.i = 20),
                           figure.c = "interactive")

   }

  grDevices::dev.off()

  }

}