From 56bab9180a44c17b4e870339aa98e1c83d76c747 Mon Sep 17 00:00:00 2001
From: aduvermy <arnaud.duvermy@ens-lyon.fr>
Date: Tue, 6 Dec 2022 18:39:31 +0100
Subject: [PATCH] init HTRfit
---
src/model.fit/HTRfit/R/model_fitting.R | 68 +++++++++++++++++++++++++
src/model.fit/HTRfit/R/preprocessing.R | 33 ++++++++++++
src/model.fit/HTRfit/R/wald_test.R | 37 ++++++++++++++
src/model.fit/HTRfit/devtools_history.R | 13 +++++
4 files changed, 151 insertions(+)
create mode 100644 src/model.fit/HTRfit/R/model_fitting.R
create mode 100644 src/model.fit/HTRfit/R/preprocessing.R
create mode 100644 src/model.fit/HTRfit/R/wald_test.R
create mode 100644 src/model.fit/HTRfit/devtools_history.R
diff --git a/src/model.fit/HTRfit/R/model_fitting.R b/src/model.fit/HTRfit/R/model_fitting.R
new file mode 100644
index 0000000..f4d9b47
--- /dev/null
+++ b/src/model.fit/HTRfit/R/model_fitting.R
@@ -0,0 +1,68 @@
+
+#' launch glm on data
+#' @param data2fit
+#' @param model2fit
+#' @param fit_by
+#' @param threads
+#' @import future
+#' @import furrr
+#' @return fit list
+#' @export
+#'
+#' @examples
+launch.glm <- function( data2fit, model2fit = k_ij ~ genotype + environment + genotype:environment,
+ fit_by = "gene_id", threads = 4 ) {
+
+ iteration_list = data2fit[fit_by] %>% unique() %>% as.vector() %>% unname() %>% .[[1]]
+ future::plan(multisession, workers = threads)
+ fit_list = iteration_list %>% furrr::future_map(.x = .,
+ ~fit.glm( data2fit[data2fit[fit_by] == .x,],
+ .x,
+ model2fit) )
+ return(fit_list)
+
+}
+
+
+#' fit model on data
+#' @param data
+#' @param id
+#' @import MASS
+#' @import dplyr
+#' @return fit
+#' @export
+#'
+#' @examples
+fit.glm <- function(data, id, model2fit) {
+ fit = MASS::glm.nb(model2fit, data = data, link = log)
+ fit.dtf = tidySummary(fit)
+ fit.dtf$inference = fit.dtf$inference %>% dplyr::mutate(gene_id = id)
+ fit.dtf$fitQuality = fit.dtf$fitQuality %>% dplyr::mutate(gene_id = id)
+ return(fit.dtf)
+}
+
+#' fit model on data
+#' @param fit
+#' @import broom
+#' @return list of element
+#' @export
+#'
+#' @examples
+tidySummary <- function(fit){
+ dtf = broom::tidy(fit)
+ gl = broom::glance(fit)
+ return( list(inference = dtf, fitQuality = gl) )
+}
+
+#' convert list to dataframe
+#' @param fit
+#' @return list of dtf
+#' @export
+#'
+#' @examples
+listFit2dtf <- function(list_fit){
+ tmp = do.call(cbind, list_fit)
+ inference.dtf = do.call(rbind, tmp[1,])
+ fitQuality.dtf = do.call(rbind, tmp[2,])
+ return( list(inference = inference.dtf, fitQuality = fitQuality.dtf) )
+}
\ No newline at end of file
diff --git a/src/model.fit/HTRfit/R/preprocessing.R b/src/model.fit/HTRfit/R/preprocessing.R
new file mode 100644
index 0000000..045ff07
--- /dev/null
+++ b/src/model.fit/HTRfit/R/preprocessing.R
@@ -0,0 +1,33 @@
+#' get df2fit for a given gene
+#' @param k_ij
+#' @param design_simulation
+#' @param gene_name
+#' @import dplyr
+#' @return dataframe
+#' @export
+#'
+#' @examples
+getdf2fit <- function(k_ij , design_simulation, gene_name){
+ df_gene_i = cbind(design_simulation, k_ij) %>%
+ dplyr::mutate(gene_id = gene_name) %>%
+ dplyr::mutate(gene_id = gene_name)
+ rownames(df_gene_i) <- NULL
+ return(df_gene_i)
+}
+
+#' get df2fit for a given gene
+#' @param countTable
+#' @param experimental_design
+#' @import purrr
+#' @return dataframe
+#' @export
+#'
+#' @examples
+reshapeCounTable<- function(countTable, experimental_design ){
+
+ gene_id_list = rownames( countTable )
+ list_df2fit = purrr::map( .x = gene_id_list, ~getdf2fit( countTable[ .x, ],
+ experimental_design, .x) )
+ data2fit = do.call( rbind, list_df2fit )
+ return(data2fit)
+}
\ No newline at end of file
diff --git a/src/model.fit/HTRfit/R/wald_test.R b/src/model.fit/HTRfit/R/wald_test.R
new file mode 100644
index 0000000..2a173a5
--- /dev/null
+++ b/src/model.fit/HTRfit/R/wald_test.R
@@ -0,0 +1,37 @@
+#' fit model on data
+#' @param W
+#' @param threshold
+#' @param alternative
+#' @import stats
+#' @return pvalue
+#' @export
+#'
+#' @examples
+wald_test <- function( W, threshold = 0, alternative = "two.sided"){
+ ## two.sided : P ( X > threshold ) + P ( X < threshold )
+ pval = 2*(stats::pnorm( W, mean = 0, sd = 1, lower.tail = FALSE))
+
+ ## greater : P ( X > threshold )
+ pval = stats::pnorm( W, mean = 0, sd = 1, lower.tail = FALSE)
+
+ ## less : P ( X <= threshold )
+ pval = stats::pnorm( W , mean = 0, sd = 1, lower.tail = TRUE)
+ return(pval)
+}
+
+
+#' fit model on data
+#' @param estimate
+#' @param threshold
+#' @param stdError
+#' @return W
+#' @export
+#'
+#' @examples
+getStatisticWaldTest<- function(estimate,stdError, threshold ){
+ ## cf https://en.wikipedia.org/wiki/Wald_test
+ return( ( estimate - threshold ) / stdError )
+}
+
+
+
diff --git a/src/model.fit/HTRfit/devtools_history.R b/src/model.fit/HTRfit/devtools_history.R
new file mode 100644
index 0000000..269a5fe
--- /dev/null
+++ b/src/model.fit/HTRfit/devtools_history.R
@@ -0,0 +1,13 @@
+usethis::use_build_ignore("devtools_history.R")
+usethis::use_package('tidyverse', type = "depends")
+usethis::use_package("dplyr")
+usethis::use_package("purrr")
+usethis::use_package("furrr")
+usethis::use_package("MASS")
+usethis::use_package("future")
+usethis::use_package("stats")
+usethis::use_package("broom")
+
+
+
+
--
GitLab