normalization and protOrigin

DESCRIPTION

@@ -26,13 +26,17 @@ Imports:
     Rsolnp,
     stringr,
     SummarizedExperiment
-Suggests: 
+Suggests:
     covr,
     rmarkdown,
     testthat (>= 3.0.0),
-    vsn, 
+    vsn,
     tidySummarizedExperiment,
-    seqinr
+    seqinr,
+    BiocStyle,
+    knitr,
+    pak
 Config/testthat/edition: 3
+VignetteBuilder: knitr
 Depends: 
     igraph (>= 2.0.0)

R/02_peptidePreprocessing.R

@@ -42,7 +42,7 @@
 #'
 #' @param D              \strong{SummarizedExperiment} \cr
 #'                       SummarizedExperiment Dataset containing peptide 
-#'                       intensities, e.g. the result of [readMqPeptideTable] or
+#'                       intensities, e.g. the result of 
 #'                       [normalizePeptideIntensities].
 #' @param group          \strong{character factor} \cr
 #'                       The groups per sample for aggregation, if not already 
@@ -65,8 +65,9 @@
 #' @examples
 #' file <- system.file("extdata", "peptides.txt", package = "bppg")
 #' D <- readMqPeptideTable(path = file, LFQ = TRUE, remove_contaminants = FALSE)
+#' D_norm <- bppg::normalizePeptideIntensities(D)
 #' group <- factor(rep(1:9, each = 3))
-#' aggregateReplicates(D, group = group)
+#' aggregateReplicates(D_norm, group = group)
 #' 
 #' @importFrom checkmate assertCharacter assertClass assertDataFrame 
 #' assertFactor assertNumber 
@@ -78,15 +79,15 @@ aggregateReplicates <- function(D,
     method = "mean",
     seq_col = "Sequence") {
     checkmate::assertClass(D, "SummarizedExperiment")
-    checkmate::assertDataFrame(SummarizedExperiment::assays(D)$intensities,
+    checkmate::assertDataFrame(SummarizedExperiment::assays(D)$intensities_norm,
         all.missing=FALSE)
     checkmate::assertFactor(group, null.ok = TRUE)
     checkmate::assertNumber(missing.limit, lower = 0, upper = 1)
     checkmate::assertCharacter(method, pattern = "mean|sum|median")
     checkmate::assertCharacter(seq_col)
 
     id <- SummarizedExperiment::rowData(D)[, seq_col]
-    intensities <- SummarizedExperiment::assays(D)$intensities
+    intensities <- SummarizedExperiment::assays(D)$intensities_norm
 
     if (is.null(group)) {
         group <- factor(SummarizedExperiment::colData(D)$group)
@@ -175,8 +176,10 @@ calculatePeptideRatios <- function(D, group_levels = NULL) {
 
 #' Normalization of peptide intensities
 #' 
-#' @inheritParams aggregateReplicates
-#'
+#' @param D              \strong{SummarizedExperiment} \cr
+#'                       SummarizedExperiment dataset containing peptide 
+#'                       intensities, e.g. 
+#'                       [readMqPeptideTable].
 #' @param method \strong{character} \cr
 #'          The method of normalization. Options are "nonorm"
 #'          (no normalization), "median", "loess",  "quantile" or "lts"
@@ -210,7 +213,8 @@ normalizePeptideIntensities <- function(D, method = "loess", lts.quantile = 0.8,
         #### choose normalization function
         fun <- limma::normalizeBetweenArrays
         args <- switch(method,
-            "loess" = list(object = log_DATA, method = "cyclicloess"),
+            "loess" = list(object = log_DATA, method = "cyclicloess",
+                adaptive.span = FALSE),
             "quantile" = list(object = log_DATA, method = "quantile"),
             "median" = list(object = log_DATA, method = "scale"))
 

R/03_FASTAPreprocessing.R

@@ -92,7 +92,7 @@
     }
     if (missed > length(stop)) {
         if (warn){
-            warning(paste("number of specified missed cleavages is greater",
+            warning(paste0("number of specified missed cleavages is greater ",
                 "than the possible maximum"))
         }
     }

R/04_graphGeneration.R

@@ -96,41 +96,44 @@
             match(igraph::V(G)$name[!igraph::V(G)$type],
                 vMapping$peptides$peptide)])
 
-    gCollapsed <- igraph::contract(G, 
+    gColl <- igraph::contract(G, 
         factor(stats::na.omit(igraph::V(G)$collSignature)),
         vertex.attr.comb = c)
 
     # remove duplicate edges
-    gCollapsed  <- igraph::simplify(gCollapsed)
+    gColl  <- igraph::simplify(gColl)
 
     # reset attributes
-    igraph::V(gCollapsed)$type <- vapply(igraph::V(gCollapsed)$type, "[", 1,
+    igraph::V(gColl)$type <- vapply(igraph::V(gColl)$type, "[", 1,
         FUN.VALUE = logical(1))
-    igraph::V(gCollapsed)$name <- vapply(igraph::V(gCollapsed)$name,
-        paste, collapse=";", FUN.VALUE = character(1)) 
+    igraph::V(gColl)$name <- vapply(igraph::V(gColl)$name, paste, collapse=";",
+        FUN.VALUE = character(1)) 
     # this is not ordered - > same ratio order
 
 
-    if (!is.null(igraph::V(gCollapsed)$pep_logRatio)) {
+    if (!is.null(igraph::V(gColl)$pep_logRatio)) {
         if (collPeptNodes) {
-            igraph::V(gCollapsed)$pep_ratio_mean[!igraph::V(gCollapsed)$type] <-
-                vapply(
-                    igraph::V(gCollapsed)$pep_logRatio[
-                        !igraph::V(gCollapsed)$type],
+            igraph::V(gColl)$pep_ratio_mean[!igraph::V(gColl)$type] <-
+                vapply(igraph::V(gColl)$pep_logRatio[!igraph::V(gColl)$type],
                     mean, FUN.VALUE = numeric(1))
         } else {
-            igraph::V(gCollapsed)$pep_logRatio <- vapply(
-                igraph::V(gCollapsed)$pep_logRatio,  "[", 1,
-                FUN.VALUE = numeric(1))
+            igraph::V(gColl)$pep_logRatio <- vapply(
+                igraph::V(gColl)$pep_logRatio,  "[", 1, FUN.VALUE = numeric(1))
         } }
 
-    if (!is.null(igraph::V(gCollapsed)$protOrigin) && collProtNodes) {
-        igraph::V(gCollapsed)$protOrigin[igraph::V(gCollapsed)$type] <-
-            vapply(igraph::V(gCollapsed)$protOrigin[igraph::V(gCollapsed)$type],
-                unique, FUN.VALUE = character(1))
+    if (!is.null(igraph::V(gColl)$protOrigin)) {
+        if (collProtNodes) {
+            igraph::V(gColl)$protOrigin[igraph::V(gColl)$type] <- vapply(
+                    igraph::V(gColl)$protOrigin[igraph::V(gColl)$type], 
+                    function(x) { paste(unique(x), collapse = ";")
+                    }, FUN.VALUE = character(1))
+        }
+        # igraph::V(gColl)$protOrigin[!igraph::V(gColl)$type] <- NA
+        igraph::V(gColl)$protOrigin <- vapply(igraph::V(gColl)$protOrigin, "[", 
+        1, FUN.VALUE = character(1))
     }
 
-    return(igraph::delete_vertex_attr(gCollapsed, "collSignature"))
+    return(igraph::delete_vertex_attr(gColl, "collSignature"))
 }
 
 
@@ -163,8 +166,8 @@ generateGraphsFromEdgelist <- function(edgelist,
     checkmate::assertFlag(collPeptNodes)
 
     if (collProtNodes || collPeptNodes) {
-        vertexMapping <- .getContractMapping(edgelist, collProtNodes,
-            collPeptNodes)
+        vertexMapping <- .getContractMapping(
+            edgelist[, c("protein", "peptide")], collProtNodes, collPeptNodes)
     }
 
     #generate graph from edge matrix
@@ -183,6 +186,17 @@ generateGraphsFromEdgelist <- function(edgelist,
                 match(igraph::V(G)$name[!igraph::V(G)$type],
                     edgelist$peptide)])
     }
+
+    if (!is.null(edgelist$protOrigin)) {
+        protOriginDF <- edgelist[, c("protein", "protOrigin")]
+        protOriginDF <- protOriginDF[!duplicated(protOriginDF), ]
+        G <- igraph::set_vertex_attr(graph = G,
+            name = "protOrigin",
+            index = igraph::V(G)[igraph::V(G)$type],
+            value = protOriginDF$protOrigin[
+                match(igraph::V(G)$name[igraph::V(G)$type],
+                    protOriginDF$protein)])
+    }
 
     if (collProtNodes || collPeptNodes) {
         G <- .contractGraph(G, vertexMapping, collProtNodes, collPeptNodes)

R/06_optimization.R

@@ -1,5 +1,6 @@
 # Functions in this file:
 # .errorEquation
+# .trackingDataFrame
 # .initializeCi
 # .initializeRi
 # .calcConstraints
@@ -63,6 +64,28 @@
 }
 
 
+#' Funktion to format Tracking information into a dataframe with columnnames.
+#' 
+#' @param i     \strong{numeric} \cr
+#'              iterrator, for when this tracking was performed
+#' @param RES   \strong{list} \cr
+#'              The result from [.errorEquation]
+#' @param RiLog \strong{numeric vector} \cr
+#'              Contains the (estimated) protein ratios (log2-scale).
+#' @param Ci    \strong{numeric vector} \cr
+#'              Contains the protein weights (estimated, sum up to 1)
+#' @return      \strong{data.frame} \cr
+#'              combined information in one data.frame.
+.trackingDataFrame <- function(i, RES, RiLog, Ci){
+    track_colnames <- c("iter", "squ_err",
+        paste0("RLog", seq_along(RiLog)), paste0("C", seq_along(Ci)))
+    Tracking <- matrix(c(i, RES$res_squ_err, RiLog, Ci), nrow = 1)
+    Tracking <- as.data.frame(Tracking)
+    colnames(Tracking) <- track_colnames
+    return(Tracking)
+}
+
+
 
 #' Calulate initial values for Ci (protein weights) for optimization
 #'
@@ -288,7 +311,7 @@
 #' \item{RES}{final result of \code{\link[bppg]{.errorEquation}}, which also
 #' contains the final, minimal error term}
 #' \item{Tracking}{Tracking of Ri, Ci and error term for the
-#'  different iterations}
+#'  different iterations (-1: no optimization)}
 #' \item{outer.iter}{Number of outer iterations needed for the optimization
 #'  algorithm to converge or stop (see also \code{\link[Rsolnp]{solnp}})}
 #' \item{convergence}{Indicates whether the solver has converged (0) or
@@ -315,6 +338,16 @@
     rjLog <- stats::na.omit(igraph::vertex_attr(G, "pep_logRatio"))
     if (is.null(rjLog)) stop("G does not contain peptide ratios.")
     checkmate::assertNumeric(rjLog)
+    if (m == 1){
+        RiLog <- mean(rjLog, na.rm = TRUE)
+        RES <- .errorEquation(RiLog = c(RiLog),
+            Ci = c(1.0), M = M, rjLog = rjLog)
+        Tracking <- .trackingDataFrame(-1, RES, c(RiLog), c(1.0))
+        result <- list(RiLog = c(RiLog), Ci = 1.0, RES = RES, 
+            Tracking = Tracking, outer.iter = 0, convergence = 0)
+        return(result)
+    }
+
     checkmate::assertNumeric(fixedCi, len = m, lower = 0, upper = 1,
         null.ok = TRUE)
     stopifnot(sum(fixedCi, na.rm = TRUE) <= 1)
@@ -331,12 +364,7 @@
     ## initial error term
     RES <- .errorEquation(RiLog = RiLog_start, Ci = Ci_start, M = M,
         rjLog = rjLog)
-
-    track_colnames <- c("iter", "squ_err",
-        paste0("RLog", seq_len(m)), paste0("C", seq_len(m)))
-    Tracking <- matrix(c(0, RES$res_squ_err, RiLog_start, Ci_start), nrow = 1)
-    Tracking <- as.data.frame(Tracking)
-    colnames(Tracking) <- track_colnames
+    Tracking <- .trackingDataFrame(0, RES, RiLog_start, Ci_start)
 
     fun <- .calcObjectiveFunction(fixedCi, M, rjLog)
     constr <- .calcConstraints(fixedCi, m)
@@ -590,7 +618,7 @@ iterateOverCi <- function(G,
 #'                                          other proteins within the same graph
 #'                                          will be used to calculate the
 #'                                          optimal solution for each protein
-#'                                          node.
+#'                                          node. Default is TRUE
 #' @param verbose                       \strong{logical} \cr
 #'                                      If \code{TRUE}, additional information
 #'                                      will be printed.
@@ -619,7 +647,7 @@ iterateOverCi <- function(G,
 #' automatedAnalysisIteratedCi(G, res)
 automatedAnalysisIteratedCi <- function(G,
                                         res,
-                                        use_results_from_other_proteins = FALSE,
+                                        use_results_from_other_proteins = TRUE,
                                         verbose = FALSE,
                                         job = NULL,
                                         error_tol = 1e-10,

R/wrapperGraph.R

@@ -97,6 +97,10 @@ generateGraphsFromFASTA <- function(fasta,
 #' @param seq_column          \strong{character} \cr
 #'                            The column name of the column of D with the
 #'                            peptide sequences.
+#' @param norm_method         \strong{character} \cr
+#'            The method of normalization. Options are "nonorm"
+#'            (no normalization), "median", "loess",  "quantile" or "lts"
+#'             normalization. Default is "loess"
 #' @param collProtNodes       \strong{logical} \cr
 #'                            If \code{TRUE}, the protein nodes
 #'                            will be collapsed.
@@ -142,6 +146,7 @@ generateGraphsFromQuantData <- function(D,
     min_aa = 6,
     max_aa = 50,
     seq_column = "Sequence",
+    norm_method = "loess",
     collProtNodes = TRUE,
     collPeptNodes = FALSE,
     suffix = "",
@@ -161,9 +166,18 @@ generateGraphsFromQuantData <- function(D,
             overwrite = TRUE, keepNA = TRUE)
     }
 
+    ## add normalization
+    D_norm <- bppg::normalizePeptideIntensities(D, method = norm_method)
+
+    if (!is.null(outpath)) {
+        openxlsx::write.xlsx(SummarizedExperiment::assays(D_norm)$intensities_norm,
+            file = paste0(outpath, "peptides_", norm_method, "_", suffix,
+                ".xlsx"), overwrite = TRUE, keepNA = TRUE)
+    }
+
     ## aggregate replicates by calculating the mean
-    group <- factor(limma::strsplit2(colnames(D), split = "_")[, 1])
-    D_aggr <- aggregateReplicates(D, method = "mean", missing.limit = 0.4,
+    group <- factor(limma::strsplit2(colnames(D_norm), split = "_")[, 1])
+    D_aggr <- aggregateReplicates(D_norm, method = "mean", missing.limit = 0.4,
         group = group, seq_col = seq_column)
 
     if (!is.null(outpath)) {