Merge pull request #19 from wvictor14/eoPred

vignette update + function addition

.gitignore

@@ -13,7 +13,9 @@ Meta
 *hexsticker.pptx
 pkgdown/
 .vscode/settings.json
+.DS_Store
 data-raw/mod.rds
 data-raw/x_test.rds
 data-raw/x_train.rds
 data-raw/y_train.rds
+R/bioc-devel.sh

DESCRIPTION

@@ -4,8 +4,10 @@ Version: 1.15.0
 Authors@R: 
     c(person("Victor", "Yuan", email = "[email protected]",
             role = c("aut", "cre")),
-    person(c("Wendy", "P."), "Robinson",
-            role = "ctb"))
+      person(c("Wendy", "P."), "Robinson",
+            role = c("aut", "ctb")),
+      person("Icíar", "Fernández-Boyano", email = "[email protected]",
+            role = c("aut", "ctb")))
 URL: 
     https://victor.rbind.io/planet, http://github.com/wvictor14/planet
 BugReports: 
@@ -22,6 +24,8 @@ Imports:
     magrittr,
     dplyr
 Suggests:
+    ExperimentHub,
+    mixOmics,
     ggplot2,
     testthat,
     tidyr,
@@ -33,7 +37,7 @@ Suggests:
 License: GPL-2
 Encoding: UTF-8
 LazyData: false
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.2
 VignetteBuilder: knitr
 biocViews: Software, DifferentialMethylation, Epigenetics, Microarray, 
     MethylationArray, DNAMethylation, CpGIsland

NAMESPACE

@@ -1,9 +1,15 @@
 # Generated by roxygen2: do not edit by hand
 
 export("%>%")
-export(pl_infer_age)
-export(pl_infer_ethnicity)
 export(predictAge)
 export(predictEthnicity)
+export(predictPreeclampsia)
+importFrom(dplyr,filter)
+importFrom(dplyr,group_by)
+importFrom(dplyr,n)
+importFrom(dplyr,row_number)
 importFrom(magrittr,"%>%")
+importFrom(methods,hasArg)
+importFrom(methods,is)
+importFrom(stats,setNames)
 importFrom(tibble,tibble)

R/planet-deprecated.R

@@ -1,14 +0,0 @@
-#' Deprecated functions in \pkg{planet}
-#' 
-#' These functions still work but will be removed (defunct) in the next version.
-#' 
-#' \itemize{
-#'  \item \code{\link{pl_infer_ethnicity}}: This function has been renamed
-#'  \code{\link{predictEthnicity}}
-#'  
-#'  \item \code{\link{pl_infer_age}}: This function has been renamed
-#'  \code{\link{predictAge}}
-#' }
-#' 
-#' @name planet-deprecated
-NULL

R/predictAge.R

@@ -31,8 +31,6 @@
 #'     mutate(inferred_ga = predictAge(plBetas, type = "RPC"))
 #'     
 #' @export predictAge
-#' @export pl_infer_age
-#' @aliases pl_infer_age
 #' 
 predictAge <- function(betas, type = "RPC") {
     data(ageCpGs, envir = environment())
@@ -89,9 +87,4 @@ predictAge <- function(betas, type = "RPC") {
         as.vector()
 
     return(age)
-}
-
-pl_infer_age <- function(betas, type = 'RPC'){
-    .Deprecated('predictAge')
-    predictAge(betas = betas, type = type)
 }

R/predictEthnicity.R

@@ -30,9 +30,6 @@
 #' predictEthnicity(plBetas)
 #' 
 #' @export predictEthnicity 
-#' @export pl_infer_ethnicity
-#' @aliases pl_infer_ethnicity
-
 predictEthnicity <- function(betas, threshold = 0.75, force = FALSE) {
   data(ethnicityCpGs, envir=environment())
   pf <- intersect(ethnicityCpGs, rownames(betas)) 
@@ -147,8 +144,3 @@ glmnet_softmax <- function(x, ignore_labels = FALSE) {
   }
   pclass
 }
-
-pl_infer_ethnicity <- function(betas, threshold = 0.75) {
-  .Deprecated("predictEthnicity")
-  predictEthnicity(betas = betas, threshold = threshold)
-}

R/predictPreeclampsia.R

@@ -0,0 +1,94 @@
+#' @title predictPreeclampsia
+#'
+#' @description Uses 45 CpGs to predict early preeclampsia (PE delivered before or at 34 weeks of gestation)
+#' on placental DNA methylation microarray data.
+#'
+#' @details Assigns the class labels "early-PE" or "normotensive" to each sample
+#' and returns a class probability.
+#'
+#' # It is recommended that users apply beta-mixture quantile normalization (BMIQ) to their data
+#' prior to prediction. This was the normalization method used on the training data.
+#'
+#' @param betas matrix or array of methylation values on the beta scale (0, 1),
+#' where the variables are arranged in rows, and samples in columns.
+#' 
+#' @param ... feeds into outersect function
+#'
+#' @return produces a list with components detailed in the `mixOmics::predict` R documentation
+#'
+#' @examples
+#'
+#' # To predict early preeclampsia on 450k/850k samples
+#'
+#' # Load data
+#' library(ExperimentHub)
+#' eh <- ExperimentHub()
+#' query(eh, "eoPredData")
+#' 
+#' # test object
+#' x_test <- eh[['EH8403']]
+#' x_test %>% predictPreeclampsia()
+#' 
+#' @export predictPreeclampsia
+#'
+predictPreeclampsia <- function(betas, ...){
+
+  # read in data to generate model
+  eh <- ExperimentHub::ExperimentHub()
+  mod <- eh[['EH8090']]
+  trainCpGs <- colnames(mod$X)
+
+  # check that there are no NAs in the predictors (or if there are, how many)
+  peCpGs <- mixOmics::selectVar(mod)$name
+  pp <- intersect(rownames(betas), peCpGs)
+
+  if(length(pp) < length(peCpGs)){
+    stop(paste(
+      "Only", length(pp), "out of 45 predictive CpGs present. All 45 predictive CpGs are needed to run the function."
+    ))
+  } else {
+    message(paste(length(pp), "of 45 predictive CpGs present."))
+    message("BMIQ normalization is recommended for best results. If choosing other method, it is recommended to compare results to predictions on BMIQ normalized data.")
+  }
+
+  # set up data for prediction
+
+  # if input data is missing any of the cpgs present in the training data, this function
+  # adds the ones that are missing as NAs
+  # necessary for `mixOmics::predict` to work
+
+  outersect <- function(x, y) {
+    sort(c(x[!x%in%y],
+           y[!y%in%x]))
+  }
+
+  if(inherits(betas, 'matrix')){
+
+  } else if (inherits(betas, 'array')) {
+
+  } else {
+    # throw an error
+    print(paste0("Input data must be a matrix or an array"))
+  }
+
+  betasSubset <- betas[rownames(betas) %in% trainCpGs,]
+
+  # order
+  betasSubset <- betasSubset[drop=FALSE,trainCpGs, ]
+
+  stopifnot(all(rownames(betasSubset) == trainCpGs))
+
+  # predict
+  out <- predict.mixo_spls(mod, t(betasSubset))
+
+  # get class probabilities
+  CP <- out$predict[,,1]
+  CP <- t(apply(as.matrix(CP), 1, function(data) exp(data)/sum(exp(data))))
+  CP <- as.data.frame(CP) %>% tibble::rownames_to_column("Sample_ID")
+  CP$PE_Status <- CP$comp1
+  CP <- CP %>%
+    dplyr::mutate(PE_Status = dplyr::case_when(EOPE > 0.55 ~ "EOPE",
+                                               EOPE < 0.55 ~ "Normotensive"))
+  return(tibble::as_tibble(CP))
+}
+