vignette update + function addition

.gitignore

@@ -14,3 +14,4 @@ Meta
 pkgdown/
 *.rds
 .vscode/settings.json
+.DS_Store

DESCRIPTION

@@ -4,8 +4,10 @@ Version: 1.7.0
 Authors@R: 
     c(person("Victor", "Yuan", email = "victor.2wy@gmail.com",
             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 = "iciarfernandez@outlook.com",
+            role = c("aut", "ctb")))
 URL: 
     https://victor.rbind.io/planet, http://github.com/wvictor14/planet
 BugReports: 
@@ -21,7 +23,9 @@ Imports:
     methods,
     tibble,
     magrittr,
-    dplyr
+    dplyr,
+    ExperimentHub,
+    mixOmics
 Suggests:
     badger,
     ggplot2,

NAMESPACE

@@ -5,5 +5,6 @@ export(pl_infer_age)
 export(pl_infer_ethnicity)
 export(predictAge)
 export(predictEthnicity)
+export(predictPreeclampsia)
 importFrom(magrittr,"%>%")
 importFrom(tibble,tibble)

R/predictEthnicity.R

@@ -3,8 +3,8 @@
 #' @description Uses 1860 CpGs to predict self-reported ethnicity on placental
 #' microarray data.
 #'
-#' @details Predicts self-reported ethnicity from 3 classes: Africans, Asians,
-#' and Caucasians, using placental DNA methylation data measured on the Infinium
+#' @details Predicts self-reported ethnicity from 3 groups: African, Asian,
+#' and European, using placental DNA methylation data measured on the Infinium
 #' 450k/EPIC methylation array. Will return membership probabilities that often
 #' reflect genetic ancestry composition.
 #'

R/predictPreeclampsia.R

@@ -0,0 +1,91 @@
+#' @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
+#' # data(peBetas)
+#' # predictPreeclampsia(peBetas, dist = "max.dist")
+#' 
+#' @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(colnames(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[,colnames(betas) %in% trainCpGs]
+
+  # order
+  betasSubset <- betasSubset[drop=FALSE,, trainCpGs]
+
+  if(all(colnames(betasSubset) == trainCpGs) == FALSE){
+    stop()
+  } else
+
+    # predict
+    out <- mixOmics:::predict.mixo_spls(mod, 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(CP)
+}