change variable names

bcm-uga · Apr 8, 2024 · ffc5aa7 · ffc5aa7
1 parent 99096d8
commit ffc5aa7
Show file tree

Hide file tree

Showing 22 changed files with 471 additions and 1,087 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -21,9 +21,7 @@ Imports:
     data.table,
     Biostrings,
     IRanges,
-    gridExtra,
-    fdrtool,
-    LEA
+    gridExtra
 Depends: 
     R (>= 2.10)
 VignetteBuilder: knitr

diff --git a/NAMESPACE b/NAMESPACE
@@ -4,6 +4,7 @@ export(AMR_build)
 export(AMR_search)
 export(combp2)
 export(estimate_effect)
+export(hdmax2_qvalue)
 export(lfmm2_med)
 export(lfmm2_med_test)
 export(plot_hdmax2)

diff --git a/R/hdmax2_helper_functions.R b/R/hdmax2_helper_functions.R
@@ -17,15 +17,15 @@
 ##' differentially methylated regions.
 ##' @author Basile Jumentier
 ##' @example 
-##'data = hdmax2::helper_ex
-##'chr = data$annotation$chr
-##'start = data$annotation$start
-##'end = data$annotation$end
-##'pval = hdmax2_step1$max2_pvalues
-##'cpg = data$annotation$cpg
-##'data_combp <- data.frame(chr, start, end, pval, cpg)
-##'colnames(data_combp) <- paste0("V", 1:5)
-##'res <- combp2(data_combp, seed = 0.6,  nCores = 2, ...)
+##' data = hdmax2::helper_ex
+##' chr = data$annotation$chr
+##' start = data$annotation$start
+##' end = data$annotation$end
+##' pval = hdmax2_step1$max2_pvalues
+##' cpg = data$annotation$cpg
+##' data_combp <- data.frame(chr, start, end, pval, cpg)
+##' colnames(data_combp) <- paste0("V", 1:5)
+##' res <- combp2(data_combp, seed = 0.6,  nCores = 2, ...)
 ##'
 combp2 <- function (data, dist.cutoff = 1000, bin.size = 310, seed = 0.01, nCores = 10) {
 
@@ -154,34 +154,35 @@ combp2 <- function (data, dist.cutoff = 1000, bin.size = 310, seed = 0.01, nCore
 ##' @export
 ##' @author Basile Jumentier
 ##' @examples
-##'data = hdmax2::helper_ex
-##'#Artificial reduction of dataset size to pass the github action check when building hdmax2 website
-##'data$methylation = data$methylation[ , 800:1000]
-##'data$annotation = data$annotation[800:1000, ]
-##'K=5
-##'## run hdmax2 step1
-##'hdmax2_step1 = hdmax2::run_AS(X = data$exposure,
-##'                              Y = data$phenotype,
+##' data = hdmax2::helper_ex
+##' #Artificial reduction of dataset size to pass the github action check when building hdmax2 website
+##' data$methylation = data$methylation[ , 800:1000]
+##' data$annotation = data$annotation[800:1000, ]
+##' K=5
+##' ## run hdmax2 step1
+##' hdmax2_step1 = hdmax2::run_AS(exposure = data$exposure,
+##'                              outcome = data$phenotype,
 ##'                              M = data$methylation,
 ##'                              K = K)
 ##'
-##'##Detecting AMR
-##'chr = data$annotation$chr
-##'start = data$annotation$start
-##'end = data$annotation$end
-##'pval = hdmax2_step1$max2_pvalues
-##'cpg = data$annotation$cpg
+##' ##Detecting AMR
+##' chr = data$annotation$chr
+##' start = data$annotation$start
+##' end = data$annotation$end
+##' pval = hdmax2_step1$max2_pvalues
+##' cpg = data$annotation$cpg
 ##'
-##'res.amr_search = hdmax2::AMR_search(
+##' res.amr_search = hdmax2::AMR_search(
 ##' chr = data$annotation$chr,
 ##' start = data$annotation$start,
 ##' end = data$annotation$end,
 ##' pval = hdmax2_step1$max2_pvalues,
 ##' cpg = data$annotation$cpg,
 ##' seed = 0.6, #Careful to change this parameter when working with real data
 ##' nCores = 2)
-##'res.amr_search$res
-##'##'
+##' res.amr_search$res
+##'
+
 AMR_search <- function(chr, start, end, pval, cpg, ...) {
 
   tmp <- data.frame(chr, start, end, pval, cpg)
@@ -217,25 +218,25 @@ AMR_search <- function(chr, start, end, pval, cpg, ...) {
 ##' @export
 ##' @author Basile Jumentier
 ##' @examples
-##'data = hdmax2::helper_ex
-##'#Artificial reduction of dataset size to pass the github action check when building hdmax2 website
-##'data$methylation = data$methylation[ , 800:1000]
-##'data$annotation = data$annotation[800:1000, ]
-##'K=5
-##'## run hdmax2 step1
-##'hdmax2_step1 = hdmax2::run_AS(X = data$exposure,
-##'                              Y = data$phenotype,
+##' data = hdmax2::helper_ex
+##' #Artificial reduction of dataset size to pass the github action check when building hdmax2 website
+##' data$methylation = data$methylation[ , 800:1000]
+##' data$annotation = data$annotation[800:1000, ]
+##' K=5
+##' ## run hdmax2 step1
+##' hdmax2_step1 = hdmax2::run_AS(exposure = data$exposure,
+##'                              outcome = data$phenotype,
 ##'                              M = data$methylation,
 ##'                              K = K)
 ##'
 ##'##Detecting AMR
-##'chr = data$annotation$chr
-##'start = data$annotation$start
-##'end = data$annotation$end
-##'pval = hdmax2_step1$max2_pvalues
-##'cpg = data$annotation$cpg
+##' chr = data$annotation$chr
+##' start = data$annotation$start
+##' end = data$annotation$end
+##' pval = hdmax2_step1$max2_pvalues
+##' cpg = data$annotation$cpg
 ##'
-##'res.amr_search = hdmax2::AMR_search(
+##' res.amr_search = hdmax2::AMR_search(
 ##' chr = data$annotation$chr,
 ##' start = data$annotation$start,
 ##' end = data$annotation$end,
@@ -244,14 +245,14 @@ AMR_search <- function(chr, start, end, pval, cpg, ...) {
 ##' seed = 0.6, #Careful to change this parameter when working with real data
 ##' nCores = 2)
 ##'
-##'res.amr_search$res
+##' res.amr_search$res
 ##'
-##'res.arm_build = hdmax2::AMR_build(res.amr_search, 
-##'methylation = data$methylation, nb_cpg = 2)
-##'#List of DMR selected
-##'head(res.arm_build$res)
-##'## CpG in the DMR
-##'res.arm_build$CpG_for_each_AMR
+##' res.arm_build = hdmax2::AMR_build(res.amr_search, 
+##' methylation = data$methylation, nb_cpg = 2)
+##' #List of DMR selected
+##' head(res.arm_build$res)
+##' ## CpG in the DMR
+##' res.arm_build$CpG_for_each_AMR
 ##'
 AMR_build <- function(res, methylation, nb_cpg = 2) {
 
@@ -312,3 +313,40 @@ AMR_build <- function(res, methylation, nb_cpg = 2) {
               CpG_for_each_AMR = DMR.select))
 }
 
+
+
+##' compute qvalues from pvalues
+##'
+##' @param pvalues result of max squared test
+##' @param theta constant defines interquartile range
+##' @return qvalues, pi0 =  H0_proportion
+##' @export
+##' @author Olivier Francois
+##' @examples
+##' data = hdmax2::simu_data
+##' K = 5
+##' hdmax2_step1 = hdmax2::run_AS(exposure = simu_data$X_binary,
+##'                               outcome = simu_data$Y_continuous,
+##'                               M = simu_data$M1, 
+##'                               K = K)
+##' 
+##'  #Select candidate mediator  
+##' qv = hdmax2::hdmax2_qvalue(hdmax2_step1$max2_pvalues)
+##' 
+##' 
+##' 
+hdmax2_qvalue <- function(pvalues, theta = 0.25){
+
+  ##Estimate the proportion of H0
+  ##relies on the interquartile range (25%-75%) of the default histogram
+  H0_proportion = 2*mean( pvalues > theta & pvalues < (theta + 0.5) )
+
+  ## compute tests number
+  test_number = length(pvalues)
+
+  ## qvalues definition
+  qvalues =  H0_proportion * test_number * pvalues / rank(pvalues)
+
+  return(list(qvalues =  qvalues , pi0 =  H0_proportion))
+}
+
diff --git a/R/hdmax2_lfmm2_functions.R b/R/hdmax2_lfmm2_functions.R
@@ -1,16 +1,18 @@
 ##' lfmm2 adapted function for hdmax2 from LEA::lfmm2
-##' @param input a response variable matrix with n rows and p columns
+##' @param input Continuous intermediary variables matrix  encompassing potential mediators with n rows and p columns.
 ##' @param env An explanatory variable matrix with n rows and d columns.
 ##' @param K latent factor number
 ##' @param lambda ridge penalization parameter
 ##' @param effect.sizes true or false to obtain effect sizes
-##' @return an object with the following attributes 
+##' @return an object with the following attributes
+##'  U
+##'  V
 ##' @export 
 ##' @author Florence Pittion, Magali Richard, Olivier Francois
 ##' @examples
 ##' data(simu_data)
 ##' K = 5
-##' mod.lfmm1 = lfmm2_med(input = simu_data$M, 
+##' mod.lfmm1 = lfmm2_med(input = simu_data$M1, 
 ##' env = simu_data$X_binary, 
 ##' K = K,
 ##' effect.sizes = FALSE)
@@ -143,20 +145,21 @@ lfmm2_med = function(input,
 ##' @param family of logistic reg
 ##' @param full compute partial regression FALSE/TRUE
 ##' @return an object with the following attributes 
+##' P-values, fscores, zscores, adj.squared, gif
 ##' @importFrom stats binomial glm lm median pchisq pf prcomp qchisq qf
 ##' @importFrom utils read.table
 ##' @export
 ##' @author Florence Pittion, Magali Richard, Olivier Francois
 ##' @examples 
-##' data(simu_data)
+##' attach(simu_data)
 ##' K = 5
-##' mod.lfmm1 = lfmm2_med(input = simu_data$M, 
+##' mod.lfmm1 = lfmm2_med(input = simu_data$M1, 
 ##' env = simu_data$X_binary, 
 ##' K = K,
 ##' effect.sizes = FALSE)
 ##' 
 ##' res_reg1 = lfmm2_med_test(mod.lfmm1, 
-##' input = simu_data$M, 
+##' input = simu_data$M1, 
 ##' env = simu_data$X_binary,
 ##' covar = cbind(simu_data$age, simu_data$gender),
 ##' genomic.control = TRUE)