mqppep_anova: mqppep_anova_script.Rmd comparison

comparison mqppep_anova_script.Rmd @ 22:61adb8801b73 draft

planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/mqppep commit 0c7ca054e77e042c8a584c9903073da064df7d8b

author	eschen42
date	Thu, 30 Jun 2022 16:16:32 +0000
parents	2c5f1a2fe16a
children	8582a9797c18

comparison

equal deleted inserted replaced

-:f3deca1a3c84
+:61adb8801b73
 ---
-title: "MaxQuant Phospho-Proteomic Enrichment Pipeline ANOVA"
+title: "MaxQuant Phosphoproteomic Enrichment Pipeline ANOVA/KSEA"
-author: "Larry Cheng; Art Eschenlauer"
+author:
-date: "May 28, 2018; Mar 16, 2022"
+- "Nick Graham^[ORCiD 0000-0002-6811-1941, University of Southern California: Los Angeles, CA, US]"
+- "Larry Cheng^[ORCiD 0000-0002-6922-6433, Rutgers School of Graduate Studies: New Brunswick, NJ, US]"
+- "Art Eschenlauer^[ORCiD 0000-0002-2882-0508, University of Minnesota: Minneapolis, Minnesota, US]"
+date:
+- "May 28, 2018"
+- "; revised June 23, 2022"
 output:
 pdf_document:
 toc: true
-latex_document:
+toc_depth: 3
-toc: true
+keep_tex: true
+header-includes:
+- \usepackage{longtable}
+- \newcommand\T{\rule{0pt}{2.6ex}}       % Top strut
+- \newcommand\B{\rule[-1.2ex]{0pt}{0pt}} % Bottom strut
 params:
-alphaFile: "test-data/alpha_levels.tabular"
+alphaFile:            "test-data/alpha_levels.tabular"
-inputFile: "test-data/UT_Phospho_ST_Sites.preproc.tabular"
+inputFile:            "test-data/test_input_for_anova.tabular"
-firstDataColumn: "^Intensity[^_]"
+preprocDb:            "test-data/test_input_for_anova.sqlite"
-imputationMethod: !r c("group-median", "median", "mean", "random")[4]
+kseaAppPrepDb:        !r c(":memory:", "test-data/mqppep.sqlite")[2]
-meanPercentile: 1
+show_toc:             true
-sdPercentile: 1.0
+firstDataColumn:      "^Intensity[^_]"
-regexSampleNames: "\\.\\d+[A-Z]$"
+imputationMethod:     !r c("group-median", "median", "mean", "random")[1]
-regexSampleGrouping: "\\d+"
+meanPercentile:       1
-imputedDataFilename: "test-data/limbo/imputedDataFilename.txt"
+sdPercentile:         1.0
-imputedQNLTDataFile: "test-data/limbo/imputedQNLTDataFile.txt"
+regexSampleNames:     "\\.\\d+[A-Z]$"
-show_toc: true
+regexSampleGrouping:  "\\d+"
+imputedDataFilename:  "test-data/limbo/imputedDataFilename.txt"
+imputedQNLTDataFile:  "test-data/limbo/imputedQNLTDataFile.txt"
+anovaKseaMetadata:    "test-data/limbo/anovaKseaMetadata.txt"
+oneWayManyCategories: !r c("aov", "kruskal.test", "oneway.test")[1]
+oneWayTwoCategories:  !r c("aov", "kruskal.test", "oneway.test")[3]
+kseaCutoffStatistic:  !r c("p.value", "FDR")[2]
+kseaCutoffThreshold:  !r c( 0.1, 0.05)[2]
+kseaMinKinaseCount:   1
+intensityHeatmapRows: 75
 ---
 <!--
-alphaFile: "test-data/alpha_levels.tabular"
+kseaCutoffStatistic:  !r c("p.value", "FDR")[2]
-inputFile: "test-data/test_input_for_anova.tabular"
+kseaCutoffThreshold:  !r c(0.05, 0.1)[1]
-inputFile: "test-data/UT_Phospho_ST_Sites.preproc.tabular"
-inputFile: "test-data/density_failure.preproc_tab.tabular"
+alphaFile:            "test-data/alpha_levels.tabular"
+inputFile:            "test-data/test_input_for_anova.tabular"
+preprocDb:            "test-data/test_input_for_anova.sqlite"
+kseaAppPrepDb:        !r c(":memory:", "test-data/mqppep.sqlite")[2]
+alphaFile:            "test-data/alpha_levels.tabular"
+inputFile:            "test-data/UT_phospho_ST_sites.preproc.tabular"
+preprocDb:            "test-data/UT_phospho_ST_sites.preproc.sqlite"
+kseaAppPrepDb:        !r c(":memory:", "test-data/UT_phospho_ST_sites.ksea.sqlite")[2]
+alphaFile:            "test-data/alpha_levels.tabular"
+inputFile:            "test-data/pY_Sites_NancyDu.txt.ppep_intensities.ppep_map.preproc.tabular"
+preprocDb:            "test-data/pY_Sites_NancyDu.txt.ppep_intensities.ppep_map.preproc.sqlite"
+kseaAppPrepDb:        !r c(":memory:", "test-data/pST_Sites_NancyDu.ksea.sqlite")[2]
+alphaFile:            "test-data/alpha_levels.tabular"
+inputFile:            "test-data/pST_Sites_NancyDu.txt.preproc.tabular"
+preprocDb:            "test-data/pST_Sites_NancyDu.txt.preproc.sqlite"
+kseaAppPrepDb:        !r c(":memory:", "test-data/pST_Sites_NancyDu.ksea.sqlite")[2]
+inputFile:            "test-data/density_failure.preproc_tab.tabular"
+kseaAppPrepDb:        !r c(":memory:", "mqppep.sqlite")[2]
 latex_document: default
 -->
 ```{r setup, include = FALSE}
+#ref for debugging: https://yihui.org/tinytex/r/#debugging
+options(tinytex.verbose = TRUE)
 # ref for parameterizing Rmd document: https://stackoverflow.com/a/37940285
+# ref for top and bottom struts: https://tex.stackexchange.com/a/50355
 knitr::opts_chunk$set(echo = FALSE, fig.dim = c(9, 10))
 # freeze the random number generator so the same results will be produced
 #  from run to run
 set.seed(28571)
+### LIBRARIES
+library(gplots)
+library(DBI)
+library(RSQLite)
+# Suppress "Warning: no DISPLAY variable so Tk is not available"
+suppressWarnings(suppressMessages(library(sqldf)))
+# required but not added to search list:
+# - DBI
+# - RSQLite
+# - ggplot2
+# - knitr
+# - latex2exp
+# - preprocessCore
+# - reshape2
+# - vioplot
 ### CONSTANTS
 const_parfin <- par("fin")
 const_boxplot_fill <- "grey94"
 const_stripchart_cex <- 0.5
 const_stripsmall_cex <-
 sqrt(const_stripchart_cex * const_stripchart_cex / 2)
 const_stripchart_jitter <- 0.3
 const_write_debug_files <- FALSE
-const_table_anchor <- "tbp"
+const_table_anchor_bp <- "bp"
+const_table_anchor_ht <- "ht"
+const_table_anchor_p <- "p"
+const_table_anchor_tbp <- "tbp"
+const_ksea_astrsk_kinases    <- 1
+const_ksea_nonastrsk_kinases <- 2
+const_ksea_all_kinases       <- 3
+const_log10_e <- log10(exp(1))
 ### FUNCTIONS
-#ANOVA filter function
+# from `demo(error.catching)`
-anova_func <- function(x, grouping_factor) {
+##' Catch *and* save both errors and warnings, and in the case of
-x_aov <- aov(as.numeric(x) ~ grouping_factor)
+##' a warning, also keep the computed result.
-pvalue <- summary(x_aov)[[1]][["Pr(>F)"]][1]
+##'
-pvalue
+##' @title tryCatch both warnings (with value) and errors
-}
+##' @param expr an \R expression to evaluate
+##' @return a list with 'value' and 'warning', where
+##'   'value' may be an error caught.
+##' @author Martin Maechler;
+##' Copyright (C) 2010-2012  The R Core Team
+try_catch_w_e <- function(expr) {
+wrn <- NULL
+# warning handler
+w_handler <- function(w) {
+wrn <<- w
+invokeRestart("muffleWarning")
+}
+list(
+value = withCallingHandlers(
+tryCatch(
+expr,
+error = function(e) e
+),
+warning = w_handler
+),
+warning = wrn
+)
+}
 write_debug_file <- function(s) {
 if (const_write_debug_files) {
 s_path <- sprintf("test-data/%s.txt", deparse(substitute(s)))
+print(sprintf("DEBUG writing file %s", spath))
 write.table(
 s,
 file = s_path,
 sep = "\t",
 col.names = TRUE,
 quote = FALSE
 )
 }
 }
-latex_collapsed_vector <- function(collapse_string, v) {
+# ref: http://adv-r.had.co.nz/Environments.html
+# "When creating your own environment, note that you should set its parent
+#   environment to be the empty environment. This ensures you don't
+#   accidentally inherit objects from somewhere else."
+# Caution: this prevents `with(my_env, expr)` from working when `expr`
+#   contains anything from the global environment, even operators!
+#   Hence, `x <- 1; get("x", new_env())` fails by design.
+new_env <- function() {
+new.env(parent = emptyenv())
+}
+### numerical/statistical helper functions
+any_nan <- function(x) {
+!any(x == "NaN")
+}
+# determine standard deviation of quantile to impute
+sd_finite <- function(x) {
+ok <- is.finite(x)
+sd(x[ok])
+}
+anova_func <- function(x, grouping_factor, one_way_f) {
+subject <- data.frame(
+intensity = x
+)
+x_aov <-
+one_way_f(
+formula = intensity ~ grouping_factor,
+data = subject
+)
+pvalue <-
+if (identical(one_way_f, aov))
+summary(x_aov)[[1]][["Pr(>F)"]][1]
+else
+pvalue <- x_aov$p.value
+pvalue
+}
+### LaTeX functions
+latex_collapsed_vector <- function(collapse_string, v, underscore_whack = TRUE) {
+v_sub <- if (underscore_whack) gsub("_", "\\\\_", v) else v
 cat(
 paste0(
-gsub("_", "\\\\_", v),
+v_sub,
 collapse = collapse_string
 )
 )
 }
-latex_itemized_collapsed <- function(collapse_string, v) {
+latex_itemized_collapsed <- function(collapse_string, v, underscore_whack = TRUE) {
 cat("\\begin{itemize}\n\\item ")
-latex_collapsed_vector(collapse_string, v)
+latex_collapsed_vector(collapse_string, v, underscore_whack)
 cat("\n\\end{itemize}\n")
 }
-latex_itemized_list <- function(v) {
+latex_itemized_list <- function(v, underscore_whack = TRUE) {
-latex_itemized_collapsed("\n\\item ", v)
+latex_itemized_collapsed("\n\\item ", v, underscore_whack)
 }
-latex_enumerated_collapsed <- function(collapse_string, v) {
+latex_enumerated_collapsed <- function(collapse_string, v, underscore_whack = TRUE) {
 cat("\\begin{enumerate}\n\\item ")
-latex_collapsed_vector(collapse_string, v)
+latex_collapsed_vector(collapse_string, v, underscore_whack)
 cat("\n\\end{enumerate}\n")
 }
 latex_enumerated_list <- function(v) {
 latex_enumerated_collapsed("\n\\item ", v)
 }
-latex_table_row <- function(v) {
+latex_table_row <- function(v, extra = "", underscore_whack = TRUE) {
-latex_collapsed_vector(" & ", v)
+latex_collapsed_vector(" & ", v, underscore_whack)
+cat(extra)
 cat(" \\\\\n")
 }
 # Use this like print.data.frame, from which it is adapted:
 data_frame_latex <-
 max = NULL,
 # string with justification of each column
 justification = NULL,
 # TRUE to center on page
 centered = TRUE,
-# optional capttion
+# optional caption
 caption = NULL,
 # h(inline); b(bottom); t (top) or p (separate page)
-anchor = "h"
+anchor = "h",
+# set underscore_whack to TRUE to escape underscores
+underscore_whack = TRUE
 ) {
 if (is.null(justification))
 justification <-
 Reduce(
 f = paste,
 n
 ),
 "\n",
 sep = ""
 )
-}
+} else if (n == 0L) {
-else if (n == 0L) {
 cat("0 rows for:\n")
-latex_itemized_list(names(x))
+latex_itemized_list(
-}
+v = names(x),
-else {
+underscore_whack = underscore_whack
+)
+} else {
 if (is.null(max))
 max <- getOption("max.print", 99999L)
 if (!is.finite(max))
 stop("invalid 'max' / getOption(\"max.print\"): ",
 max)
 justification,
 "}\n",
 sep = ""
 )
 )
+# ref: https://tex.stackexchange.com/a/50353
+#   Describes use of \rule{0pt}{3ex}
 if (!is.null(caption))
-cat(" \\hline\\hline\n")
+cat("\\B \\\\ \\hline\\hline\n")
-latex_table_row(colnames(m))
+# ref for top and bottom struts: https://tex.stackexchange.com/a/50355
+latex_table_row(
+v = colnames(m),
+extra = "\\T\\B",
+underscore_whack = underscore_whack
+)
 cat("\\hline\n")
 for (i in seq_len(length(m[, 1]))) {
-latex_table_row(m[i, ])
+latex_table_row(
+v = m[i, ],
+underscore_whack = underscore_whack
+)
 }
 cat(
 paste(
 " \\end{tabular}",
 "\\end{table}",
 n - n0, "rows ]\n")
 }
 invisible(x)
 }
-```
+hypersub <-
+function(s) {
-## Purpose
+hyper <- tolower(s)
+hyper <- gsub("[^a-z0-9]+", "-", hyper)
-Perform imputation of missing values, quantile normalization, and ANOVA.
+hyper <- gsub("[-]+",       "-", hyper)
+hyper <- sub("^[-]",        "",  hyper)
+hyper <- sub("[-]$",        "",  hyper)
+return(hyper)
+}
+subsection_header <-
+function(s) {
+hyper <- hypersub(s)
+cat(
+sprintf(
+"\\hypertarget{%s}\n{\\subsection{%s}\\label{%s}}\n",
+hyper, s, hyper
+)
+)
+}
+subsubsection_header <-
+function(s) {
+hyper <- hypersub(s)
+cat(
+sprintf(
+"\\hypertarget{%s}\n{\\subsubsection{%s}\\label{%s}}\n",
+hyper, s, hyper
+)
+)
+}
+### SQLite functions
+ddl_exec <- function(db, sql) {
+discard <- DBI::dbExecute(conn = db, statement = sql)
+if (FALSE && discard != 0) {
+need_newpage <- TRUE
+if (need_newpage) {
+need_newpage <<- FALSE
+cat("\\newpage\n")
+}
+o_file <- stdout()
+cat("\n\\begin{verbatim}\n")
+cat(sql, file = o_file)
+cat(sprintf("\n%d rows affected by DDL\n", discard), file = o_file)
+cat("\n\\end{verbatim}\n")
+}
+}
+dml_no_rows_exec <- function(db, sql) {
+discard <- DBI::dbExecute(conn = db, statement = sql)
+if (discard != 0) {
+need_newpage <- TRUE
+if (need_newpage) {
+need_newpage <<- FALSE
+cat("\\newpage\n")
+}
+cat("\n\\begin{verbatim}\n")
+o_file <- stdout()
+cat(sql, file = o_file)
+cat(sprintf("\n%d rows affected by DML\n", discard), file = o_file)
+cat("\n\\end{verbatim}\n")
+}
+}
+### KSEA functions and helpers
+# Adapted from KSEAapp::KSEA.Scores to allow retrieval of:
+# - maximum log2(FC)
+ksea_scores <- function(
+# For human data, typically, ksdata = KSEAapp::ksdata
+ksdata,
+# Input data file having columns:
+# - Protein     : abbreviated protein name
+# - Gene        : HUGO gene name
+# - Peptide     : peptide sequence without indications of phosphorylation
+# - Reside.Both : position(s) of phosphorylation within Gene sequence
+#                 - First letter designates AA that is modified
+#                 - Numbers indicate position within Gene
+#                 - Multiple values are separated by semicolons
+#   - p         : p-value
+#   - FC        : fold-change
+px,
+# A binary input of TRUE or FALSE, indicating whether or not to include
+#   NetworKIN predictions
+networkin,
+# A numeric value between 1 and infinity setting the minimum NetworKIN
+#   score (can be left out if networkin = FALSE)
+networkin_cutoff
+) {
+if (length(grep(";", px$Residue.Both)) == 0) {
+# There are no Residue.Both entries having semicolons, so new is
+#   simply px except two columns are renamed and a column is added
+#   for log2(abs(fold-change))
+new <- px
+colnames(new)[c(2, 4)] <- c("SUB_GENE", "SUB_MOD_RSD")
+new$log2_fc <- log2(abs(as.numeric(as.character(new$FC))))
+new <- new[complete.cases(new$log2_fc), ]
+} else {
+# Split each row having semicolons in Residue.Both into rows that are
+#   duplicated in all respects except that each row has a single
+#   member of the set "split-on-semicolon-Residue.Both"
+px_double <- px[grep(";", px$Residue.Both), ]
+residues <- as.character(px_double$Residue.Both)
+residues <- as.matrix(residues, ncol = 1)
+split <- strsplit(residues, split = ";")
+# x gets count of residues in each row,
+#   i.e., 1 + count of semicolons
+x <- sapply(split, length)
+# Here is the set of split rows
+px_single <- data.frame(
+Protein      = rep(px_double$Protein, x),
+Gene         = rep(px_double$Gene,    x),
+Peptide      = rep(px_double$Peptide, x),
+Residue.Both = unlist(split),
+p            = rep(px_double$p,       x),
+FC           = rep(px_double$FC,      x)
+)
+# new first gets the split rows
+new <- px[-grep(";", px$Residue.Both), ]
+# to new, append the rows that didn't need splitting in the first place
+new <- rbind(new, px_single)
+# map Gene         to SUB_GENE
+# map Residue.Both to SUB_MOD_RSD
+colnames(new)[c(2, 4)] <- c("SUB_GENE", "SUB_MOD_RSD")
+# Eliminate any non-positive values to prevent introduction of
+#   infinite or NaN values
+new[(0 <= new$log2_fc), "log2_fc"] <- NA
+# Because of preceding step, there is no need for abs in the next line
+new$log2_fc <- log2(as.numeric(as.character(new$FC)))
+# Convert any illegal values from NaN to NA
+new[is.nan(new$log2_fc), "log2_fc"] <- NA
+# Eliminate rows having missing values (e.g., non-imputed data)
+new <- new[complete.cases(new$log2_fc), ]
+}
+if (networkin == TRUE) {
+# When NetworKIN is true, filter on NetworKIN.cutoff which includes
+#   PhosphoSitePlus data *because its networkin_score is set to Inf*
+ksdata_filtered <- ksdata[grep("[a-z]", ksdata$Source), ]
+ksdata_filtered <- ksdata_filtered[
+(ksdata_filtered$networkin_score >= networkin_cutoff), ]
+} else {
+# Otherwise, simply use PhosphSitePlus rows
+ksdata_filtered <- ksdata[
+grep("PhosphoSitePlus", ksdata$Source), ]
+}
+# Join the two data.frames on common columns SUB_GENE and SUB_MOD_RSD
+#   colnames of ksdata_filtered:
+#     "KINASE" "KIN_ACC_ID" "GENE" "KIN_ORGANISM" "SUBSTRATE" "SUB_GENE_ID"
+#     "SUB_ACC_ID" "SUB_GENE" "SUB_ORGANISM" "SUB_MOD_RSD" "SITE_GRP_ID"
+#     "SITE_...7_AA" "networkin_score" "Source"
+#   colnames of new:
+#     "Protein" "SUB_GENE" "Peptide" "SUB_MOD_RSD" "p" "FC" "log2_fc"
+# Equivalent to:
+#   SELECT a.*. b.Protein, b.Peptide, b.p, b.FC, b.log2_fc
+#     FROM ksdata_filtered a
+#       INNER JOIN new b
+#         ON a.SUB_GENE = b.SUB_GENE
+#           AND a.SUB_MOD_RSD = b.SUB_MOD_RSD
+ksdata_dataset <- base::merge(ksdata_filtered, new)
+#   colnames of ksdata_dataset:
+#     "KINASE"      "KIN_ACC_ID"   "GENE"       "KIN_ORGANISM" "SUBSTRATE"
+#     "SUB_GENE_ID" "SUB_ACC_ID"   "SUB_GENE"   "SUB_ORGANISM" "SUB_MOD_RSD"
+#     "SITE_GRP_ID" "SITE_...7_AA" "networkin_score"  "Source" "Protein"
+#     "Peptide"     "p"            "FC"         "log2_fc" (uniprot_no_isoform)
+# Re-order dataset; prior to accounting for isoforms
+ksdata_dataset <- ksdata_dataset[order(ksdata_dataset$GENE), ]
+# Extract non-isoform accession in UniProtKB
+ksdata_dataset$uniprot_no_isoform <- sapply(
+ksdata_dataset$KIN_ACC_ID,
+function(x) unlist(strsplit(as.character(x), split = "-"))[1]
+)
+# Discard previous results while selecting interesting columns ...
+ksdata_dataset_abbrev <- ksdata_dataset[, c(5, 1, 2, 16:19, 14)]
+# Column names are now:
+#   "GENE"       "SUB_GENE"        "SUB_MOD_RSD"    "Peptide" "p"
+#   "FC" "log2_fc" "Source"
+# Make column names human-readable
+colnames(ksdata_dataset_abbrev) <- c(
+"Kinase.Gene", "Substrate.Gene", "Substrate.Mod", "Peptide", "p",
+"FC", "log2FC", "Source"
+)
+# SELECT * FROM ksdata_dataset_abbrev
+#   ORDER BY Kinase.Gene, Substrate.Gene, Substrate.Mod, p
+ksdata_dataset_abbrev <-
+ksdata_dataset_abbrev[
+order(
+ksdata_dataset_abbrev$Kinase.Gene,
+ksdata_dataset_abbrev$Substrate.Gene,
+ksdata_dataset_abbrev$Substrate.Mod,
+ksdata_dataset_abbrev$p),
+]
+# First aggregation step to account for multiply phosphorylated peptides
+#   and differing peptide sequences; the goal here is to combine results
+#   for all measurements of the same substrate.
+# SELECT  `Kinase.Gene`, `Substrate.Gene`, `Substrate.Mod`,
+#         `Source`, avg(log2FC) AS log2FC
+#   FROM  ksdata_dataset_abbrev
+#   GROUP BY `Kinase.Gene`, `Substrate.Gene`, `Substrate.Mod`,
+#         `Source`
+#   ORDER BY `Kinase.Gene`;
+# in two steps:
+# (1) compute average log_2(fold-change)
+ksdata_dataset_abbrev <- aggregate(
+log2FC ~ Kinase.Gene + Substrate.Gene + Substrate.Mod + Source,
+data = ksdata_dataset_abbrev,
+FUN = mean
+)
+# (2) order by Kinase.Gene
+ksdata_dataset_abbrev <-
+ksdata_dataset_abbrev[order(ksdata_dataset_abbrev$Kinase.Gene), ]
+# SELECT  `Kinase.Gene`, count(*)
+#   FROM  ksdata_dataset_abbrev
+#  GROUP BY `Kinase.Gene`;
+# in two steps:
+# (1) Extract the list of Kinase.Gene names
+kinase_list <- as.vector(ksdata_dataset_abbrev$Kinase.Gene)
+# (2) Convert to a named list of counts of kinases in ksdata_dataset_abrev,
+#   named by Kinase.Gene
+kinase_list <- as.matrix(table(kinase_list))
+# Second aggregation step to account for all substrates per kinase
+# CREATE TABLE mean_fc
+#   AS
+# SELECT  `Kinase.Gene`, avg(log2FC) AS log2FC
+#   FROM  ksdata_dataset_abbrev
+#   GROUP BY `Kinase.Gene`
+mean_fc <- aggregate(
+log2FC ~ Kinase.Gene,
+data = ksdata_dataset_abbrev,
+FUN = mean
+)
+# mean_fc columns: "Kinase.Gene", "log2FC"
+if (FALSE) {
+# I need to re-think this; I was trying to find the most-represented
+#   peptide, but that horse has already left the barn
+# SELECT  `Kinase.Gene`, max(abs(log2FC)) AS log2FC
+#   FROM  ksdata_dataset_abbrev
+#   GROUP BY `Kinase.Gene`
+max_fc <- aggregate(
+log2FC ~ Kinase.Gene,
+data = ksdata_dataset_abbrev,
+FUN = function(r) max(abs(r))
+)
+}
+# Create column 3: mS
+mean_fc$m_s <- mean_fc[, 2]
+# Create column 4: Enrichment
+mean_fc$enrichment <- mean_fc$m_s / abs(mean(new$log2_fc, na.rm = TRUE))
+# Create column 5: m, count of substrates
+mean_fc$m <- kinase_list
+# Create column 6: z-score
+mean_fc$z_score <- (
+(mean_fc$m_s - mean(new$log2_fc, na.rm = TRUE)) *
+sqrt(mean_fc$m)) / sd(new$log2_fc, na.rm = TRUE)
+# Create column 7: p-value, deduced from z-score
+mean_fc$p_value <- pnorm(-abs(mean_fc$z_score))
+# Create column 8: FDR, deduced by Benjamini-Hochberg adustment from p-value
+mean_fc$fdr <- p.adjust(mean_fc$p_value, method = "fdr")
+# Remove log2FC column, which is duplicated as mS
+mean_fc <- mean_fc[order(mean_fc$Kinase.Gene), -2]
+# Correct the column names which we had to hack because of the linter...
+colnames(mean_fc) <- c(
+"Kinase.Gene", "mS", "Enrichment", "m", "z.score", "p.value", "FDR"
+)
+return(mean_fc)
+}
+low_fdr_barplot <- function(
+rslt,
+i_cntrst,
+i,
+a_level,
+b_level,
+fold_change,
+caption
+) {
+rslt_score_list_i <- rslt$score_list[[i]]
+if (!is.null(rslt_score_list_i)) {
+rslt_score_list_i_nrow <- nrow(rslt_score_list_i)
+k <- data.frame(
+contrast = as.integer(i_cntrst),
+a_level = rep.int(a_level, rslt_score_list_i_nrow),
+b_level = rep.int(b_level, rslt_score_list_i_nrow),
+kinase_gene = rslt_score_list_i$Kinase.Gene,
+mean_log2_fc = rslt_score_list_i$mS,
+enrichment = rslt_score_list_i$Enrichment,
+substrate_count = rslt_score_list_i$m,
+z_score = rslt_score_list_i$z.score,
+p_value = rslt_score_list_i$p.value,
+fdr = rslt_score_list_i$FDR
+)
+selector <- switch(
+ksea_cutoff_statistic,
+"FDR" = {
+k$fdr
+},
+"p.value" = {
+k$p_value
+},
+stop(
+sprintf(
+"Unexpected cutoff statistic %s rather than 'FDR' or 'p.value'",
+ksea_cutoff_statistic
+)
+)
+)
+k <- k[selector < ksea_cutoff_threshold, ]
+if (nrow(k) > 1) {
+op <- par(mai = c(1, 1.5, 0.4, 0.4))
+numeric_z_score <- as.numeric(k$z_score)
+z_score_order <- order(numeric_z_score)
+kinase_name <- k$kinase_gene
+long_caption <-
+sprintf(
+"Kinase z-score, %s < %s, %s",
+ksea_cutoff_statistic,
+ksea_cutoff_threshold,
+caption
+)
+my_cex_caption <- 65.0 / max(65.0, nchar(long_caption))
+cat("\n\\clearpage\n")
+barplot(
+height = numeric_z_score[z_score_order],
+border = NA,
+xpd = FALSE,
+cex.names = 1.0,
+cex.axis = 1.0,
+main = long_caption,
+cex.main = my_cex_caption,
+names.arg = kinase_name[z_score_order],
+horiz = TRUE,
+srt = 45,
+las = 1)
+par(op)
+}
+}
+}
+# note that this adds elements to the global variable `ksea_asterisk_hash`
+low_fdr_print <- function(
+rslt,
+i_cntrst,
+i,
+a_level,
+b_level,
+fold_change,
+caption
+) {
+rslt_score_list_i <- rslt$score_list[[i]]
+if (!is.null(rslt_score_list_i)) {
+rslt_score_list_i_nrow <- nrow(rslt_score_list_i)
+k <- contrast_ksea_scores <- data.frame(
+contrast = as.integer(i_cntrst),
+a_level = rep.int(a_level, rslt_score_list_i_nrow),
+b_level = rep.int(b_level, rslt_score_list_i_nrow),
+kinase_gene = rslt_score_list_i$Kinase.Gene,
+mean_log2_fc = rslt_score_list_i$mS,
+enrichment = rslt_score_list_i$Enrichment,
+substrate_count = rslt_score_list_i$m,
+z_score = rslt_score_list_i$z.score,
+p_value = rslt_score_list_i$p.value,
+fdr = rslt_score_list_i$FDR
+)
+selector <- switch(
+ksea_cutoff_statistic,
+"FDR" = {
+k$fdr
+},
+"p.value" = {
+k$p_value
+},
+stop(
+sprintf(
+"Unexpected cutoff statistic %s rather than 'FDR' or 'p.value'",
+ksea_cutoff_statistic
+)
+)
+)
+k <- k[selector < ksea_cutoff_threshold, ]
+# save kinase names to ksea_asterisk_hash
+for (kinase_name in k$kinase_gene) {
+ksea_asterisk_hash[[kinase_name]] <- 1
+}
+db_write_table_overwrite <- (i_cntrst < 2)
+db_write_table_append    <- !db_write_table_overwrite
+RSQLite::dbWriteTable(
+conn = db,
+name = "contrast_ksea_scores",
+value = contrast_ksea_scores,
+append = db_write_table_append
+)
+selector <- switch(
+ksea_cutoff_statistic,
+"FDR" = {
+contrast_ksea_scores$fdr
+},
+"p.value" = {
+contrast_ksea_scores$p_value
+},
+stop(
+sprintf(
+"Unexpected cutoff statistic %s rather than 'FDR' or 'p.value'",
+ksea_cutoff_statistic
+)
+)
+)
+output_df <- contrast_ksea_scores[
+selector < ksea_cutoff_threshold,
+c("kinase_gene", "mean_log2_fc", "enrichment", "substrate_count",
+"z_score", "p_value", "fdr")
+]
+output_order <- with(output_df, order(mean_log2_fc, kinase_gene))
+output_df <- output_df[output_order, ]
+colnames(output_df) <-
+c(
+colnames(output_df)[1],
+colnames(output_df)[2],
+"enrichment",
+"m_s",
+"z_score",
+"p_value",
+"fdr"
+)
+output_df$fdr <- sprintf("%0.4f", output_df$fdr)
+output_df$p_value <- sprintf("%0.2e", output_df$p_value)
+output_df$z_score <- sprintf("%0.2f", output_df$z_score)
+output_df$m_s <- sprintf("%d", output_df$m_s)
+output_df$enrichment <- sprintf("%0.2f", output_df$enrichment)
+output_ncol <- ncol(output_df)
+colnames(output_df) <-
+c(
+"Kinase",
+"\\(\\overline{\\log_2 (|\\text{fold-change}|)}\\)",
+"Enrichment",
+"Substrates",
+"z-score",
+"p-value",
+"FDR"
+)
+selector <- switch(
+ksea_cutoff_statistic,
+"FDR" = {
+rslt$score_list[[i]]$FDR
+},
+"p.value" = {
+rslt$score_list[[i]]$p.value
+},
+stop(
+sprintf(
+"Unexpected cutoff statistic %s rather than 'FDR' or 'p.value'",
+ksea_cutoff_statistic
+)
+)
+)
+if (sum(selector < ksea_cutoff_threshold) > 0) {
+math_caption <- gsub("{", "\\{", caption,      fixed = TRUE)
+math_caption <- gsub("}", "\\}", math_caption, fixed = TRUE)
+data_frame_latex(
+x = output_df,
+justification = "l c c c c c c",
+centered = TRUE,
+caption = sprintf(
+"\\text{%s}, %s < %s",
+math_caption,
+ksea_cutoff_statistic,
+ksea_cutoff_threshold
+),
+anchor = const_table_anchor_p
+)
+} else {
+cat(
+sprintf(
+"\\break
+No kinases had
+\\(\\text{%s}_\\text{enrichment} < %s\\)
+for contrast %s\\hfill\\break\n",
+ksea_cutoff_statistic,
+ksea_cutoff_threshold,
+caption
+)
+)
+}
+}
+}
+# create_breaks is a helper for ksea_heatmap
+create_breaks <- function(merged_scores) {
+if (min(merged_scores, na.rm = TRUE) < -1.6) {
+breaks_neg <- seq(-1.6, 0, length.out = 30)
+breaks_neg <-
+append(
+seq(min(merged_scores, na.rm = TRUE), -1.6, length.out = 10),
+breaks_neg
+)
+breaks_neg <- sort(unique(breaks_neg))
+} else {
+breaks_neg <- seq(-1.6, 0, length.out = 30)
+}
+if (max(merged_scores, na.rm = TRUE) > 1.6) {
+breaks_pos <- seq(0, 1.6, length.out = 30)
+breaks_pos <-
+append(
+breaks_pos,
+seq(1.6, max(merged_scores, na.rm = TRUE),
+length.out = 10)
+)
+breaks_pos <- sort(unique(breaks_pos))
+} else {
+breaks_pos <- seq(0, 1.6, length.out = 30)
+}
+breaks_all <- unique(append(breaks_neg, breaks_pos))
+mycol_neg <-
+gplots::colorpanel(n = length(breaks_neg),
+low = "blue",
+high = "white")
+mycol_pos <-
+gplots::colorpanel(n = length(breaks_pos) - 1,
+low = "white",
+high = "red")
+mycol <- unique(append(mycol_neg, mycol_pos))
+color_breaks <- list(breaks_all, mycol)
+return(color_breaks)
+}
+# draw_kseaapp_summary_heatmap is a helper function for ksea_heatmap
+draw_kseaapp_summary_heatmap <- function(
+x,
+sample_cluster,
+merged_asterisk,
+my_cex_row,
+color_breaks,
+margins,
+...
+) {
+merged_scores <- x
+if (!is.matrix(x)) {
+cat(
+paste0(
+"No plot because \\texttt{typeof(x)} is '",
+typeof(x),
+"' rather than 'matrix'.\n\n"
+)
+)
+} else if (nrow(x) < 2) {
+cat("No plot because matrix x has ", nrow(x), " rows.\n\n")
+cat("\\begin{verbatim}\n")
+str(x)
+cat("\\end{verbatim}\n")
+} else if (ncol(x) < 2) {
+cat("No plot because matrix x has ", ncol(x), " columns.\n\n")
+cat("\\begin{verbatim}\n")
+str(x)
+cat("\\end{verbatim}\n")
+} else {
+gplots::heatmap.2(
+x            = merged_scores,
+Colv         = sample_cluster,
+scale        = "none",
+cellnote     = merged_asterisk,
+notecol      = "white",
+cexCol       = 0.9,
+# Heuristically assign size of row labels
+cexRow       = min(1.0, ((3 * my_cex_row) ^ 1.7) / 2.25),
+srtCol       = 45,
+srtRow       = 45,
+notecex      = 3 * my_cex_row,
+col          = color_breaks[[2]],
+density.info = "none",
+trace        = "none",
+breaks       = color_breaks[[1]],
+lmat         = rbind(c(0, 3), c(2, 1), c(0, 4)),
+lhei         = c(0.4, 8.0, 1.1),
+lwid         = c(0.5, 3),
+key          = FALSE,
+margins      = margins,
+...
+)
+}
+}
+# Adapted from KSEAapp::KSEA.Heatmap
+ksea_heatmap <- function(
+# the data frame outputs from the KSEA.Scores() function, in list format
+score_list,
+# a character vector of all the sample names for heatmap annotation:
+# - the names must be in the same order as the data in score_list
+# - please avoid long names, as they may get cropped in the final image
+sample_labels,
+# character string of either "p.value" or "FDR" indicating the data column
+#   to use for marking statistically significant scores
+stats,
+# a numeric value between 0 and infinity indicating the min. number of
+#   substrates a kinase must have to be included in the heatmap
+m_cutoff,
+# a numeric value between 0 and 1 indicating the p-value/FDR cutoff
+#   for indicating significant kinases in the heatmap
+p_cutoff =
+stop("argument 'p_cutoff' is required for function 'ksea_heatmap'"),
+# a binary input of TRUE or FALSE, indicating whether or not to perform
+#   hierarchical clustering of the sample columns
+sample_cluster,
+# a binary input of TRUE or FALSE, indicating whether or not to export
+#   the heatmap as a .png image into the working directory
+export = FALSE,
+# bottom and right margins; adjust as needed if contrast names are too long
+margins = c(6, 20),
+# print which kinases?
+# - Mandatory argument, must be one of const_ksea_.*_kinases
+which_kinases,
+# additional arguments to gplots::heatmap.2, such as:
+# - main: main title of plot
+# - xlab: x-axis label
+# - ylab: y-axis label
+...
+) {
+filter_m <- function(dataset, m_cutoff) {
+filtered <- dataset[(dataset$m >= m_cutoff), ]
+return(filtered)
+}
+score_list_m <- lapply(score_list, function(...) filter_m(..., m_cutoff))
+for (i in seq_len(length(score_list_m))) {
+names <- colnames(score_list_m[[i]])[c(2:7)]
+colnames(score_list_m[[i]])[c(2:7)] <-
+paste(names, i, sep = ".")
+}
+master <-
+Reduce(
+f = function(...) {
+base::merge(..., by = "Kinase.Gene", all = FALSE)
+},
+x = score_list_m
+)
+row.names(master) <- master$Kinase.Gene
+columns <- as.character(colnames(master))
+merged_scores <- as.matrix(master[, grep("z.score", columns), drop = FALSE])
+colnames(merged_scores) <- sample_labels
+merged_stats <- as.matrix(master[, grep(stats, columns)])
+asterisk <- function(mtrx, p_cutoff) {
+new <- data.frame()
+for (i in seq_len(nrow(mtrx))) {
+for (j in seq_len(ncol(mtrx))) {
+my_value <- mtrx[i, j]
+if (!is.na(my_value) && my_value < p_cutoff) {
+new[i, j] <- "*"
+} else {
+new[i, j] <- ""
+}
+}
+}
+return(new)
+}
+merged_asterisk <- as.matrix(asterisk(merged_stats, p_cutoff))
+# begin hack to print only significant rows
+asterisk_rows <- rowSums(merged_asterisk == "*") > 0
+all_rows <- rownames(merged_stats)
+names(asterisk_rows) <- all_rows
+non_asterisk_rows <- names(asterisk_rows[asterisk_rows == FALSE])
+asterisk_rows <- names(asterisk_rows[asterisk_rows == TRUE])
+merged_scores_asterisk <- merged_scores[names(asterisk_rows), ]
+merged_scores_non_asterisk <- merged_scores[names(non_asterisk_rows), ]
+# end hack to print only significant rows
+row_list <- list()
+row_list[[const_ksea_astrsk_kinases]] <- asterisk_rows
+row_list[[const_ksea_all_kinases]] <- all_rows
+row_list[[const_ksea_nonastrsk_kinases]] <- non_asterisk_rows
+i <- which_kinases
+my_row_names <- row_list[[i]]
+scrs <- merged_scores[my_row_names, ]
+stts <- merged_stats[my_row_names, ]
+merged_asterisk <- as.matrix(asterisk(stts, p_cutoff))
+color_breaks <- create_breaks(scrs)
+plot_height <- nrow(scrs) ^ 0.55
+plot_width <- ncol(scrs) ^ 0.7
+my_cex_row <- 0.25 * 16 / plot_height
+if (export == "TRUE") {
+png(
+"KSEA.Merged.Heatmap.png",
+width = plot_width * 300,
+height = 2 * plot_height * 300,
+res = 300,
+pointsize = 14
+)
+}
+draw_kseaapp_summary_heatmap(
+x               = scrs,
+sample_cluster  = sample_cluster,
+merged_asterisk = merged_asterisk,
+my_cex_row      = my_cex_row,
+color_breaks    = color_breaks,
+margins         = margins
+)
+if (export == "TRUE") {
+dev.off()
+}
+return(my_row_names)
+}
+# helper for heatmaps of phosphopeptide intensities
+draw_intensity_heatmap <-
+function(
+m,                              # matrix with rownames already formatted
+cutoff,                         # cutoff used by hm_heading_function
+hm_heading_function,            # construct and cat heading from m and cutoff
+hm_main_title,                  # main title for plot (drawn below heading)
+suppress_row_dendrogram = TRUE, # set to false to show dendrogram
+max_peptide_count               # experimental:
+= intensity_hm_rows,      #   values of 50 and 75 worked well
+...                             # passthru parameters for heatmap
+) {
+peptide_count <- 0
+# emit the heading for the heatmap
+if (hm_heading_function(m, cutoff)) {
+peptide_count <- min(max_peptide_count, nrow(m))
+if (nrow(m) > 1) {
+m_margin <- m[peptide_count:1, ]
+# Margin setting was heuristically derived
+margins <-
+c(0.5, # col
+max(80, sqrt(nchar(rownames(m_margin)))) * 5 / 16  # row
+)
+}
+if (nrow(m) > 1) {
+tryCatch(
+{
+old_oma <- par("oma")
+par(cex.main = 0.6)
+# Heuristically determined character size adjustment formula
+char_contractor <-
+250000 / (
+max(4500, (nchar(rownames(m_margin)))^2) * intensity_hm_rows
+)
+heatmap(
+m[peptide_count:1, ],
+Rowv = if (suppress_row_dendrogram) NA else NULL,
+Colv = NA,
+cexRow = char_contractor,
+cexCol = char_contractor * 50 / max_peptide_count,
+scale = "row",
+margins = margins,
+main =
+"Unimputed, unnormalized log(intensities)",
+xlab = "",
+las = 1,
+...
+)
+},
+error = function(e) {
+cat(
+sprintf(
+"\nCould not draw heatmap, possibly because of too many missing values.  Internal message: %s\n",
+e$message
+)
+)
+},
+finally = par(old_oma)
+)
+}
+}
+return(peptide_count)
+}
+```
+```{r, echo = FALSE, fig.dim = c(9, 10), results = 'asis'}
+cat("\\listoftables\n")
+```
+# Purpose
+To perform for phosphopeptides:
+- imputation of missing values,
+- quantile normalization,
+- ANOVA (using the R stats::`r params$oneWayManyCategories` function), and
+- KSEA (Kinase-Substrate Enrichment Analysis) using code adapted from the CRAN `KSEAapp` package to search for kinase substrates from the following databases:
+- PhosphoSitesPlus [https://www.phosphosite.org](https://www.phosphosite.org)
+- The Human Proteome Database [http://hprd.org](http://hprd.org)
+- NetworKIN [http://networkin.science/](http://networkin.science/)
+- Phosida [http://pegasus.biochem.mpg.de/phosida/help/motifs.aspx](http://pegasus.biochem.mpg.de/phosida/help/motifs.aspx)
 ```{r include = FALSE}
+### GLOBAL VARIABLES
+# parameters for KSEA
+ksea_cutoff_statistic <- params$kseaCutoffStatistic
+ksea_cutoff_threshold <- params$kseaCutoffThreshold
+ksea_min_kinase_count <- params$kseaMinKinaseCount
+ksea_heatmap_titles <- list()
+ksea_heatmap_titles[[const_ksea_astrsk_kinases]] <-
+sprintf(
+"Summary for all kinases enriched in one or more contrasts at %s < %s",
+ksea_cutoff_statistic,
+ksea_cutoff_threshold
+)
+ksea_heatmap_titles[[const_ksea_all_kinases]] <-
+"Summary figure for all contrasts and all kinases"
+ksea_heatmap_titles[[const_ksea_nonastrsk_kinases]] <-
+sprintf(
+"Summary for all kinases not enriched at %s < %s in any contrast",
+ksea_cutoff_statistic,
+ksea_cutoff_threshold
+)
+# hash to hold names of significantly enriched kinases
+ksea_asterisk_hash <- new_env()
+# READ PARAMETERS (mostly)
+intensity_hm_rows <- params$intensityHeatmapRows
 # Input Filename
 input_file <- params$inputFile
 # First data column - ideally, this could be detected via regexSampleNames,
 #   but for now leave it as is.
 }
 # False discovery rate adjustment for ANOVA
 #  Since pY abundance is low, set to 0.10 and 0.20 in addition to 0.05
 val_fdr <-
-read.table(file = params$alphaFile, sep = "\t", header = F, quote = "")
+read.table(file = params$alphaFile, sep = "\t", header = FALSE, quote = "")
 if (
 ncol(val_fdr) != 1 ||
 sum(!is.numeric(val_fdr[, 1])) ||
 sum(val_fdr[, 1] < 0) ||
 val_fdr <- val_fdr[, 1]
 #Imputed Data filename
 imputed_data_filename <- params$imputedDataFilename
 imp_qn_lt_data_filenm <- params$imputedQNLTDataFile
+anova_ksea_mtdt_file  <- params$anovaKseaMetadata
-#ANOVA data filename
 ```
 ```{r echo = FALSE}
 # Imputation method, should be one of
 #   "random", "group-median", "median", or "mean"
 # Regular expression of Sample Names, e.g., "\\.(\\d+)[A-Z]$"
 regex_sample_names <- params$regexSampleNames
 # Regular expression to extract Sample Grouping from Sample Name;
-#   if error occurs, compare sample_factor_levels and sample_name_matches
+#   if error occurs, compare sample_treatment_levels vs. sample_name_matches
 #   to see if groupings/pairs line up
 #   e.g., "(\\d+)"
 regex_sample_grouping <- params$regexSampleGrouping
+one_way_all_categories_fname <- params$oneWayManyCategories
+one_way_all_categories <- try_catch_w_e(
+match.fun(one_way_all_categories_fname))
+if (!is.function(one_way_all_categories$value)) {
+write("fatal error for parameter oneWayManyCategories:", stderr())
+write(one_way_all_categories$value$message,             stderr())
+if (sys.nframe() > 0) quit(save = "no", status = 1)
+stop("Cannot continue. Goodbye.")
+}
+one_way_all_categories <- one_way_all_categories$value
+one_way_two_categories_fname <- params$oneWayManyCategories
+one_way_two_categories <- try_catch_w_e(
+match.fun(one_way_two_categories_fname))
+if (!is.function(one_way_two_categories$value)) {
+cat("fatal error for parameter oneWayTwoCategories: \n")
+cat(one_way_two_categories$value$message, fill = TRUE)
+if (sys.nframe() > 0) quit(save = "no", status = 1)
+stop("Cannot continue. Goodbye.")
+}
+one_way_two_categories <- one_way_two_categories$value
+preproc_db       <- params$preprocDb
+ksea_app_prep_db <- params$kseaAppPrepDb
+result <- file.copy(
+from      = preproc_db,
+to        = ksea_app_prep_db,
+overwrite = TRUE
+)
+if (!result) {
+write(
+sprintf(
+"fatal error copying initial database '%s' to output '%s'",
+preproc_db,
+ksea_app_prep_db,
+),
+stderr()
+)
+if (sys.nframe() > 0) quit(save = "no", status = 1)
+stop("Cannot continue. Goodbye.")
+}
 ```
 ```{r echo = FALSE}
 ### READ DATA
-library(data.table)
 # read.table reads a file in table format and creates a data frame from it.
 #   - note that `quote = ""` means that quotation marks are treated literally.
 full_data <- read.table(
 file = input_file,
 sep = "\t",
-header = T,
+header = TRUE,
 quote = "",
 check.names = FALSE
 )
 ```
-## Extract Sample Names and Factor Levels
+# Extract Sample Names and Treatment Levels
-Column names parsed from input file are shown in Table 1; sample names and factor levels, in Table 2.
+Column names parsed from input file are shown in Table 1; sample names and treatment levels, in Table 2.
 ```{r echo = FALSE, results = 'asis'}
 data_column_indices <- grep(first_data_column, names(full_data), perl = TRUE)
 }
 cat(
 sprintf(
 paste(
-"\n\nPeptide-intensity data for each sample is",
+"\n\nThe input data file has peptide-intensity data for each sample",
 "in one of columns %d through %d.\n\n"
 ),
 min(data_column_indices),
 max(data_column_indices)
 )
 )
 # Write column names as a LaTeX enumerated list.
 column_name_df <- data.frame(
 column = seq_len(length(colnames(full_data))),
-name = colnames(full_data)
+name = paste0("\\verb@", colnames(full_data), "@")
 )
 data_frame_latex(
 x = column_name_df,
 justification = "l l",
 centered = TRUE,
-caption = "Input data column name",
+caption = "Input data column names",
-anchor = const_table_anchor
+anchor = const_table_anchor_bp,
+underscore_whack = FALSE
 )
 ```
 ```{r echo = FALSE, results = 'asis'}
-#```{r echo = FALSE, results = 'asis'}
 quant_data <- full_data[first_data_column:length(full_data)]
 quant_data[quant_data == 0] <- NA
-rownames(quant_data) <- full_data$Phosphopeptide
+rownames(quant_data) <- rownames(full_data) <- full_data$Phosphopeptide
-# Get factors -> group replicates (as indicated by terminal letter)
+# Extract factors and trt-replicates using regular expressions.
-#   by the preceding digits;
+# Typically:
-# Assuming that regex_sample_names <- "\\.(\\d+)[A-Z]$"
+#   regex_sample_names    is "\\.\\d+[A-Z]$"
-#   get factors ->
+#   regex_sample_grouping is "\\d+"
-#      group runs (samples) by ignoring terminal [A-Z] in sample names
+# This would distinguish trt-replicates by terminal letter [A-Z]
-#   e.g.
+#   in sample names and group them into trts by the preceding digits.
+#   e.g.:
 #     group .1A .1B .1C into group 1;
 #     group .2A .2B .2C, into group 2;
 #     etc.
 m <- regexpr(regex_sample_names, colnames(quant_data), perl = TRUE)
 sample_name_matches <- regmatches(names(quant_data), m)
 colnames(quant_data) <- sample_name_matches
 write_debug_file(quant_data)
-m2 <- regexpr(regex_sample_grouping, sample_name_matches, perl = TRUE)
+rx_match <- regexpr(regex_sample_grouping, sample_name_matches, perl = TRUE)
-sample_factor_levels <- as.factor(regmatches(sample_name_matches, m2))
+sample_treatment_levels <- as.factor(regmatches(sample_name_matches, rx_match))
 number_of_samples <- length(sample_name_matches)
-sample_factor_df <- data.frame(
+sample_treatment_df <- data.frame(
-sample = sample_name_matches,
+level = sample_treatment_levels,
-level = sample_factor_levels
+sample = sample_name_matches
 )
 data_frame_latex(
-x = sample_factor_df,
+x = sample_treatment_df,
-justification = "c c",
+justification = "rp{0.2\\linewidth} lp{0.3\\linewidth}",
 centered = TRUE,
-caption = "Factor level",
+caption = "Treatment levels",
-anchor = const_table_anchor
+anchor = const_table_anchor_tbp,
-)
+underscore_whack = FALSE
-```
+)
+```
 ```{r echo = FALSE, results = 'asis'}
 cat("\\newpage\n")
 ```
-### Are the log-transformed sample distributions similar?
+## Are the log-transformed sample distributions similar?
 ```{r echo = FALSE, fig.dim = c(9, 5.5), results = 'asis'}
 quant_data[quant_data == 0] <- NA  #replace 0 with NA
 quant_data_log <- log10(quant_data)
 cat("\n\n\n")
 ```
 ```{r echo = FALSE, fig.align = "left", fig.dim = c(9, 4), results = 'asis'}
-library(ggplot2)
 if (nrow(quant_data_log) > 1) {
 quant_data_log_stack <- stack(quant_data_log)
-ggplot(
+ggplot2::ggplot(quant_data_log_stack, ggplot2::aes(x = values)) +
-quant_data_log_stack,
+ggplot2::xlab(latex2exp::TeX("$log_{10}$(peptide intensity)")) +
-aes(x = values)
+ggplot2::ylab("Probability density") +
-) + xlab(latex2exp::TeX("$log_{10}$(peptide intensity)")) +
+ggplot2::geom_density(ggplot2::aes(group = ind, colour = ind), na.rm = TRUE)
-ylab("Probability density") +
-geom_density(
-aes(group = ind, colour = ind),
-na.rm = TRUE
-)
 } else {
 cat("No density plot because there are too few peptides.\n\n")
 }
 ```
-### Globally, are peptide intensities are approximately unimodal?
+## Globally, are peptide intensities are approximately unimodal?
 <!--
 # bquote could be used as an alternative to latex2exp::TeX below particularly
 #   and when plotting math expressions generally, at the expense of mastering
 #   another syntax, which hardly seems worthwhile when I need to use TeX
 main = latex2exp::TeX("Frequency vs. $log_{10}$(peptide intensity)"),
 xlab = latex2exp::TeX("$log_{10}$(peptide intensity)")
 )
 ```
-### Distribution of standard deviations of $log_{10}(\text{intensity})$, ignoring missing values:
+## Distribution of standard deviations of $log_{10}(\text{intensity})$, ignoring missing values
 ```{r echo = FALSE, fig.align = "left", fig.dim = c(9, 5), results = 'asis'}
 # determine quantile
 q1 <- quantile(logvalues, probs = mean_percentile)[1]
-# determine standard deviation of quantile to impute
-sd_finite <- function(x) {
-ok <- is.finite(x)
-sd(x[ok])
-}
 # 1 = row of matrix (ie, phosphopeptide)
 sds <- apply(quant_data_log, 1, sd_finite)
 if (sum(!is.na(sds)) > 2) {
 plot(
-density(sds, na.rm = T)
+density(sds, na.rm = TRUE)
 , main = "Smoothed estimated probability density vs. std. deviation"
 , sub = "(probability estimation made with Gaussian smoothing)"
 , ylab = "Probability density"
 )
 } else {
 ```{r echo = FALSE}
 # prep for trt-median based imputation
 ```
-## Impute Missing Values
+# Impute Missing Values
 ```{r echo = FALSE}
-imp_smry_potential_before <- length(logvalues)
+imp_smry_pot_peptides_before <- nrow(quant_data_log)
-imp_smry_missing_before   <- number_to_impute
+imp_smry_missing_values_before   <- number_to_impute
 imp_smry_pct_missing      <- pct_missing_values
 ```
 ```{r echo = FALSE}
 #Determine number of cells to impute
 ```
 ```{r echo = FALSE}
-#Impute data
-quant_data_imp <- quant_data
 # Identify which values are missing and need to be imputed
-ind <- which(is.na(quant_data_imp), arr.ind = TRUE)
+ind <- which(is.na(quant_data), arr.ind = TRUE)
 ```
 ```{r echo = FALSE, results = 'asis'}
 # Apply imputation
 switch(
 imputation_method
 , "group-median" = {
+quant_data_imp <- quant_data
 imputation_method_description <-
 paste("Substitute missing value with",
-"median peptide intensity for sample group.\n"
+"median peptide-intensity for sample group.\n"
 )
-sample_level_integers <- as.integer(sample_factor_levels)
+sample_level_integers <- as.integer(sample_treatment_levels)
-for (i in seq_len(length(levels(sample_factor_levels)))) {
+# Take the accurate ln(x+1) because the data are log-normally distributed
+#   and because median can involve an average of two measurements.
+quant_data_imp <- log1p(quant_data_imp)
+for (i in seq_len(length(levels(sample_treatment_levels)))) {
+# Determine the columns for this factor-level
 level_cols <- i == sample_level_integers
-ind <- which(is.na(quant_data_imp[, level_cols]), arr.ind = TRUE)
+# Extract those columns
-quant_data_imp[ind, level_cols] <-
+lvlsbst <- quant_data_imp[, level_cols, drop = FALSE]
-apply(quant_data_imp[, level_cols], 1, median, na.rm = T)[ind[, 1]]
+# assign to ind the row-column pairs corresponding to each NA
+ind <- which(is.na(lvlsbst), arr.ind = TRUE)
+# No group-median exists if there is only one sample
+#   a given ppep has no measurement; otherwise, proceed.
+if (ncol(lvlsbst) > 1) {
+the_centers <-
+apply(lvlsbst, 1, median, na.rm = TRUE)
+for (j in seq_len(nrow(lvlsbst))) {
+for (k in seq_len(ncol(lvlsbst))) {
+if (is.na(lvlsbst[j, k])) {
+lvlsbst[j, k] <- the_centers[j]
+}
+}
+}
+quant_data_imp[, level_cols] <- lvlsbst
+}
 }
+# Take the accurate e^x - 1 to match scaling of original input.
+quant_data_imp <- round(expm1(quant_data_imp_ln <- quant_data_imp))
 good_rows <- !is.na(rowMeans(quant_data_imp))
 }
 , "median" = {
+quant_data_imp <- quant_data
 imputation_method_description <-
 paste("Substitute missing value with",
-"median peptide intensity across all sample classes.\n"
+"median peptide-intensity across all sample classes.\n"
 )
-quant_data_imp[ind] <- apply(quant_data_imp, 1, median, na.rm = T)[ind[, 1]]
+# Take the accurate ln(x+1) because the data are log-normally distributed
-good_rows <- !is.na(rowMeans(quant_data_imp))
+#   and because median can involve an average of two measurements.
+quant_data_imp <- log1p(quant_data_imp)
+quant_data_imp[ind] <- apply(quant_data_imp, 1, median, na.rm = TRUE)[ind[, 1]]
+# Take the accurate e^x - 1 to match scaling of original input.
+quant_data_imp <- round(expm1(quant_data_imp_ln <- quant_data_imp))
+good_rows <- !is.nan(rowMeans(quant_data_imp))
 }
 , "mean" = {
+quant_data_imp <- quant_data
 imputation_method_description <-
 paste("Substitute missing value with",
-"mean peptide intensity across all sample classes.\n"
+"geometric-mean peptide intensity across all sample classes.\n"
 )
-quant_data_imp[ind] <- apply(quant_data_imp, 1, mean, na.rm = T)[ind[, 1]]
+# Take the accurate ln(x+1) because the data are log-normally distributed,
-good_rows <- !is.na(rowMeans(quant_data_imp))
+#   so arguments to mean should be previously transformed.
+#   this will have to be
+quant_data_imp <- log1p(quant_data_imp)
+# Assign to NA cells the mean for the row
+quant_data_imp[ind] <- apply(quant_data_imp, 1, mean, na.rm = TRUE)[ind[, 1]]
+# Take the accurate e^x - 1 to match scaling of original input.
+quant_data_imp <- round(expm1(quant_data_imp_ln <- quant_data_imp))
+good_rows <- !is.nan(rowMeans(quant_data_imp))
 }
 , "random" = {
-m1 <- median(sds, na.rm = T) * sd_percentile #sd to be used is the median sd
+quant_data_imp <- quant_data
+m1 <- median(sds, na.rm = TRUE) * sd_percentile #sd to be used is the median sd
 # If you want results to be reproducible, you will want to call
 #   base::set.seed before calling stats::rnorm
 imputation_method_description <-
 paste("Substitute each missing value with random intensity",
 sprintf(
 100 * mean_percentile))
 cat(sprintf("sd_percentile (from input parameter) is %0.2f\n\n",
 sd_percentile))
 quant_data_imp[ind] <-
 10 ^ rnorm(number_to_impute, mean = q1, sd = m1)
-good_rows <- !is.na(rowMeans(quant_data_imp))
+quant_data_imp_ln <- log1p(quant_data_imp)
+good_rows <- !is.nan(rowMeans(quant_data_imp))
 }
 )
+quant_data_imp_log10 <- quant_data_imp_ln * const_log10_e
 if (length(good_rows) < 1) {
 print("ERROR: Cannot impute data; there are no good rows!")
 return(-1)
 }
 cat("\\quad\n\nImputation method:\n\n\n", imputation_method_description)
 ```
 ```{r echo = FALSE}
-imp_smry_potential_after <- sum(good_rows)
+imp_smry_pot_peptides_after <- sum(good_rows)
 imp_smry_rejected_after  <- sum(!good_rows)
-imp_smry_missing_after   <- sum(is.na(quant_data_imp[good_rows, ]))
+imp_smry_missing_values_after   <- sum(is.na(quant_data_imp[good_rows, ]))
 ```
 ```{r echo = FALSE, results = 'asis'}
 # ref: http://www1.maths.leeds.ac.uk/latex/TableHelp1.pdf
-tabular_lines <- paste(
+tabular_lines_fmt <- paste(
 "\\begin{table}[hb]", # h(inline); b(bottom); t (top) or p (separate page)
 " \\caption{Imputation Results}",
 " \\centering", # \centering centers the table on the page
 " \\begin{tabular}{l c c c}",
 "  \\hline\\hline",
 "\\end{table}",
 sep = "\n"
 )
 tabular_lines <-
 sprintf(
-tabular_lines,
+tabular_lines_fmt,
-imp_smry_potential_before,
+imp_smry_pot_peptides_before,
-imp_smry_missing_before,
+imp_smry_missing_values_before,
 imp_smry_pct_missing,
 "%",
-imp_smry_potential_after,
+imp_smry_pot_peptides_after,
-imp_smry_missing_after,
+imp_smry_missing_values_after,
 imp_smry_rejected_after
 )
 cat(tabular_lines)
 ```
 ```{r echo = FALSE}
 # Zap rows where imputation was ineffective
 full_data         <- full_data        [good_rows, ]
 quant_data        <- quant_data       [good_rows, ]
+quant_data_imp <- quant_data_imp[good_rows, ]
 write_debug_file(quant_data_imp)
-quant_data_imp <- quant_data_imp[good_rows, ]
 quant_data_imp_good_rows <- quant_data_imp
 write_debug_file(quant_data_imp_good_rows)
 ```
 as.matrix(
 log10(quant_data_imp[!is.na(quant_data)])
 )
 )
-if (sum(!is.na(d_combined)) > 2 && sum(!is.na(d_original)) > 2) {
+if (sum(!is.na(d_original)) > 2) {
+d_original <- density(d_original)
+} else {
+can_plot_before_after_imp <- FALSE
+}
+if (can_plot_before_after_imp && sum(is.na(d_combined)) < 1) {
 d_combined <- density(d_combined)
-d_original <- density(d_original)
 } else {
 can_plot_before_after_imp <- FALSE
 }
 if (sum(is.na(quant_data)) > 0) {
 as.numeric(
 as.matrix(
 log10(quant_data_imp[is.na(quant_data)])
 )
 )
-if (sum(!is.na(d_combined)) > 2) {
+if (can_plot_before_after_imp && sum(is.na(d_imputed)) < 1) {
 d_imputed <- (density(d_imputed))
 } else {
 can_plot_before_after_imp <- FALSE
 }
 } else {
 }
 ```
 ```{r echo = FALSE, fig.dim = c(9, 5.5), results = 'asis'}
+zero_sd_rownames <-
+rownames(quant_data_imp)[
+is.na((apply(quant_data_imp, 1, sd, na.rm = TRUE)) == 0)
+]
+if (length(zero_sd_rownames) >= nrow(quant_data_imp)) {
+stop("All peptides have zero standard deviation.  Cannot continue.")
+}
+if (length(zero_sd_rownames) > 0) {
+cat(
+sprintf("%d peptides with zero variance were removed from statistical consideration",
+length(zero_sd_rownames)
+)
+)
+zap_named_rows <- function(df, nms) {
+return(df[!(row.names(df) %in% nms), ])
+}
+quant_data_imp <- zap_named_rows(quant_data_imp, zero_sd_rownames)
+quant_data     <- zap_named_rows(quant_data,     zero_sd_rownames)
+full_data      <- zap_named_rows(full_data,      zero_sd_rownames)
+}
 if (sum(is.na(quant_data)) > 0) {
 cat("\\leavevmode\\newpage\n")
 # data visualization
 old_par <- par(
 mai = par("mai") + c(0.5, 0, 0, 0)
 )
+# Copy quant data to x
 x <- quant_data
-x <- blue_dots <- x / x
+# x gets to have values of:
-blue_dots <- log10(blue_dots * quant_data)
+#  - NA for observed values
+#  - 1 for missing values
 x[is.na(x)] <- 0
-x[x == 1] <- NA
+x[x > 1] <- NA
 x[x == 0] <- 1
-quant_data_imp_log10 <-
+# Log-transform imputed data
-log10(quant_data_imp)
+# update variable because rows may have been eliminated from quant_data_imp
+quant_data_imp_log10 <- log10(quant_data_imp)
 write_debug_file(quant_data_imp_log10)
+# Set blue_dots to log of quant data or NA for NA quant data
+blue_dots <- log10(quant_data)
+# Set red_dots to log of imputed data or NA for observed quant data
 red_dots <- quant_data_imp_log10 * x
 count_red <- sum(!is.na(red_dots))
 count_blue <- sum(!is.na(blue_dots))
 ylim_save <- ylim <- c(
 min(red_dots, blue_dots, na.rm = TRUE),
 max(red_dots, blue_dots, na.rm = TRUE)
 boxplot(
 blue_dots
 , las = 1 # "always horizontal"
 , col = const_boxplot_fill
 , ylim = ylim
-, main = "Peptide intensities before and after imputation"
+, main = "Peptide intensities after eliminating unusable peptides"
 , sub = boxplot_sub
 , xlab = "Sample"
 , ylab = latex2exp::TeX("$log_{10}$(peptide intensity)")
 )
 col = "red",                # Color of the symbol
 vertical = TRUE,            # Vertical mode
 add = TRUE                  # Add it over
 )
-} else {
+}
-# violin plot
+if (TRUE) {
+# show measured values in blue on left half-violin plot
 cat("\\leavevmode\n\\quad\n\n\\quad\n\n")
 vioplot::vioplot(
 x = lapply(blue_dots, function(x) x[!is.na(x)]),
 col = "lightblue1",
 side = "left",
 main = "Distributions of observed and imputed data",
 sub = "Light blue = observed data; Pink = imputed data",
 xlab = "Sample",
 ylab = latex2exp::TeX("$log_{10}$(peptide intensity)")
 )
+red_violins <- lapply(red_dots, function(x) x[!is.na(x)])
+cols_to_delete <- c()
+for (ix in seq_len(length(red_violins))) {
+if (length(red_violins[[ix]]) < 1) {
+cols_to_delete <- c(cols_to_delete, ix)
+}
+}
+# destroy any unimputable columns
+if (!is.null(cols_to_delete)) {
+red_violins <- red_violins[-cols_to_delete]
+}
+# plot imputed values in red on right half-violin plot
 vioplot::vioplot(
-x = lapply(red_dots, function(x) x[!is.na(x)]),
+x = red_violins,
 col = "lightpink1",
 side = "right",
 plotCentre = "line",
-add = T
+add = TRUE
 )
 }
 par(old_par)
 # density plot
 cat("\\leavevmode\n\n\n\n\n\n\n")
 if (can_plot_before_after_imp) {
-ylim <- c(0, max(d_combined$y, d_original$y, d_imputed$y))
+ylim <- c(
+0,
+max(
+if (is.list(d_combined)) d_combined$y else d_combined,
+if (is.list(d_original)) d_original$y else d_original,
+if (is.list(d_imputed)) d_imputed$y else d_imputed,
+na.rm = TRUE
+)
+)
 plot(
 d_combined,
 ylim = ylim,
 sub =
 paste(
 }
 cat("\\leavevmode\\newpage\n")
 }
 ```
-## Perform Quantile Normalization
+# Perform Quantile Normalization
+The excellent `normalize.quantiles` function from
+*[the `preprocessCore` Bioconductor package](http://bioconductor.org/packages/release/bioc/html/preprocessCore.html)*
+performs "quantile normalization" as described Bolstad *et al.* (2003),
+DOI *[10.1093/bioinformatics/19.2.185](https://doi.org/10.1093%2Fbioinformatics%2F19.2.185)*
+and *its supplementary material [http://bmbolstad.com/misc/normalize/normalize.html](http://bmbolstad.com/misc/normalize/normalize.html)*,
+i.e., it assumes that the goal is to detect
+subtle differences among grossly similar samples (having similar distributions)
+by equailzing intra-quantile quantitations.
+Unfortunately, one software library upon which it depends
+*[suffers from a concurrency defect](https://support.bioconductor.org/p/122925/#9135989)*
+that requires that a specific, non-concurrent version of the library be
+installed.  The installation command equivalent to what was used to install the library to produce the results presented here is:
+```
+conda install bioconductor-preprocesscore openblas=0.3.3
+```
 <!--
 # Apply quantile normalization using preprocessCore::normalize.quantiles
 # ---
 # tool repository: http://bioconductor.org/packages/release/bioc/html/preprocessCore.html
 #     ```
 # conda installation (necessary because of a bug in recent openblas):
 #   conda install bioconductor-preprocesscore openblas=0.3.3
 # ...
 # ---
-# normalize.quantiles {preprocessCore}	--  Quantile Normalization
+# normalize.quantiles {preprocessCore} --  Quantile Normalization
 #
 # Description:
-#   Using a normalization based upon quantiles, this function normalizes a matrix of probe level intensities.
+#   Using a normalization based upon quantiles, this function normalizes a
+#     matrix of probe level intensities.
+#
+#   THIS FUNCTIONS WILL HANDLE MISSING DATA (ie NA values), based on the
+#     assumption that the data is missing at random.
 #
 # Usage:
 #   normalize.quantiles(x, copy = TRUE, keep.names = FALSE)
 #
 # Arguments:
 #       and Variance. Bioinformatics 19(2), pp 185-193. DOI 10.1093/bioinformatics/19.2.185
 #       http://bmbolstad.com/misc/normalize/normalize.html
 # ...
 -->
 ```{r echo = FALSE, results = 'asis'}
-library(preprocessCore)
 if (nrow(quant_data_imp) > 0) {
-quant_data_imp_qn <- normalize.quantiles(as.matrix(quant_data_imp))
+quant_data_imp_qn <- preprocessCore::normalize.quantiles(
+as.matrix(quant_data_imp), keep.names = TRUE
+)
 } else {
 quant_data_imp_qn <- as.matrix(quant_data_imp)
 }
 quant_data_imp_qn <- as.data.frame(quant_data_imp_qn)
-names(quant_data_imp_qn) <- names(quant_data_imp)
 write_debug_file(quant_data_imp_qn)
 quant_data_imp_qn_log <- log10(quant_data_imp_qn)
-rownames(quant_data_imp_qn_log) <- rownames(quant_data_imp)
 write_debug_file(quant_data_imp_qn_log)
 quant_data_imp_qn_ls <- t(scale(t(log10(quant_data_imp_qn))))
-any_nan <- function(x) {
-!any(x == "NaN")
-}
 sel <- apply(quant_data_imp_qn_ls, 1, any_nan)
-quant_data_imp_qn_ls2 <- quant_data_imp_qn_ls[which(sel), ]
+quant_data_imp_qn_ls2 <- quant_data_imp_qn_ls
+quant_data_imp_qn_ls2 <- quant_data_imp_qn_ls2[which(sel), ]
 quant_data_imp_qn_ls2 <- as.data.frame(quant_data_imp_qn_ls2)
-rownames(quant_data_imp_qn_ls2) <- rownames(quant_data_imp)[sel]
 quant_data_imp_qn_ls <- as.data.frame(quant_data_imp_qn_ls)
-rownames(quant_data_imp_qn_ls) <- rownames(quant_data_imp)
 write_debug_file(quant_data_imp_qn_ls)
 write_debug_file(quant_data_imp_qn_ls2)
-#output quantile normalized data
+# Create data.frame used by ANOVA analysis
 data_table_imp_qn_lt <- cbind(full_data[1:9], quant_data_imp_qn_log)
-write.table(
-data_table_imp_qn_lt,
-file = imp_qn_lt_data_filenm,
-sep = "\t",
-col.names = TRUE,
-row.names = FALSE,
-quote = FALSE
-)
 ```
 <!-- ACE insertion begin -->
-### Checking that normalized, imputed, log-transformed sample distributions are similar:
+## Are normalized, imputed, log-transformed sample distributions similar?
 ```{r echo = FALSE, fig.dim = c(9, 5.5), results = 'asis'}
 # Save unimputed quant_data_log for plotting below
 unimputed_quant_data_log <- quant_data_log
 # log10 transform (after preparing for zero values,
 #   which should never happen...)
 quant_data_imp_qn[quant_data_imp_qn == 0] <- .000000001
 quant_data_log <- log10(quant_data_imp_qn)
-# Output with imputed, un-normalized data
-data_table_imputed <- cbind(full_data[1:9], quant_data_imp)
-write.table(
-data_table_imputed
-, file = imputed_data_filename
-, sep = "\t"
-, col.names = TRUE
-, row.names = FALSE
-, quote = FALSE
-)
 how_many_peptides <- nrow(quant_data_log)
 if ((how_many_peptides) > 0) {
 cat(
 ```
 ```{r echo = FALSE, fig.align = "left", fig.dim = c(9, 4), results = 'asis'}
 if (nrow(quant_data_log) > 1) {
 quant_data_log_stack <- stack(quant_data_log)
-ggplot(
+ggplot2::ggplot(quant_data_log_stack, ggplot2::aes(x = values)) +
-quant_data_log_stack,
+ggplot2::xlab(latex2exp::TeX("$log_{10}$(peptide intensity)")) +
-aes(x = values)
+ggplot2::ylab("Probability density") +
-) + xlab(latex2exp::TeX("$log_{10}$(peptide intensity)")) +
+ggplot2::geom_density(ggplot2::aes(group = ind, colour = ind), na.rm = TRUE)
-ylab("Probability density") +
-geom_density(
-aes(group = ind, colour = ind),
-na.rm = TRUE
-)
 } else {
 cat("No density plot because there are fewer than two peptides to plot.\n\n")
 }
 ```
 ```{r echo = FALSE, results = 'asis'}
 cat("\\leavevmode\\newpage\n")
 ```
-## Perform ANOVA Filters
+# ANOVA Analysis
 ```{r, echo = FALSE}
 # Make new data frame containing only Phosphopeptides
 #   to connect preANOVA to ANOVA (connect_df)
 connect_df <- data.frame(
 )
 colnames(connect_df) <- c("Phosphopeptide", "Intensity")
 ```
 ```{r echo = FALSE, fig.dim = c(9, 10), results = 'asis'}
-count_of_factor_levels <- length(levels(sample_factor_levels))
+count_of_treatment_levels <- length(levels(sample_treatment_levels))
-if (count_of_factor_levels < 2) {
+if (count_of_treatment_levels < 2) {
 nuke_control_sequences <-
 function(s) {
 s <- gsub("[\\]", "xyzzy_plugh", s)
 s <- gsub("[$]", "\\\\$", s)
 s <- gsub("xyzzy_plugh", "$\\\\backslash$", s)
 )
 cat("\n\n\n")
 cat("Sample group assignments are:\n",
 "\\begin{quote}\n",
-paste(regmatches(sample_name_matches, m2), collapse = "\n\n\n"),
+paste(regmatches(sample_name_matches, rx_match), collapse = "\n\n\n"),
 "\n\\end{quote}\n\n")
 } else {
 p_value_data_anova_ps <-
 apply(
 quant_data_imp_qn_log,
 1,
 anova_func,
-grouping_factor = sample_factor_levels
+grouping_factor = sample_treatment_levels,
+one_way_f = one_way_all_categories
 )
 p_value_data_anova_ps_fdr <-
 p.adjust(p_value_data_anova_ps, method = "fdr")
 p_value_data <- data.frame(
 # output ANOVA file to constructed filename,
 #   e.g.    "Outputfile_pST_ANOVA_STEP5.txt"
 #   becomes "Outpufile_pST_ANOVA_STEP5_FDR0.05.txt"
 # Re-output datasets to include p-values
+metadata_plus_p <- cbind(full_data[1:9], p_value_data[, 2:3])
 write.table(
-cbind(full_data[1:9], p_value_data[, 2:3], quant_data_imp),
+cbind(metadata_plus_p, quant_data_imp),
 file = imputed_data_filename,
 sep = "\t",
 col.names = TRUE,
 row.names = FALSE,
 quote = FALSE
 )
 write.table(
-cbind(full_data[1:9], p_value_data[, 2:3], quant_data_imp_qn_log),
+cbind(metadata_plus_p, quant_data_imp_qn_log),
 file = imp_qn_lt_data_filenm,
 sep = "\t",
 col.names = TRUE,
 row.names = FALSE,
 quote = FALSE
 # <!-- ACE insertion start -->
 if (nrow(filtered_p) && nrow(filtered_data_filtered) > 0) {
 if (first_page_suppress == 1) {
 first_page_suppress <- 0
-}
+} else {
-else {
 cat("\\newpage\n")
 }
-cat(sprintf(
+if (nrow(filtered_data_filtered) > 1) {
-"Intensities for %d peptides whose adjusted p-value < %0.2f:\n",
+subsection_header(sprintf(
-nrow(filtered_data_filtered),
+"Intensity distributions for %d phosphopeptides whose adjusted p-value < %0.2f\n",
-cutoff
+nrow(filtered_data_filtered),
-))
+cutoff
+))
+} else {
+subsection_header(sprintf(
+"Intensity distribution for one phosphopeptide (%s) whose adjusted p-value < %0.2f\n",
+rownames(filtered_data_filtered)[1],
+cutoff
+))
+}
 cat("\n\n\n")
 cat("\n\n\n")
 old_oma <- par("oma")
 old_par <- par(
 fin = c(9, 7.25)
 )
 # ref: https://r-charts.com/distribution/add-points-boxplot/
 # Vertical plot
 colnames(filtered_data_filtered) <- sample_name_matches
-boxplot(
+tryCatch(
-filtered_data_filtered,
+boxplot(
-main = "Imputed, normalized intensities", # no line plot
+filtered_data_filtered,
-las = 1,
+main = "Imputed, normalized intensities", # no line plot
-col = const_boxplot_fill,
+las = 1,
-ylab = latex2exp::TeX("$log_{10}$(peptide intensity)")
+col = const_boxplot_fill,
+ylab = latex2exp::TeX("$log_{10}$(peptide intensity)")
+),
+error = function(e) print(e)
 )
 par(old_par)
 } else {
 cat(sprintf(
 "%s < %0.2f\n\n\n\n\n",
 cutoff
 ))
 }
 if (nrow(filtered_data_filtered) > 0) {
-#Add Phosphopeptide column to anova_filtered table
+# Add Phosphopeptide column to anova_filtered table
-anova_filtered_merge <- merge(
+# The assumption here is that the first intensity is unique;
+#   this is a hokey assumption but almost definitely will
+#   be true in the real world unless there is a computation
+#   error upstream.
+anova_filtered_merge <- base::merge(
 x = connect_df,
 y = filtered_data_filtered,
 by.x = "Intensity",
 by.y = 1
 )
 anova_filtered_merge_order <- rownames(filtered_p)
+anova_filtered <- data.frame(
+ppep    = anova_filtered_merge$Phosphopeptide,
+intense = anova_filtered_merge$Intensity,
+data    = anova_filtered_merge[, 2:number_of_samples + 1]
+)
+colnames(anova_filtered) <-
+c("Phosphopeptide", colnames(filtered_data_filtered))
+# Merge qualitative columns into the ANOVA data
+output_table <- data.frame(anova_filtered$Phosphopeptide)
+output_table <- base::merge(
+x = output_table,
+y = data_table_imp_qn_lt,
+by.x = "anova_filtered.Phosphopeptide",
+by.y = "Phosphopeptide"
+)
+# Produce heatmap to visualize significance and the effect of imputation
 anova_filtered_merge_format <- sapply(
 X = filtered_p$fdr_adjusted_anova_p
 ,
 FUN = function(x) {
 else
 paste0("(%0.4e) %s")
 }
 )
-anova_filtered <- data.table(
+cat_hm_heading <- function(m, cutoff) {
-anova_filtered_merge$Phosphopeptide
+cat("\\newpage\n")
-,
+if (nrow(m) > intensity_hm_rows) {
-anova_filtered_merge$Intensity
+subsection_header(
-,
+paste(
-anova_filtered_merge[, 2:number_of_samples + 1]
+sprintf("Heatmap for the %d most-significant peptides",
-)
+intensity_hm_rows),
-colnames(anova_filtered) <-
+sprintf("whose adjusted p-value < %0.2f\n", cutoff)
-c("Phosphopeptide", colnames(filtered_data_filtered))
+)
+)
-# Merge qualitative columns into the ANOVA data
+} else {
-output_table <- data.frame(anova_filtered$Phosphopeptide)
+if (nrow(m) == 1) {
-output_table <- merge(
+return(FALSE)
-x = output_table,
+} else {
-y = data_table_imp_qn_lt,
+subsection_header(
-by.x = "anova_filtered.Phosphopeptide",
+paste(
-by.y = "Phosphopeptide"
+sprintf("Heatmap for %d usable peptides whose", nrow(m)),
-)
+sprintf("adjusted p-value < %0.2f\n", cutoff)
+)
-# Produce heatmap to visualize significance and the effect of imputation
+)
+}
+}
+cat("\n\n\n")
+cat("\n\n\n")
+return(TRUE)
+}
+# construct matrix with appropriate rownames
 m <-
 as.matrix(unimputed_quant_data_log[anova_filtered_merge_order, ])
-m_nan_rows <- rowSums(
-matrix(
-as.integer(is.na(m)),
-nrow = nrow(m)
-)
-)
-m <- m[!m_nan_rows, , drop = FALSE]
 if (nrow(m) > 0) {
 rownames_m <- rownames(m)
 rownames(m) <- sapply(
 X = seq_len(nrow(m))
 ,
 filtered_p$fdr_adjusted_anova_p[i],
 rownames_m[i]
 )
 }
 )
-how_many_peptides <- min(50, nrow(m))
+}
-number_of_peptides_found <- how_many_peptides
+# draw the heading and heatmap
-if (nrow(m) > 1) {
+if (nrow(m) > 0) {
-m_margin <- m[how_many_peptides:1, ]
+number_of_peptides_found <-
-margins <-
+draw_intensity_heatmap(
-c(max(nchar(colnames(m_margin))) * 10 / 16 # col
+m                       = m,
-, max(nchar(rownames(m_margin))) * 5 / 16 # row
+cutoff                  = cutoff,
-)
+hm_heading_function     = cat_hm_heading,
-}
+hm_main_title           = "Unimputed, unnormalized log(intensities)",
+suppress_row_dendrogram = FALSE
-cat("\\newpage\n")
-if (nrow(m) > 50) {
-cat("Heatmap for the 50 most-significant peptides",
-sprintf(
-"whose adjusted p-value < %0.2f\n",
-cutoff)
-)
-} else {
-if (nrow(m) == 1) {
-next
-} else {
-cat(
-sprintf("Heatmap for %d usable peptides whose", nrow(m)),
-sprintf("adjusted p-value < %0.2f\n", cutoff)
 )
-}
-}
-cat("\n\n\n")
-cat("\n\n\n")
-try(
-if (nrow(m) > 1) {
-old_oma <- par("oma")
-par(cex.main = 0.6)
-heatmap(
-m[how_many_peptides:1, ],
-Rowv = NA,
-Colv = NA,
-cexRow = 0.7,
-cexCol = 0.8,
-scale = "row",
-margins = margins,
-main =
-"Unimputed, unnormalized intensities",
-xlab = "",
-las = 1 #, fin = c(9, 5.5)
-)
-}
-)
 }
 }
 }
 }
 cat("\\leavevmode\n\n\n")
 ```
+```{r sqlite, echo = FALSE, fig.dim = c(9, 10), results = 'asis'}
+if (count_of_treatment_levels > 1) {
+# Prepare two-way contrasts with adjusted p-values
+# Strategy:
+# - use imputed, log-transformed data:
+#   - remember this when computing log2(fold-change)
+# - each contrast is between a combination of trt levels
+# - for each contrast, compute samples that are members
+# - compute one-way test:
+#   - use `oneway.test` (Welch test) if numbers of samples
+#     are not equivalent between trt levels
+#   - otherwise, aov is fine but offers no advantage
+# - adjust p-value, assuming that
+#   (# of pppeps)*(# of contrasts) tests were performed
+# Each contrast is between a combination of trt levels
+m2 <- combn(
+x = seq_len(length(levels(sample_treatment_levels))),
+m = 2,
+simplify = TRUE
+)
+contrast_count <- ncol(m2)
+# For each contrast, compute samples that are members
+# - local function to construct a data.frame for each contrast
+#   - the contrast in the first "column"
+f_m2 <-
+function(cntrst, lvl1, lvl2) {
+return(
+data.frame(
+contrast = cntrst,
+level = sample_treatment_levels[
+sample_treatment_levels %in%
+levels(sample_treatment_levels)[c(lvl1, lvl2)]
+],
+label = sample_name_matches[
+sample_treatment_levels %in%
+levels(sample_treatment_levels)[c(lvl1, lvl2)]
+]
+)
+)
+}
+# - compute a df for each contrast
+sample_level_dfs <- lapply(
+X = 1:contrast_count,
+FUN = function(i) f_m2(i, m2[1, i], m2[2, i])
+)
+# - compute a single df for all contrasts
+combined_contrast_df <- Reduce(f = rbind, x = sample_level_dfs)
+# - dispose objects to free resources
+rm(sample_level_dfs)
+# - write the df to a DB for later join-per-contrast
+db <- RSQLite::dbConnect(RSQLite::SQLite(), ksea_app_prep_db)
+RSQLite::dbWriteTable(
+conn = db,
+name = "contrast",
+value = combined_contrast_df,
+overwrite = TRUE
+)
+# Create UK for insert
+ddl_exec(db, "
+CREATE UNIQUE INDEX IF NOT EXISTS contrast__uk__idx
+ON contrast(contrast, label);
+"
+)
+# Create indexes for join
+ddl_exec(db, "
+-- index for join in contrast_ppep_smpl_qnlt on a.label < b.label
+CREATE INDEX IF NOT EXISTS contrast__label__idx
+ON contrast(label);
+"
+)
+ddl_exec(db, "
+-- index for joining two contrast_lvl_ppep_avg_quant on contrast
+CREATE INDEX IF NOT EXISTS contrast__contrast__idx
+ON contrast(contrast);
+"
+)
+ddl_exec(db, "
+-- index for joining two contrast_lvl_ppep_avg_quant on phophospep
+CREATE INDEX IF NOT EXISTS contrast__level__idx
+ON contrast(level);
+"
+)
+# - dispose objects to free resources
+rm(combined_contrast_df)
+# Use imputed, log-transformed data
+# - remember that this was donoe when computing log2(fold-change)
+# - melt data matrix for use in later join-per-contrast
+casted <- cbind(
+data.frame(vrbl = rownames(quant_data_imp_qn_log)),
+quant_data_imp_qn_log
+)
+quant_data_imp_qn_log_melted <- reshape2::melt(
+casted,
+id.vars = "vrbl"
+)
+colnames(quant_data_imp_qn_log_melted) <-
+c("phosphopep", "sample", "quant")
+# - dispose objects to free resources
+rm(casted)
+# - write the df to a DB for use in later join-per-contrast
+RSQLite::dbWriteTable(
+conn = db,
+name = "ppep_smpl_qnlt",
+value = quant_data_imp_qn_log_melted,
+overwrite = TRUE
+)
+# Create UK for insert
+ddl_exec(db, "
+CREATE UNIQUE INDEX IF NOT EXISTS ppep_smpl_qnlt__uk__idx
+ON ppep_smpl_qnlt(phosphopep, sample);
+"
+)
+# Create index for join
+ddl_exec(db, "
+-- index for join in contrast_ppep_smpl_qnlt
+CREATE INDEX IF NOT EXISTS ppep_smpl_qnlt__sample__idx
+ON ppep_smpl_qnlt(sample);
+"
+)
+ddl_exec(db, "
+-- index for joining two contrast_lvl_ppep_avg_quant on phopho.pep
+CREATE INDEX IF NOT EXISTS ppep_smpl_qnlt__phosphopep__idx
+ON ppep_smpl_qnlt(phosphopep);
+"
+)
+# - dispose objects to free resources
+rm(quant_data_imp_qn_log_melted)
+# - drop views if exist
+ddl_exec(db, "
+-- drop view dependent on contrast_lvl_ppep_avg_quant
+DROP VIEW IF EXISTS v_contrast_log2_fc;
+"
+)
+ddl_exec(db, "
+-- drop table dependent on contrast_ppep_smpl_qnlt
+DROP TABLE IF EXISTS contrast_lvl_ppep_avg_quant;
+"
+)
+ddl_exec(db, "
+DROP TABLE IF EXISTS contrast_lvl_lvl_metadata;
+"
+)
+ddl_exec(db, "
+DROP VIEW IF EXISTS v_contrast_lvl_metadata;
+"
+)
+ddl_exec(db, "
+-- drop view dependent on contrast_ppep_smpl_qnlt
+DROP VIEW IF EXISTS v_contrast_lvl_ppep_avg_quant;
+"
+)
+ddl_exec(db, "
+DROP VIEW IF EXISTS v_contrast_lvl_lvl;
+"
+)
+ddl_exec(db, "
+-- drop view upon which other views depend
+DROP VIEW IF EXISTS contrast_ppep_smpl_qnlt;
+"
+)
+# - create view
+dml_no_rows_exec(db, "
+-- view contrast_ppep_smpl_qnlt is used for each phopshopep to
+--   compute p-value for test of trt effect for two trt levels
+CREATE VIEW contrast_ppep_smpl_qnlt
+AS
+SELECT  contrast,
+level,
+phosphopep,
+sample,
+quant
+FROM  contrast AS c,
+ppep_smpl_qnlt AS q
+WHERE q.sample = c.label
+ORDER BY contrast, level, phosphopep
+;
+"
+)
+# - create simplification views
+dml_no_rows_exec(db, "
+CREATE VIEW v_contrast_lvl_metadata
+AS
+SELECT  contrast,
+level,
+group_concat(label, ';') AS samples
+FROM  contrast
+GROUP BY contrast, level
+/* view v_contrast_lvl_metadata is used
+to simplify creation of table contrast_lvl_lvl_metadata */
+;
+"
+)
+dml_no_rows_exec(db, "
+CREATE VIEW v_contrast_lvl_ppep_avg_quant
+AS
+SELECT  contrast,
+level,
+phosphopep,
+avg(quant)                AS avg_quant
+FROM  contrast_ppep_smpl_qnlt
+GROUP BY contrast, level, phosphopep
+/* view v_contrast_lvl_ppep_avg_quant is used
+to simplify view v_contrast_log2_fc */
+;
+"
+)
+# - create contrast-metadata table
+dml_no_rows_exec(db, "
+CREATE TABLE contrast_lvl_lvl_metadata
+AS
+SELECT DISTINCT
+a.contrast              AS ab_contrast,
+a.level                 AS a_level,
+b.level                 AS b_level,
+a.samples               AS a_samples,
+b.samples               AS b_samples,
+'log2(level_'||a.level||'/level_'||b.level||')'
+AS fc_description
+FROM  v_contrast_lvl_metadata AS a,
+v_contrast_lvl_metadata AS b
+WHERE a.contrast = b.contrast
+AND a.level > b.level
+/* view v_contrast_lvl_lvl is used
+to simplify view v_contrast_log2_fc */
+;
+"
+)
+# - create pseudo-materialized view table
+dml_no_rows_exec(db, "
+CREATE VIEW v_contrast_lvl_lvl
+AS
+SELECT DISTINCT
+a.contrast AS ab_contrast,
+a.level    AS a_level,
+b.level    AS b_level
+FROM  contrast   AS a,
+contrast   AS b
+WHERE a.contrast = b.contrast
+AND a.level > b.level
+/* view v_contrast_lvl_lvl is used
+to simplify view v_contrast_log2_fc */
+;
+"
+)
+# - create view to compute log2(fold-change)
+dml_no_rows_exec(db, "
+CREATE VIEW v_contrast_log2_fc
+AS
+SELECT  ab.ab_contrast                         AS contrast,
+m.a_level                              AS a_level,
+c.avg_quant                            AS a_quant,
+m.a_samples                            AS a_samples,
+ab.b_level                             AS b_level,
+d.avg_quant                            AS b_quant,
+m.b_samples                            AS b_samples,
+m.fc_description                       AS fc_description,
+3.32193 * ( d.avg_quant - c.avg_quant) AS log2_fc,
+d.phosphopep                           AS phosphopep
+FROM  contrast_lvl_lvl_metadata                  AS m,
+v_contrast_lvl_ppep_avg_quant              AS d,
+v_contrast_lvl_lvl AS ab
+INNER JOIN v_contrast_lvl_ppep_avg_quant AS c
+ON c.contrast = ab.ab_contrast
+AND c.level   = ab.a_level
+WHERE d.contrast    = ab.ab_contrast
+AND m.ab_contrast = ab.ab_contrast
+AND d.level       = ab.b_level
+AND d.phosphopep  = c.phosphopep
+/* view to compute log2(fold-change) */
+;
+"
+)
+# For each contrast, compute samples that are members
+# compute one-way test:
+# - use `oneway.test` (Welch test) if numbers of samples
+#   are not equivalent between trt levels
+# - otherwise, aov is fine but offers no advantage
+for (contrast in contrast_count:2) {
+invisible(contrast)
+}
+for (contrast in 1:contrast_count) {
+contrast_df <- sqldf::sqldf(
+x = paste0("
+SELECT level, phosphopep, sample, quant
+FROM contrast_ppep_smpl_qnlt
+WHERE contrast = ", contrast, "
+ORDER BY phosphopep, level, sample
+"),
+connection = db
+)
+contrast_cast <- reshape2::dcast(
+data = contrast_df,
+formula = phosphopep ~ sample,
+value.var = "quant"
+)
+contrast_cast_ncol <- ncol(contrast_cast)
+contrast_cast_data <- contrast_cast[, 2:contrast_cast_ncol]
+contrast_cast_samples <- colnames(contrast_cast_data)
+# - order grouping_factor by order of sample columns of contrast_cast_data
+grouping_factor <- sqldf::sqldf(
+x = paste0("
+SELECT sample, level
+FROM contrast_ppep_smpl_qnlt
+WHERE contrast = ", contrast, "
+ORDER BY phosphopep, level, sample
+LIMIT ", contrast_cast_ncol - 1
+),
+connection = db
+)
+rownames(grouping_factor) <- grouping_factor$sample
+grouping_factor <- grouping_factor[, "level", drop = FALSE]
+# - run the two-level (one-way) test
+p_value_data_contrast_ps <-
+apply(
+X = contrast_cast_data,
+MARGIN = 1, # apply to rows
+FUN = anova_func,
+grouping_factor =
+as.factor(as.numeric(grouping_factor$level)), # anova_func arg2
+one_way_f = one_way_two_categories, # anova_func arg3
+simplify = TRUE # TRUE is the default for simplify
+)
+contrast_data_adj_p_values <- p.adjust(
+p = p_value_data_contrast_ps,
+method = "fdr",
+n = length(p_value_data_contrast_ps) # this is the default, length(p)
+)
+# - compute the fold-change
+contrast_p_df <-
+data.frame(
+contrast = contrast,
+phosphopep = contrast_cast$phosphopep,
+p_value_raw = p_value_data_contrast_ps,
+p_value_adj = contrast_data_adj_p_values
+)
+db_write_table_overwrite <- (contrast < 2)
+db_write_table_append <- !db_write_table_overwrite
+RSQLite::dbWriteTable(
+conn = db,
+name = "contrast_ppep_p_val",
+value = contrast_p_df,
+append = db_write_table_append
+)
+# Create UK for insert
+ddl_exec(db, "
+CREATE UNIQUE INDEX IF NOT EXISTS contrast_ppep_p_val__uk__idx
+ON contrast_ppep_p_val(phosphopep, contrast);
+"
+)
+}
+# Perhaps this could be done more elegantly using unique keys
+#   or creating the tables before saving data to them, but this
+#   is fast and, if the database exists on disk rather than in
+#   memory, it doesn't stress memory.
+dml_no_rows_exec(db, "
+CREATE TEMP table contrast_log2_fc
+AS
+SELECT *
+FROM  v_contrast_log2_fc
+ORDER BY contrast, phosphopep
+;
+"
+)
+dml_no_rows_exec(db, "
+CREATE TEMP table ppep_p_val
+AS
+SELECT  p_value_raw,
+p_value_adj,
+contrast   AS p_val_contrast,
+phosphopep AS p_val_ppep
+FROM  contrast_ppep_p_val
+ORDER BY contrast, phosphopep
+;
+"
+)
+dml_no_rows_exec(db, "
+DROP TABLE IF EXISTS contrast_log2_fc_p_val
+;
+"
+)
+dml_no_rows_exec(db, "
+CREATE TABLE contrast_log2_fc_p_val
+AS
+SELECT  a.*,
+b.p_value_raw,
+b.p_value_adj,
+b.p_val_contrast,
+b.p_val_ppep
+FROM  contrast_log2_fc a, ppep_p_val b
+WHERE a.rowid = b.rowid
+AND a.phosphopep = b.p_val_ppep
+;
+"
+)
+# Create UK
+ddl_exec(db, "
+CREATE UNIQUE INDEX IF NOT EXISTS contrast_log2_fc_p_val__uk__idx
+ON contrast_log2_fc_p_val(phosphopep, contrast);
+"
+)
+# Create indices for future queries
+ddl_exec(db, "
+CREATE INDEX IF NOT EXISTS contrast_log2_fc_p_val__contrast__idx
+ON contrast_log2_fc_p_val(contrast);
+"
+)
+ddl_exec(db, "
+CREATE INDEX IF NOT EXISTS contrast_log2_fc_p_val__phosphopep__idx
+ON contrast_log2_fc_p_val(phosphopep);
+"
+)
+ddl_exec(db, "
+CREATE INDEX IF NOT EXISTS contrast_log2_fc_p_val__p_value_raw__idx
+ON contrast_log2_fc_p_val(p_value_raw);
+"
+)
+ddl_exec(db, "
+CREATE INDEX IF NOT EXISTS contrast_log2_fc_p_val__p_value_adj__idx
+ON contrast_log2_fc_p_val(p_value_adj);
+"
+)
+dml_no_rows_exec(db, "
+DROP VIEW IF EXISTS v_contrast_log2_fc_p_val
+;
+"
+)
+dml_no_rows_exec(db, "
+CREATE VIEW v_contrast_log2_fc_p_val
+AS
+SELECT  contrast,
+a_level,
+a_samples,
+b_level,
+b_samples,
+a_quant,
+b_quant,
+fc_description,
+log2_fc,
+p_value_raw,
+p_value_adj,
+phosphopep
+FROM  contrast_log2_fc_p_val
+ORDER BY contrast, phosphopep
+;
+"
+)
+ddl_exec(db, "
+DROP TABLE IF EXISTS kseaapp_metadata
+;
+"
+)
+dml_no_rows_exec(db, "
+CREATE TABLE kseaapp_metadata
+AS
+WITH extended(deppep, ppep, gene_name, uniprot_id, phosphoresidue) AS (
+SELECT DISTINCT
+deppep.seq,
+ppep.seq,
+GeneName||';',
+UniProtID||';',
+PhosphoResidue||';'
+FROM
+ppep, deppep, mrgfltr_metadata
+WHERE
+mrgfltr_metadata.ppep_id = ppep.id
+AND
+ppep.deppep_id = deppep.id
+)
+SELECT
+ppep                                                  AS `ppep`,
+SUBSTR(uniprot_id,     1, INSTR(uniprot_id,';') - 1 ) AS `Protein`,
+SUBSTR(gene_name,      1, INSTR(gene_name,';')  - 1 ) AS `Gene`,
+deppep                                                AS `Peptide`,
+REPLACE(
+REPLACE(
+SUBSTR(phosphoresidue, 1, INSTR(phosphoresidue,';') - 1 ),
+'p',
+''
+),
+', ',
+';'
+)                                                   AS `Residue.Both`
+FROM extended
+;
+"
+)
+# Create indexes for join
+ddl_exec(db, "
+CREATE INDEX IF NOT EXISTS kseaapp_metadata__ppep__idx
+ON kseaapp_metadata(ppep);
+"
+)
+ddl_exec(db, "
+DROP VIEW IF EXISTS v_kseaapp_contrast
+;
+"
+)
+dml_no_rows_exec(db, "
+CREATE VIEW v_kseaapp_contrast
+AS
+SELECT  a.*, b.Protein, b.Gene, b.Peptide, b.`Residue.Both`
+FROM  v_contrast_log2_fc_p_val a, kseaapp_metadata b
+WHERE b.ppep = a.phosphopep
+;
+"
+)
+ddl_exec(db, "
+DROP VIEW IF EXISTS v_kseaapp_input
+;
+"
+)
+dml_no_rows_exec(db, "
+CREATE VIEW v_kseaapp_input
+AS
+SELECT  v.contrast,
+v.phosphopep,
+m.`Protein`,
+m.`Gene`,
+m.`Peptide`,
+m.`Residue.Both`,
+v.p_value_raw AS `p`,
+v.log2_fc AS `FC`
+FROM  kseaapp_metadata AS m,
+v_contrast_log2_fc_p_val AS v
+WHERE m.ppep = v.phosphopep
+AND NOT m.`Gene` = 'No_Gene_Name'
+AND NOT v.log2_fc = 0
+;
+"
+)
+}
+```
+```{r echo = FALSE, results = 'asis'}
+cat("\\newpage\n")
+```
+# KSEA Analysis
+Results of Kinase-Substrate Enrichment Analysis are presented here, if the substrates for any kinases are relatively enriched.   Enrichments are found by the CRAN `KSEAapp` package:
+- The package is available on CRAN, at https:/cran.r-project.org/package=KSEAapp
+- The method used is described in Casado et al. (2013) [doi:10.1126/scisignal.2003573](https:/doi.org/10.1126/scisignal.2003573) and Wiredja et al (2017) [doi:10.1093/bioinformatics/btx415](https:/doi.org/10.1093/bioinformatics/btx415).
+- An online alternative (supporting only analysis of human data) is available at [https:/casecpb.shinyapps.io/ksea/](https:/casecpb.shinyapps.io/ksea/).
+For each kinase, $i$, and each two-way contrast of treatments, $j$, an enrichment $z$-score is computed as:
+$$
+\text{kinase enrichment score}_{j,i} = \frac{(\overline{s}_{j,i} - \overline{p}_j)\sqrt{m_{j,i}}}{\delta_j}
+$$
+and fold-enrichment is computed as:
+$$
+\text{Enrichment}_{j,i} = \frac{\overline{s}_{j,i}}{\overline{p}_j}
+$$
+where:
+- $\overline{s}_{j,i}$ is the mean $\log_2 (|\text{fold-change|})$ in intensities (for contrast $j$) of known substrates of the kinase $i$,
+- $\overline{p}_j$ is the mean $\log_2 (|\text{fold-change}|)$ of all phosphosites identified in contrast $j$, and
+- $m_{j,i}$ is the total number of phosphosite substrates of kinase $i$ identified in contrast $j$,
+- $\delta_j$ is the standard deviation of the $\log_2 (|\text{fold-change}|)$ for contrast $j$ across all phosphosites in the dataset.
+- Note that the absolute value of fold-change is used so that both increased and decreased substrates of a kinase will contribute to its enrichment score.
+$\text{FDR}_{j,i}$ is computed from the $p$-value for the z-score using the R `stats::p.adjust` function, applying the False Discovery Rate correction from Benjamini and Hochberg (1995) [doi:10.1111/j.2517-6161.1995.tb02031.x](https:/doi.org/10.1111/j.2517-6161.1995.tb02031.x)
+Color intensity in heatmaps reflects magnitude of $z$-score for enrichment of respective kinase in respective contrast; hue reflects the sign of the $z$-score (blue, negative; red, positive).
+Asterisks in heatmaps reflect enrichments that are significant at `r ksea_cutoff_statistic` < `r ksea_cutoff_threshold`.
+- Kinase names are generally as presented at Phospho.ELM [http://phospho.elm.eu.org/kinases.html](http://phospho.elm.eu.org/kinases.html) (when available), although Phospho.ELM data are not yet incorporated into this analysis.
+- Kinase names having the suffix '(HPRD)' are as presented at [http://hprd.org/serine_motifs](http://hprd.org/serine_motifs) and [http://hprd.org/tyrosine_motifs](http://hprd.org/tyrosine_motifs) and are as originally reported in the Amanchy et al., 2007 (doi: [10.1038/nbt0307-285](https://doi.org/10.1038/nbt0307-285)).
+- Kinase-strate deata were also taken from [http://networkin.science/download.shtml](http://networkin.science/download.shtml) and from PhosphoSitePlus [https://www.phosphosite.org/staticDownloads](https://www.phosphosite.org/staticDownloads).
+```{r ksea, echo = FALSE, fig.dim = c(9, 10), results = 'asis'}
+db <- RSQLite::dbConnect(RSQLite::SQLite(), ksea_app_prep_db)
+# -- eliminate the table that's about to be defined
+ddl_exec(db, "
+DROP TABLE IF EXISTS site_metadata;
+")
+# -- define the site_metadata table
+ddl_exec(db, "
+CREATE TABLE site_metadata(
+id            INTEGER PRIMARY KEY
+, site_type_id  INTEGER REFERENCES site_type(id)
+, full          TEXT    UNIQUE ON CONFLICT IGNORE
+, abbrev        TEXT
+, pattern       TEXT
+, motif         TEXT
+);
+")
+# -- populate the table with initial values
+ddl_exec(db, "
+INSERT INTO site_metadata(full, abbrev, site_type_id)
+SELECT  DISTINCT kinase_map, kinase_map, site_type_id
+FROM  ppep_gene_site
+ORDER BY kinase_map;
+")
+# -- drop bogus KSData view if exists
+ddl_exec(db, "
+DROP VIEW IF EXISTS ks_data_v;
+")
+# -- create view to serve as an impostor for  KSEAapp::KSData
+ddl_exec(db, "
+CREATE VIEW IF NOT EXISTS ks_data_v
+AS
+SELECT
+'NA' AS KINASE,
+'NA' AS KIN_ACC_ID,
+kinase_map AS GENE,
+'NA' AS KIN_ORGANISM,
+'NA' AS SUBSTRATE,
+0 AS SUB_GENE_ID,
+'NA' AS SUB_ACC_ID,
+gene_names AS SUB_GENE,
+'NA' AS SUB_ORGANISM,
+phospho_peptide AS SUB_MOD_RSD,
+0 AS SITE_GROUP_ID,
+'NA' AS 'SITE_7AA',
+2 AS networkin_score,
+type_name AS Source
+FROM ppep_gene_site_view;
+")
+contrast_metadata_df <-
+sqldf::sqldf("select * from contrast_lvl_lvl_metadata", connection = db)
+rslt <- new_env()
+rslt$score_list <- list()
+rslt$name_list  <- list()
+rslt$longname_list  <- list()
+ddl_exec(db, "
+DROP TABLE IF EXISTS contrast_ksea_scores;
+"
+)
+next_index <- 0
+err_na_subscr_df_const <-
+"missing values are not allowed in subscripted assignments of data frames"
+for (i_cntrst in seq_len(nrow(contrast_metadata_df))) {
+cntrst_a_level <- contrast_metadata_df[i_cntrst, "a_level"]
+cntrst_b_level <- contrast_metadata_df[i_cntrst, "b_level"]
+cntrst_fold_change <- contrast_metadata_df[i_cntrst, 6]
+contrast_label <- sprintf("%s -> %s", cntrst_b_level, cntrst_a_level)
+contrast_longlabel <- (
+sprintf(
+"Trt %s {%s} -> Trt %s {%s}",
+contrast_metadata_df[i_cntrst, "b_level"],
+gsub(
+pattern = ";",
+replacement = ", ",
+x = contrast_metadata_df[i_cntrst, "b_samples"],
+fixed = TRUE
+),
+contrast_metadata_df[i_cntrst, "a_level"],
+gsub(
+pattern = ";",
+replacement = ", ",
+x = contrast_metadata_df[i_cntrst, "a_samples"],
+fixed = TRUE
+)
+)
+)
+kseaapp_input <-
+sqldf::sqldf(
+x = sprintf("
+SELECT `Protein`, `Gene`, `Peptide`, phosphopep AS `Residue.Both`, `p`, `FC`
+FROM v_kseaapp_input
+WHERE contrast = %d
+",
+i_cntrst
+),
+connection = db
+)
+pseudo_ksdata <- dbReadTable(db, "ks_data_v")
+# This hack is because SQL table has the log2-transformed values
+kseaapp_input[, "FC"] <- 2 ** kseaapp_input[, "FC", drop = TRUE]
+main_title <- (
+sprintf(
+"Change from treatment %s to treatment %s",
+contrast_metadata_df[i_cntrst, "b_level"],
+contrast_metadata_df[i_cntrst, "a_level"]
+)
+)
+sub_title <- contrast_longlabel
+tryCatch(
+expr = {
+ksea_scores_rslt <-
+ksea_scores(
+ksdata           = pseudo_ksdata, # KSEAapp::KSData,
+px               = kseaapp_input,
+networkin        = TRUE,
+networkin_cutoff = 2
+)
+if (0 < sum(!is.nan(ksea_scores_rslt$FDR))) {
+next_index <- 1 + next_index
+rslt$score_list[[next_index]] <- ksea_scores_rslt
+rslt$name_list[[next_index]] <- contrast_label
+rslt$longname_list[[next_index]] <- contrast_longlabel
+low_fdr_print(
+rslt = rslt,
+i_cntrst = i_cntrst,
+i = next_index,
+a_level = cntrst_a_level,
+b_level = cntrst_b_level,
+fold_change = cntrst_fold_change,
+caption = contrast_longlabel
+)
+}
+},
+error = function(e) str(e)
+)
+}
+plotted_kinases <- NULL
+if (length(rslt$score_list) > 1) {
+for (i in seq_len(length(ksea_heatmap_titles))) {
+hdr <- ksea_heatmap_titles[[i]]
+which_kinases <- i
+cat("\\clearpage\n\\begin{center}\n")
+if (i == const_ksea_astrsk_kinases) {
+subsection_header(hdr)
+} else {
+subsection_header(hdr)
+}
+cat("\\end{center}\n")
+plotted_kinases <- ksea_heatmap(
+# the data frame outputs from the KSEA.Scores() function, in list format
+score_list = rslt$score_list,
+# a character vector of all the sample names for heatmap annotation:
+# - the names must be in the same order as the data in score_list
+# - please avoid long names, as they may get cropped in the final image
+sample_labels = rslt$name_list,
+# character string of either "p.value" or "FDR" indicating the data column
+#   to use for marking statistically significant scores
+stats = c("p.value", "FDR")[2],
+# a numeric value between 0 and infinity indicating the min. number of
+#   substrates a kinase must have to be included in the heatmap
+m_cutoff = 1,
+# a numeric value between 0 and 1 indicating the p-value/FDR cutoff
+#   for indicating significant kinases in the heatmap
+p_cutoff = 0.05,
+# a binary input of TRUE or FALSE, indicating whether or not to perform
+#   hierarchical clustering of the sample columns
+sample_cluster = TRUE,
+# a binary input of TRUE or FALSE, indicating whether or not to export
+#   the heatmap as a .png image into the working directory
+export = FALSE,
+# additional arguments to gplots::heatmap.2, such as:
+# - main: main title of plot
+# - xlab: x-axis label
+# - ylab: y-axis label
+xlab = "Contrast",
+ylab = "Kinase",
+# print which kinases:
+# - 1 : all kinases
+# - 2 : significant kinases
+# - 3 : non-significant kinases
+which_kinases = which_kinases
+)
+cat("\\begin{center}\n")
+cat("Color intensities reflects $z$-score magnitudes; hue reflects $z$-score sign.  Asterisks reflect significance.\n")
+cat("\\end{center}\n")
+} # end for (i in ...
+} # end if (length ...
+for (i_cntrst in seq_len(length(rslt$score_list))) {
+next_index <- i_cntrst
+cntrst_a_level <- contrast_metadata_df[i_cntrst, "a_level"]
+cntrst_b_level <- contrast_metadata_df[i_cntrst, "b_level"]
+cntrst_fold_change <- contrast_metadata_df[i_cntrst, 6]
+contrast_label <- sprintf("%s -> %s", cntrst_b_level, cntrst_a_level)
+contrast_longlabel <- (
+sprintf(
+"Trt %s {%s} -> Trt %s {%s}",
+contrast_metadata_df[i_cntrst, "b_level"],
+gsub(
+pattern = ";",
+replacement = ", ",
+x = contrast_metadata_df[i_cntrst, "b_samples"],
+fixed = TRUE
+),
+contrast_metadata_df[i_cntrst, "a_level"],
+gsub(
+pattern = ";",
+replacement = ", ",
+x = contrast_metadata_df[i_cntrst, "a_samples"],
+fixed = TRUE
+)
+)
+)
+main_title <- (
+sprintf(
+"Change from treatment %s to treatment %s",
+contrast_metadata_df[i_cntrst, "b_level"],
+contrast_metadata_df[i_cntrst, "a_level"]
+)
+)
+sub_title <- contrast_longlabel
+tryCatch(
+expr = {
+ksea_scores_rslt <- rslt$score_list[[next_index]]
+if (0 < sum(!is.nan(ksea_scores_rslt$FDR))) {
+low_fdr_barplot(
+rslt = rslt,
+i_cntrst = i_cntrst,
+i = next_index,
+a_level = cntrst_a_level,
+b_level = cntrst_b_level,
+fold_change = cntrst_fold_change,
+caption = contrast_longlabel
+)
+}
+},
+error = function(e) str(e)
+)
+}
+```
+```{r enriched, echo = FALSE, fig.dim = c(9, 10), results = 'asis'}
+# Use enriched kinases to find enriched kinase-substrate pairs
+enriched_kinases <- data.frame(kinase = ls(ksea_asterisk_hash))
+all_enriched_substrates <- sqldf("
+SELECT
+gene AS kinase,
+ppep,
+'('||group_concat(gene||'-'||sub_gene)||') '||ppep AS label
+FROM (
+SELECT DISTINCT gene, sub_gene, SUB_MOD_RSD AS ppep
+FROM pseudo_ksdata
+WHERE GENE IN (SELECT kinase FROM enriched_kinases)
+)
+GROUP BY ppep
+")
+# helper used to label per-kinase substrate enrichment figure
+cat_enriched_heading <- function(m, cut_args) {
+cutoff <- cut_args$cutoff
+kinase <- cut_args$kinase
+statistic <- cut_args$statistic
+threshold <- cut_args$threshold
+cat("\\newpage\n")
+if (nrow(m) > intensity_hm_rows) {
+subsection_header(
+paste(
+sprintf(
+"Lowest p-valued %d (of %d) enriched %s-substrates,",
+intensity_hm_rows,
+nrow(m),
+kinase
+),
+sprintf(" KSEA %s < %0.2f\n", statistic, threshold)
+)
+)
+} else {
+if (nrow(m) == 1) {
+return(FALSE)
+} else {
+subsection_header(
+paste(
+sprintf(
+"%d enriched %s-substrates,",
+nrow(m),
+kinase
+),
+sprintf(
+" KSEA %s < %0.2f\n",
+statistic,
+threshold
+)
+)
+)
+}
+}
+cat("\n\n\n")
+cat("\n\n\n")
+return(TRUE)
+}
+# Disabling heatmaps for substrates pending decision whether to eliminate them altogether
+if (FALSE)
+for (kinase_name in sort(enriched_kinases$kinase)) {
+enriched_substrates <-
+all_enriched_substrates[
+all_enriched_substrates$kinase == kinase_name,
+,
+drop = FALSE
+]
+# Get the intensity values for the heatmap
+enriched_intensities <-
+as.matrix(unimputed_quant_data_log[enriched_substrates$ppep, , drop = FALSE])
+# Remove rows having too many NA values to be relevant
+na_counter <- is.na(enriched_intensities)
+na_counts <- apply(na_counter, 1, sum)
+enriched_intensities <-
+enriched_intensities[na_counts < ncol(enriched_intensities) / 2, , drop = FALSE]
+# Rename the rows with the display-name for the heatmap
+rownames(enriched_intensities) <-
+sapply(
+X = rownames(enriched_intensities),
+FUN = function(rn) {
+enriched_substrates[enriched_substrates$ppep == rn, "label"]
+}
+)
+# Format as matrix for heatmap
+m <- as.matrix(enriched_intensities)
+# Draw the heading and heatmap
+if (nrow(m) > 0) {
+cut_args <- new_env()
+cut_args$cutoff <- cutoff
+cut_args$kinase <- kinase_name
+cut_args$statistic <- ksea_cutoff_statistic
+cut_args$threshold <- ksea_cutoff_threshold
+number_of_peptides_found <-
+draw_intensity_heatmap(
+m                       = m,
+cutoff                  = cut_args,
+hm_heading_function     = cat_enriched_heading,
+hm_main_title
+= "Unnormalized (zero-imputed) intensities of enriched kinase-substrates",
+suppress_row_dendrogram = FALSE
+)
+}
+}
+# Write output tabular files
+# get kinase, ppep, concat(kinase) tuples for enriched kinases
+kinase_ppep_label <- sqldf("
+WITH
+t(ppep, label) AS
+(
+SELECT DISTINCT
+SUB_MOD_RSD AS ppep,
+group_concat(gene, '; ') AS label
+FROM pseudo_ksdata
+WHERE GENE IN (SELECT kinase FROM enriched_kinases)
+GROUP BY ppep
+),
+k(kinase, ppep_join) AS
+(
+SELECT DISTINCT gene AS kinase, SUB_MOD_RSD AS ppep_join
+FROM pseudo_ksdata
+WHERE GENE IN (SELECT kinase FROM enriched_kinases)
+)
+SELECT k.kinase, t.ppep, t.label
+FROM  t,  k
+WHERE t.ppep = k.ppep_join
+ORDER BY k.kinase, t.ppep
+")
+# extract what we need from full_data
+impish <- cbind(rownames(quant_data_imp), quant_data_imp)
+colnames(impish)[1] <- "Phosphopeptide"
+data_table_imputed_sql <- "
+SELECT
+f.*,
+k.label AS KSEA_enrichments,
+q.*
+FROM
+metadata_plus_p f
+LEFT JOIN kinase_ppep_label k
+ON f.Phosphopeptide = k.ppep,
+impish q
+WHERE
+f.Phosphopeptide = q.Phosphopeptide
+"
+data_table_imputed <- sqldf(data_table_imputed_sql)
+# Zap the duplicated 'Phosphopeptide' column named 'ppep'
+data_table_imputed <-
+data_table_imputed[, c(1:12, 14:ncol(data_table_imputed))]
+# Output with imputed, un-normalized data
+write.table(
+data_table_imputed
+, file = imputed_data_filename
+, sep = "\t"
+, col.names = TRUE
+, row.names = FALSE
+, quote = FALSE
+)
+#output quantile normalized data
+impish <- cbind(rownames(quant_data_imp_qn_log), quant_data_imp_qn_log)
+colnames(impish)[1] <- "Phosphopeptide"
+data_table_imputed <- sqldf(data_table_imputed_sql)
+# Zap the duplicated 'Phosphopeptide' column named 'ppep'
+data_table_imputed <-
+data_table_imputed[, c(1:12, 14:ncol(data_table_imputed))]
+write.table(
+data_table_imputed,
+file = imp_qn_lt_data_filenm,
+sep = "\t",
+col.names = TRUE,
+row.names = FALSE,
+quote = FALSE
+)
+ppep_kinase <- sqldf("
+SELECT DISTINCT k.ppep, k.kinase
+FROM (
+SELECT DISTINCT gene AS kinase, SUB_MOD_RSD AS ppep
+FROM pseudo_ksdata
+WHERE GENE IN (SELECT kinase FROM enriched_kinases)
+) k
+ORDER BY k.ppep, k.kinase
+")
+RSQLite::dbWriteTable(
+conn = db,
+name = "ksea_enriched_ks",
+value = ppep_kinase,
+append = FALSE
+)
+RSQLite::dbWriteTable(
+conn = db,
+name = "anova_signif",
+value = p_value_data,
+append = FALSE
+)
+ddl_exec(db, "
+DROP VIEW IF EXISTS stats_metadata_v;
+"
+)
+dml_no_rows_exec(db, "
+CREATE VIEW stats_metadata_v
+AS
+SELECT DISTINCT  m.*,
+p.raw_anova_p,
+p.fdr_adjusted_anova_p,
+kek.kinase AS ksea_enrichments
+FROM
+mrgfltr_metadata_view m
+LEFT JOIN anova_signif p
+ON m.phospho_peptide = p.phosphopeptide
+LEFT JOIN ksea_enriched_ks kek
+ON m.phospho_peptide = kek.ppep
+;
+"
+)
+write.table(
+dbReadTable(db, "stats_metadata_v"),
+file = anova_ksea_mtdt_file,
+sep = "\t",
+col.names = TRUE,
+row.names = FALSE,
+quote = FALSE
+)
+```
+```{r parmlist, echo = FALSE, fig.dim = c(9, 10), results = 'asis'}
+cat("\\leavevmode\n\n\n")
+# write parameters to report
+param_unlist <- unlist(as.list(params))
+param_df <- data.frame(
+parameter = paste0("\\verb@", names(param_unlist), "@"),
+value = paste0("\\verb@", gsub("$", "\\$", param_unlist, fixed = TRUE), "@")
+)
+data_frame_latex(
+x = param_df,
+justification = "p{0.35\\linewidth} p{0.6\\linewidth}",
+centered = TRUE,
+caption = "Input parameters",
+anchor = const_table_anchor_bp,
+underscore_whack = FALSE
+)
+# write parameters to SQLite output
+mqppep_anova_script_param_df <- data.frame(
+script    = "mqppep_anova_script.Rmd",
+parameter = names(param_unlist),
+value     = param_unlist
+)
+ddl_exec(db, "
+DROP TABLE IF EXISTS script_parameter;
+"
+)
+ddl_exec(db, "
+CREATE TABLE IF NOT EXISTS script_parameter(
+script    TEXT,
+parameter TEXT,
+value     ANY,
+UNIQUE (script, parameter) ON CONFLICT REPLACE
+)
+;
+"
+)
+RSQLite::dbWriteTable(
+conn = db,
+name = "script_parameter",
+value = mqppep_anova_script_param_df,
+append = TRUE
+)
+# We are done with output
+RSQLite::dbDisconnect(db)
+```
+<!--
+There's gotta be a better way...
+loaded_packages_df <-  sessioninfo::package_info("loaded")
+loaded_packages_df[, "library"] <- as.character(loaded_packages_df$library)
+loaded_packages_df <- data.frame(
+package = loaded_packages_df$package,
+version = loaded_packages_df$loadedversion,
+date    = loaded_packages_df$date
+)
+data_frame_latex(
+x = loaded_packages_df,
+justification = "l | l l",
+centered = FALSE,
+caption = "Loaded R packages",
+anchor = const_table_anchor_bp
+)
+-->

Mercurial > repos > eschen42 > mqppep_anova

comparison mqppep_anova_script.Rmd @ 22:61adb8801b73 draft