New parameter for limma analysis: possibility to include gender as a co-factor.

DESCRIPTION

@@ -32,4 +32,4 @@ Imports:
     methods,
     illuminaHumanv3.db,
     topconfects
-RoxygenNote: 7.0.2
+RoxygenNote: 7.1.0

R/create_design.R

@@ -13,6 +13,9 @@
 #' @param correct_for_batch A boolean indicating whether we should correct for a poential batch
 #' effect in the limma design matrix. This expects that the actual batch information exists in
 #' the object pheno_data (as field Batch). Default to FALSE.
+#' @param correct_for_gender A boolean indicating whether we should correct for sample gender
+#' in the limma design matrix. This expects that a clinical factor termed Gender exists. Default
+#' to FALSE.
 #' @param cofactor_name A string indicating a potential cofactor (included in the phenotypic data).
 #' The default is NULL, and no cofactor is considered. This will replace the age correction if
 #' both were deemed necessary.
@@ -20,17 +23,18 @@
 #' is FALSE by default.
 #' @return The limma fit.
 create_design <- function(pheno_data, limma_parameters,
-                          correct_for_age   = FALSE,
-                          correct_for_batch = FALSE,
-                          cofactor_name     = NULL,
-                          verbose           = FALSE) {
-  
+                          correct_for_age    = FALSE,
+                          correct_for_batch  = FALSE,
+                          correct_for_gender = FALSE,
+                          cofactor_name      = NULL,
+                          verbose            = FALSE) {
+
   # We extract the limma analysis parameters
   clinical_factor <- limma_parameters[[1]]
   has_block       <- !is.null(cofactor_name)
   clinical        <- factor(pheno_data[[clinical_factor]])
   rm(clinical_factor)
-  
+
   # We create the design matrix (regular case - no block).
   # Depending on the co-factors, we have several possibilities.
   # We sill decide based on (i) a block based on a co-factor (paired samples),
@@ -40,122 +44,273 @@ create_design <- function(pheno_data, limma_parameters,
   if (has_block == FALSE) {
     if (correct_for_age == FALSE) {
       if (correct_for_batch == FALSE) {
-        # We do not have a block, we do not include age, nor batch.
-        design           <- stats::model.matrix(~ 0 + clinical)
-        colnames(design) <- levels(clinical)
-        
-        # We log information.
-        if (verbose == TRUE) {
-          message(paste0("[", Sys.time(), "] Design matrix created."))
+        if (correct_for_gender == FALSE) {
+
+          # We do not have a block, we do not include age, nor batch, nor gender.
+          design           <- stats::model.matrix(~ 0 + clinical)
+          colnames(design) <- levels(clinical)
+
+          # We log information.
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created."))
+          }
+        } else {
+
+          # We do not have a block, we do not include age, nor batch, but we include gender.
+          genders          <- factor(pheno_data[["Gender"]])
+          design           <- stats::model.matrix(~ 0 + clinical + genders)
+          colnames(design) <- c(levels(clinical), levels(genders)[2:length(levels(genders))])
+
+          # We log information.
+          rm(genders)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (genders)."))
+          }
         }
       } else {
+
         # Here we have correct_for_batch = TRUE
-        # We do not have a block, we do not include age but we do include batch.
-        batches          <- factor(pheno_data[["Batch"]])
-        design           <- stats::model.matrix(~ 0 + clinical + batches)
-        colnames(design) <- c(levels(clinical), levels(batches)[2:length(levels(batches))])
-        
-        # We log information.
-        rm(batches)
-        if (verbose == TRUE) {
-          message(paste0("[", Sys.time(), "] Design matrix created (batches)."))
+        if (correct_for_gender == FALSE) {
+
+          # We do not have a block, we do not include age, nor gender but we do include batch.
+          batches          <- factor(pheno_data[["Batch"]])
+          design           <- stats::model.matrix(~ 0 + clinical + batches)
+          colnames(design) <- c(levels(clinical), levels(batches)[2:length(levels(batches))])
+
+          # We log information.
+          rm(batches)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (batches)."))
+          }
+        } else {
+
+          # We do not have a block, we do not include age but we do include batch and gender.
+          batches          <- factor(pheno_data[["Batch"]])
+          genders          <- factor(pheno_data[["Gender"]])
+          design           <- stats::model.matrix(~ 0 + clinical + batches + genders)
+          colnames(design) <- c(levels(clinical), levels(batches)[2:length(levels(batches))],
+                                levels(genders)[2:length(levels(genders))])
+
+          # We log information.
+          rm(batches, genders)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (batches, genders)."))
+          }
         }
       }
     } else {
+
       # Here we have correct_for_age = TRUE
       if (correct_for_batch == FALSE) {
-        # We do not have a block, we do include age, but not batch.
-        ages             <- as.numeric(pheno_data[["Age"]])
-        design           <- stats::model.matrix(~ 0 + clinical + ages)
-        colnames(design) <- c(levels(clinical), "Age")
-        
-        # We log information.
-        rm(ages)
-        if (verbose == TRUE) {
-          message(paste0("[", Sys.time(), "] Design matrix created (ages)."))
+        if (correct_for_gender == FALSE) {
+
+          # We do not have a block, we do include age, but not batch, nor gender.
+          ages             <- as.numeric(pheno_data[["Age"]])
+          design           <- stats::model.matrix(~ 0 + clinical + ages)
+          colnames(design) <- c(levels(clinical), "Age")
+
+          # We log information.
+          rm(ages)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (ages)."))
+          }
+        } else {
+
+          # We do not have a block, we do include age and gender but not batch.
+          ages             <- as.numeric(pheno_data[["Age"]])
+          genders          <- factor(pheno_data[["Gender"]])
+          design           <- stats::model.matrix(~ 0 + clinical + ages + genders)
+          colnames(design) <- c(levels(clinical), "Age", levels(genders)[2:length(levels(genders))])
+
+          # We log information.
+          rm(ages, genders)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (ages, genders)."))
+          }
         }
       } else {
+
         # Here we have correct_for_batch = TRUE
-        # We do not have a block, we do include age and batch.
-        ages             <- as.numeric(pheno_data[["Age"]])
-        batches          <- factor(pheno_data[["Batch"]])
-        design           <- stats::model.matrix(~ 0 + clinical + ages + batches)
-        colnames(design) <- c(levels(clinical), "Age", levels(batches)[2:length(levels(batches))])
-        
-        # We log information.
-        rm(ages, batches)
-        if (verbose == TRUE) {
-          message(paste0("[", Sys.time(), "] Design matrix created (ages, batches)."))
+        if (correct_for_gender == FALSE) {
+
+          # We do not have a block, we do include age and batch but not gender.
+          ages             <- as.numeric(pheno_data[["Age"]])
+          batches          <- factor(pheno_data[["Batch"]])
+          design           <- stats::model.matrix(~ 0 + clinical + ages + batches)
+          colnames(design) <- c(levels(clinical), "Age", levels(batches)[2:length(levels(batches))])
+
+          # We log information.
+          rm(ages, batches)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (ages, batches)."))
+          }
+        } else {
+
+          # We do not have a block, we do include age, batch and gender.
+          ages             <- as.numeric(pheno_data[["Age"]])
+          batches          <- factor(pheno_data[["Batch"]])
+          genders          <- factor(pheno_data[["Gender"]])
+          design           <- stats::model.matrix(~ 0 + clinical + ages + batches + genders)
+          colnames(design) <- c(levels(clinical), "Age", levels(batches)[2:length(levels(batches))],
+                                levels(genders)[2:length(levels(genders))])
+
+          # We log information.
+          rm(ages, batches, genders)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (ages, batches, genders)."))
+          }
         }
       }
     }
   } else {
+
     # Here we have has_block = TRUE
     if (correct_for_age == FALSE) {
       if (correct_for_batch == FALSE) {
-        # We do have a block, we do not include age, nor batch.
-        clinical_cofactor <- factor(pheno_data[[cofactor_name]]) # nolint
-        design            <- stats::model.matrix(~ 0 + clinical + clinical_cofactor)
-        colnames(design)  <- c(levels(clinical), cofactor_name)
-        
-        # We log information.
-        rm(clinical_cofactor)
-        if (verbose == TRUE) {
-          message(paste0("[", Sys.time(), "] Design matrix created (pairs)."))
+        if (correct_for_gender == FALSE) {
+
+          # We do have a block, we do not include age, nor batch, nor gender.
+          clinical_cofactor <- factor(pheno_data[[cofactor_name]]) # nolint
+          design            <- stats::model.matrix(~ 0 + clinical + clinical_cofactor)
+          colnames(design)  <- c(levels(clinical), cofactor_name)
+
+          # We log information.
+          rm(clinical_cofactor)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (pairs)."))
+          }
+        } else {
+
+          # We do have a block, we do not include age, nor batch, but we include gender.
+          clinical_cofactor <- factor(pheno_data[[cofactor_name]]) # nolint
+          genders           <- factor(pheno_data[["Gender"]])
+          design            <- stats::model.matrix(~ 0 + clinical + genders + clinical_cofactor)
+          colnames(design)  <- c(levels(clinical), levels(genders)[2:length(levels(genders))],
+                                 cofactor_name)
+
+          # We log information.
+          rm(clinical_cofactor, genders)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (pairs, genders)."))
+          }
         }
       } else {
+
         # Here we have correct_for_batch = TRUE
-        # We do have a block, we do not include age but we do include batch.
-        clinical_cofactor <- factor(pheno_data[[cofactor_name]]) # nolint
-        batches           <- factor(pheno_data[["Batch"]])
-        design            <- stats::model.matrix(~ 0 + clinical + batches + clinical_cofactor)
-        colnames(design)  <- c(levels(clinical),
-                               levels(batches)[2:length(levels(batches))],
-                               cofactor_name)
-        
-        # We log information.
-        rm(clinical_cofactor, batches)
-        if (verbose == TRUE) {
-          message(paste0("[", Sys.time(), "] Design matrix created (pairs, batches)."))
+        if (correct_for_gender == FALSE) {
+
+          # We do have a block, we do not include age, nor gender but we do include batch.
+          clinical_cofactor <- factor(pheno_data[[cofactor_name]]) # nolint
+          batches           <- factor(pheno_data[["Batch"]])
+          design            <- stats::model.matrix(~ 0 + clinical + batches + clinical_cofactor)
+          colnames(design)  <- c(levels(clinical),
+                                 levels(batches)[2:length(levels(batches))],
+                                 cofactor_name)
+
+          # We log information.
+          rm(clinical_cofactor, batches)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (pairs, batches)."))
+          }
+        } else {
+
+          # We do have a block, we do not include age, nor gender but we do include batch.
+          clinical_cofactor <- factor(pheno_data[[cofactor_name]]) # nolint
+          batches           <- factor(pheno_data[["Batch"]])
+          genders           <- factor(pheno_data[["Gender"]])
+          design            <- stats::model.matrix(~ 0 + clinical + batches + genders + clinical_cofactor)
+          colnames(design)  <- c(levels(clinical),
+                                 levels(batches)[2:length(levels(batches))],
+                                 levels(genders)[2:length(levels(genders))],
+                                 cofactor_name)
+
+          # We log information.
+          rm(clinical_cofactor, batches, genders)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (pairs, batches, genders)."))
+          }
         }
       }
     } else {
+
       # Here we have correct_for_age = TRUE
       if (correct_for_batch == FALSE) {
-        # We do have a block, we do include age, but not batch.
-        clinical_cofactor <- factor(pheno_data[[cofactor_name]]) # nolint
-        ages              <- as.numeric(pheno_data[["Age"]])
-        design            <- stats::model.matrix(~ 0 + clinical + ages + clinical_cofactor)
-        colnames(design)  <- c(levels(clinical), "Age", cofactor_name)
-        
-        # We log information.
-        rm(clinical_cofactor, ages)
-        if (verbose == TRUE) {
-          message(paste0("[", Sys.time(), "] Design matrix created (pairs, ages)."))
+        if (correct_for_gender == FALSE) {
+
+          # We do have a block, we do include age, but not batch, nor gender.
+          clinical_cofactor <- factor(pheno_data[[cofactor_name]]) # nolint
+          ages              <- as.numeric(pheno_data[["Age"]])
+          design            <- stats::model.matrix(~ 0 + clinical + ages + clinical_cofactor)
+          colnames(design)  <- c(levels(clinical), "Age", cofactor_name)
+
+          # We log information.
+          rm(clinical_cofactor, ages)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (pairs, ages)."))
+          }
+        } else {
+
+          # We do have a block, we do include age and gender but not batch.
+          clinical_cofactor <- factor(pheno_data[[cofactor_name]]) # nolint
+          ages              <- as.numeric(pheno_data[["Age"]])
+          genders           <- factor(pheno_data[["Gender"]])
+          design            <- stats::model.matrix(~ 0 + clinical + ages + genders + clinical_cofactor)
+          colnames(design)  <- c(levels(clinical), "Age",
+                                 levels(genders)[2:length(levels(genders))], cofactor_name)
+
+          # We log information.
+          rm(clinical_cofactor, ages, genders)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (pairs, ages, genders)."))
+          }
         }
       } else {
+
         # Here we have correct_for_batch = TRUE
-        # We do have a block, we do include age and batch.
-        clinical_cofactor <- factor(pheno_data[[cofactor_name]]) # nolint
-        ages              <- as.numeric(pheno_data[["Age"]])
-        batches           <- factor(pheno_data[["Batch"]])
-        design            <- stats::model.matrix(~ 0 + clinical + ages + batches +
-                                                   clinical_cofactor)
-        colnames(design)  <- c(levels(clinical),
-                               "Age",
-                               levels(batches)[2:length(levels(batches))],
-                               cofactor_name)
-        
-        # We log information.
-        rm(clinical_cofactor, ages, batches)
-        if (verbose == TRUE) {
-          message(paste0("[", Sys.time(), "] Design matrix created (pairs, ages, batches)."))
+        if (correct_for_gender == FALSE) {
+
+          # We do have a block, we do include age and batch but not gender.
+          clinical_cofactor <- factor(pheno_data[[cofactor_name]]) # nolint
+          ages              <- as.numeric(pheno_data[["Age"]])
+          batches           <- factor(pheno_data[["Batch"]])
+          design            <- stats::model.matrix(~ 0 + clinical + ages + batches +
+                                                     clinical_cofactor)
+          colnames(design)  <- c(levels(clinical),
+                                 "Age",
+                                 levels(batches)[2:length(levels(batches))],
+                                 cofactor_name)
+
+          # We log information.
+          rm(clinical_cofactor, ages, batches)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (pairs, ages, batches)."))
+          }
+        } else {
+
+          # We do have a block, we do include age, batch and gender.
+          clinical_cofactor <- factor(pheno_data[[cofactor_name]]) # nolint
+          ages              <- as.numeric(pheno_data[["Age"]])
+          batches           <- factor(pheno_data[["Batch"]])
+          genders           <- factor(pheno_data[["Gender"]])
+          design            <- stats::model.matrix(~ 0 + clinical + ages + batches + genders +
+                                                     clinical_cofactor)
+          colnames(design)  <- c(levels(clinical),
+                                 "Age",
+                                 levels(batches)[2:length(levels(batches))],
+                                 levels(genders)[2:length(levels(genders))],
+                                 cofactor_name)
+
+          # We log information.
+          rm(clinical_cofactor, ages, batches, genders)
+          if (verbose == TRUE) {
+            message(paste0("[", Sys.time(), "] Design matrix created (pairs, ages, batches, genders)."))
+          }
         }
       }
     }
   }
   rm(clinical)
-  
+
   # We return the design.
   return(design)
-}
+}
\ No newline at end of file

R/run_limma.R

@@ -14,9 +14,12 @@
 #' ignored for the time being.
 #' @param correct_for_age A boolean indicating whether we should correct for age in the limma
 #' design matrix. This expects that a clinical factor termed Age exists. Default to FALSE.
-#' @param correct_for_batch A boolean indicating whether we should correct for a poential batch
+#' @param correct_for_batch A boolean indicating whether we should correct for a potential batch
 #' effect in the limma design matrix. This expects that the actual batch information exists in
 #' the object pheno_data (as field Batch). Default to FALSE.
+#' @param correct_for_gender A boolean indicating whether we should correct for sample gender
+#' in the limma design matrix. This expects that a clinical factor termed Gender exists. Default
+#' to FALSE.
 #' @param cofactor_name A string indicating a potential cofactor (included in the phenotypic data).
 #' The default is NULL, and no cofactor is considered. This will replace the age correction if
 #' both were deemed necessary.
@@ -36,15 +39,16 @@
 #' is FALSE by default.
 #' @return The limma fit.
 run_limma <- function(pheno_data, eset, limma_parameters,
-                      correct_for_age   = FALSE,
-                      correct_for_batch = FALSE,
-                      cofactor_name     = NULL,
-                      block_key         = NULL,
-                      run_topconfects   = FALSE,
-                      output_data_dir   = "",
-                      file_prefix       = "",
-                      file_suffix       = "",
-                      verbose           = FALSE) {
+                      correct_for_age    = FALSE,
+                      correct_for_batch  = FALSE,
+                      correct_for_gender = FALSE,
+                      cofactor_name      = NULL,
+                      block_key          = NULL,
+                      run_topconfects    = FALSE,
+                      output_data_dir    = "",
+                      file_prefix        = "",
+                      file_suffix        = "",
+                      verbose            = FALSE) {
 
   # We extract the limma analysis parameters
   clinical_factor <- limma_parameters[[1]]
@@ -66,12 +70,13 @@ run_limma <- function(pheno_data, eset, limma_parameters,
   # We sill decide based on (i) a block based on a co-factor (paired samples),
   # (ii) the age clincal descriptor and (iii) the experimental batches.
   # If we need to add a block based on a co-factor.
-  design <- create_design(pheno_data        = pheno_data,
-                          limma_parameters  = limma_parameters,
-                          correct_for_age   = correct_for_age,
-                          correct_for_batch = correct_for_batch,
-                          cofactor_name     = cofactor_name,
-                          verbose           = verbose)
+  design <- create_design(pheno_data         = pheno_data,
+                          limma_parameters   = limma_parameters,
+                          correct_for_age    = correct_for_age,
+                          correct_for_batch  = correct_for_batch,
+                          correct_for_gender = correct_for_gender,
+                          cofactor_name      = cofactor_name,
+                          verbose            = verbose)
 
   # We build the contrast (only one so far).
   contrast <- eval(substitute(limma::makeContrasts(coeff1 = c1, levels = design),
@@ -95,7 +100,7 @@ run_limma <- function(pheno_data, eset, limma_parameters,
   if (has_block) {
     corfit <- limma::duplicateCorrelation(eset, design, block = pheno_data[[block_key]])
     lmfit  <- limma::lmFit(eset, design, block = pheno_data[[block_key]],
-                                 correlation = corfit$consensus)
+                                 correlation   = corfit$consensus)
     if (run_topconfects) {
       lmfit_topconfects <- limma::lmFit(eset, design_topconfects,
                                         block       = pheno_data[[block_key]],
@@ -124,6 +129,7 @@ run_limma <- function(pheno_data, eset, limma_parameters,
 
   # We now finish the TopConfects analysis if necesary.
   if (run_topconfects) {
+
     # We first run TopConfects on the dedicated lmfit object.
     confects <- topconfects::limma_confects(lmfit_topconfects, coef = 1, fdr = 0.05)
 
@@ -170,6 +176,7 @@ run_limma <- function(pheno_data, eset, limma_parameters,
       # a rank rank plot.
       limma_top <- limma::topTable(fit, coef = 1, p.value = 0.05, n = Inf)
       if (nrow(limma_top) > 0) {
+
         # Log info.
         if (verbose == TRUE) {
           message(paste0("[", Sys.time(), "] [TopConfects] ", nrow(limma_top),
@@ -199,12 +206,14 @@ run_limma <- function(pheno_data, eset, limma_parameters,
         grDevices::dev.off()
         rm(p, mdplot_fn)
       } else {
+
         # Log info.
         if (verbose == TRUE) {
           message(paste0("[", Sys.time(), "] [TopConfects] no gene identified with Limma."))
         }
       } # End of the plotting (if there are Limma genes).
     } else {
+
       # Log info.
       if (verbose == TRUE) {
         message(paste0("[", Sys.time(), "] [TopConfects] no gene identified with TopConfects."))