removed parallel; added cpg_islands na filtering

Author: romagnolid

DESCRIPTION

@@ -4,14 +4,13 @@ Type: Package
 Title: Purity Assessment through clonal MEthylation Sites
 Description: Exploiting data from DNA methylation, this package provides the operations
     required to evaluate the purity of tumor samples. 
-Version: 0.1.0
+Version: 0.2.0
 Authors@R: c(person("Dario", "Romagnoli", email="[email protected]",
                     role=c("aut", "cre")),
              person("Matteo", "Benelli", role="aut"))
 Depends: R (>= 3.3.1)
 biocViews: ROC
 Imports:
-    assertthat,
     ROC,
     parallel
 Suggests:

NAMESPACE

@@ -2,6 +2,6 @@
 
 export(compute_AUC)
 export(compute_purity)
-export(map_to_islands)
+export(reduce_to_islands)
 export(select_informative_islands)
 export(select_informative_sites)

PAMES.Rproj

@@ -6,7 +6,7 @@ AlwaysSaveHistory: Default
 
 EnableCodeIndexing: Yes
 UseSpacesForTab: Yes
-NumSpacesForTab: 4
+NumSpacesForTab: 2
 Encoding: UTF-8
 
 RnwWeave: Sweave

R/PAMES.R

@@ -6,9 +6,9 @@
 #' The two main functions of PAMES are \code{\link{select_informative_sites}}
 #' (to retreive sites of interest), and \code{\link{compute_purity}} (to compute
 #' purity of tumor samples providing a list of sites generated by
-#' \code{\link{select_informative_sites}}). Addionally, if working with
+#' \code{\link{select_informative_sites}}). Additionally, if working with
 #' Bisulphite Sequencing data, users must first
-#' map the data with \code{\link{map_to_islands}}, than select informative
+#' reduce the data with \code{\link{reduce_to_islands}}, than select informative
 #' CpG islands with \code{\link{select_informative_islands}}, and finally
 #' estimate purity with \code{\link{compute_purity}}.
 #

R/compute_AUC.R

@@ -5,30 +5,31 @@
 #'
 #' @param tumor a matrix of beta scores generated by an Illumina BeadChip.
 #' @param control a matrix of beta scores generated by an Illumina BeadChip.
-#' @param max_NAs_fraction fraction of NAs (considered independently in tumor and
-#' control samples) above which a site will not be selected (default=0.5).
-#' @param n_threads number of cores to use
+#' @param max_NAs_frac a value between 0 and 1. Sites with a fraction of NAs
+#' major than or equal to it are excluded (default = 1; no NAs allowed).
 #' @return a vector of AUC scores
 #' @export
-compute_AUC <- function(tumor, control, max_NAs_frac=0.5){
-    stopifnot(max_NAs_frac >= 0 | max_NAs_frac <= 1)
-    stopifnot(nrow(tumor) == nrow(control))
+compute_AUC <- function(tumor, control, max_NAs_frac=1){
+  tumor <- as.matrix(tumor)
+  control <- as.matrix(control)
+  stopifnot(max_NAs_frac >= 0 | max_NAs_frac <= 1)
+  stopifnot(nrow(tumor) == nrow(control))
 
-    message(sprintf("[%s] Checking per-row NAs fraction...",  Sys.time()))
-    auc <- rep(NA, nrow(tumor))
-    valid_tumor_row_logi <- apply(tumor, 1, function(x) sum(is.na(x)) / length(x) < max_NAs_frac)
-    valid_control_row_logi <- apply(control, 1, function(x) sum(is.na(x)) / length(x) < max_NAs_frac)
-    valid_row_idx <- which(valid_tumor_row_logi & valid_control_row_logi)
+  message(sprintf("[%s] Checking per-row NAs fraction...",  Sys.time()))
+  auc <- rep(NA, nrow(tumor))
+  valid_tumor_row_logi <- apply(tumor, 1, function(x) sum(is.na(x)) / length(x) < max_NAs_frac)
+  valid_control_row_logi <- apply(control, 1, function(x) sum(is.na(x)) / length(x) < max_NAs_frac)
+  valid_row_idx <- which(valid_tumor_row_logi & valid_control_row_logi)
 
-    full_table <- cbind(tumor[valid_row_idx, ], control[valid_row_idx, ])
-    state <- c(rep(1, ncol(tumor)), rep(0, ncol(control)))
-    message(sprintf("[%s] Computing AUC ", Sys.time()))
-    valid_auc <- apply(full_table, 1, function(x){
-        non_NA_idx <- which(!is.na(x))
-        roc <- ROC::rocdemo.sca(truth=state[non_NA_idx], data=x[non_NA_idx], cutpts=seq(0, 1, .01))
-        ROC::AUC(roc)
-    })
-    auc[valid_row_idx] <- valid_auc
-    message(sprintf("[%s] Done.",  Sys.time()))
-    return(auc)
+  full_table <- cbind(tumor[valid_row_idx, ], control[valid_row_idx, ])
+  state <- c(rep(1, ncol(tumor)), rep(0, ncol(control)))
+  message(sprintf("[%s] Computing AUC ", Sys.time()))
+  valid_auc <- apply(full_table, 1, function(x){
+    non_NA_idx <- which(!is.na(x))
+    roc <- ROC::rocdemo.sca(truth=state[non_NA_idx], data=x[non_NA_idx], cutpts=seq(0, 1, .01))
+    ROC::AUC(roc)
+  })
+  auc[valid_row_idx] <- valid_auc
+  message(sprintf("[%s] Done.",  Sys.time()))
+  return(auc)
 }

R/compute_islands_indexes.R

@@ -0,0 +1,31 @@
+#' Associated to each cpg islands the corresponding cpg site
+#' @param cpg_sites a data.frame reporting the location of CpG sites from Bisulphite
+#' Sequencing data and requiring two columns: chromosome and 1-based coordinate of CpG sites.
+#' @param cpg_islands_df a data.frame reporting the location of CpG islands. By default it
+#' loads the `cpg_islands` data.frame that comes with the package (based on data
+#' downloaded from [UCSC Genome Browser](https://genome.ucsc.edu/) on date 2016-03-01.
+#' @examples
+#' \dontrun{
+#' cpg_indexes <- compute_islands_indexes(bs_toy_sites)
+#' }
+#' \donttest{
+#' cpg_indexes_updated <- compute_islands_indexes(bs_toy_sites, cpg_islands_df = cpg_islands_v2)
+#' }
+compute_island_indexes <- function(cpg_sites, cpg_islands_df = cpg_islands){
+  if (!is.data.frame(cpg_sites) | !is.data.frame(cpg_islands_df))
+    stop("'cpg_sites' and 'cpg_islands_df' must be data.frames")
+  if (!is.character(cpg_sites[[1]]) | !is.numeric(cpg_sites[[2]]))
+    stop("First to colums of 'cpg_sites' must be 'character' and 'numeric', respectively")
+  if (any(!startsWith(cpg_sites[[1]], "chr")))
+    stop("Please be sure that chromosome symbols in 'cpg_sites' has the format 'chrN'")
+  if (!is.character(cpg_islands_df[[1]]) | !is.numeric(cpg_islands_df[[2]])| !is.numeric(cpg_islands_df[[3]]))
+    stop("First column of 'cpg_islands_df' must be 'character', second and third 'numeric'")
+
+  message(sprintf("[%s] Computing indexes for 'CpG sites to CpG islands' mapping...",  Sys.time()))
+  cpg_indexes <- apply(cpg_islands_df, 1, function(x) {
+    same_chromosome  <- cpg_sites[[1]] == x[["chr"]]
+    is_between <- cpg_sites[[2]] >= x[["start"]] & cpg_sites[[2]] <= x[["end"]]
+    return(which(same_chromosome & is_between))
+  })
+  return(cpg_indexes)
+}

R/compute_purity.R

@@ -16,9 +16,9 @@
 #' }
 compute_purity <- function(tumor, list_of_sites) {
     tumor <- as.matrix(tumor)
-    assertthat::assert_that(is.list(list_of_sites),
-                length(list_of_sites) == 2,
-                all(names(list_of_sites) %in% c("hyper", "hypo")) )
+    stopifnot(is.list(list_of_sites))
+    stopifnot(length(list_of_sites) == 2)
+    stopifnot(all(names(list_of_sites) %in% c("hyper", "hypo")))
 
     beta_values <- rbind(tumor[list_of_sites[["hyper"]],],
                          1 - tumor[list_of_sites[["hypo"]],])

R/data.R

@@ -24,26 +24,24 @@
 
 #' Curated set of CpG islands
 #'
-#' Automatically loaded by \code{\link{map_to_islands}}.
-#'
 #' Downloaded on 2016-03-01.
 #' @source \url{https://genome.ucsc.edu/}
 "cpg_islands"
 
 
 #' Toy data
 #'
-#' Test with \code{\link{map_to_islands}}.
+#' Test with \code{\link{reduce_to_islands}}.
 "bs_tumor_toy_data"
 
 #' Toy data
 #'
-#' Test with \code{\link{map_to_islands}}.
+#' Test with \code{\link{reduce_to_islands}}.
 "bs_control_toy_data"
 
 #' Toy data
 #'
-#' Test with \code{\link{map_to_islands}}.
+#' Test with \code{\link{reduce_to_islands}}.
 "bs_toy_sites"
 
 

R/map_to_islands.R

@@ -0,0 +1,50 @@
+#' Convert beta values from Bisulphite Sequencing to single CpG island beta values
+#'
+#' Reduce beta values from different CpG sites to single beta values
+#' associated to one CpG island where CpG sites are located.
+#'
+#' Bisulphite Sequencing is used to retrieve level of DNA methylation at
+#' base resolution. Different experiments may retrieve
+#' different CpG sites. This function makes a direct comparison possible.
+#'
+#' @param beta_values a matrix of beta values.
+#' @param cpg_indexes a list of indexes generated by \code{compute_islands_indexes}.
+#' @param min_CpGs an integer (default to 3). Minimum number of CpG sites within
+#' a single CpG island required to compute the reduced beta value (return NA otherwise).
+#' @return a matrix of beta values (one row per CpG island).
+#' @examples
+#' \donttest{
+#' reduced <- reduce_to_islands(bs_control_toy_data, bs_toy_sites)
+#' }
+#' @export
+reduce_to_islands <- function(cpg_sites_matrix, cpg_indexes, min_CpGs = 3){
+  cpg_sites_matrix <- as.matrix(cpg_sites_matrix)
+  message(sprintf("[%s] Reducing beta values...",  Sys.time()))
+  cpg_islands_matrix <- do.call("rbind",
+    lapply(cpg_indexes, function(x)
+      reduce_island(cpg_sites_matrix[x, , drop = F], min_CpGs)))
+  message(sprintf("[%s] Done",  Sys.time()))
+  return(cpg_islands_matrix)
+}
+
+
+#' Reduce a matrix subset
+#'
+#' Takes into consideration the eventual presence of NAs
+#' @param cpg_sites_matrix a matrix
+#' @param n minimum required number of sites per island (return NA otherwise)
+#' @keywords internal
+reduce_island <- function(cpg_sites_matrix, n) {
+  if (ncol(cpg_sites_matrix) < n) {
+    beta_values <- NA
+  } else {
+    # identify sites without any beta-score
+    all_NAs_sites <- rowSums(is.na(cpg_sites_matrix)) <= ncol(cpg_sites_matrix)
+    if (sum(!all_NAs_sites) < n) {
+      beta_values <- NA
+    } else {
+      beta_values <- apply(cpg_sites_matrix, 2, median, na.rm = T)
+    }
+  }
+  return(beta_values)
+}

R/reduce_to_islands.R

@@ -18,7 +18,8 @@
 #' info_sites <- select_informative_islands(mapped_bs_tumor_toy_data, auc)
 #' }
 #' @export
-select_informative_islands <- function(tumor, auc, max_sites = 20){
+select_informative_islands <- function(tumor, auc, hyper_range = c(min = .40, max = .90),
+  hypo_range = c(min = .10, max = .60), max_sites = 20){
     # check parameters ---------------------------------------------------------
     if (max_sites %% 2 != 0 & !is.integer(max_sites))
         stop("'max_sites' must be even")
@@ -30,8 +31,8 @@ select_informative_islands <- function(tumor, auc, max_sites = 20){
     # minimum and maximum beta per site ----------------------------------------
     beta_max <- suppressWarnings(apply(tumor, 1, max, na.rm=T))
     beta_min <- suppressWarnings(apply(tumor, 1, min, na.rm=T))
-    hyper_idx <- which(beta_min < .40 & beta_max > .90 & auc > .80)
-    hypo_idx  <- which(beta_min < .10 & beta_max > .60 & auc < .20)
+    hyper_idx <- which(beta_min < hyper_range[1] & beta_max > hyper_range[2] & auc > .80)
+    hypo_idx  <- which(beta_min < hypo_range[1]  & beta_max > hypo_range[2]  & auc < .20)
 
     message(sprintf("[%s] Total hyper-methylated sites = %i", Sys.time(), length(hyper_idx)))
     message(sprintf("[%s] Total hypo-methylated sites = %i",  Sys.time(), length(hypo_idx)))