update package with fusen

DESCRIPTION

 Package: kmerclust
 Title: Kmerclust A Small Package To Test If The Kmer Count Correspond To A
     XY Or XO Genome
-Version: 0.0.3
+Version: 0.0.4
 Authors@R: 
     person("Laurent", "Modolo", , "laurent.modolo@ens-lyon.fr", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-7606-4110"))
@@ -30,7 +30,7 @@ Suggests:
 VignetteBuilder: 
     knitr
 Config/fusen/version: 0.5.2
-Date/Publication: 2023/08/05
+Date/Publication: 2023/08/30
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3

R/boostrap_em_constraint.R

@@ -17,16 +17,25 @@
 #' @noRd
 boostrap_EM_constraint <- function(x, sex = "XY", threshold = 1, nboot = 100, bootsize = 1000, core = 6) {
     if (sex == "XY") {
-      kappa_prior <- c(.8, .15, .05)
-      mu_prior <- list(c(0, 0), c(-0.5, 2), c(0,-2))
+      kappa_prior <- c(.7, .2, .1)
+      mu_prior <- list(c(1, 1), c(1, 1.5), c(1, .5))
+      kappa_weight <- c(.9, .9, .9)
+      mu_weight <- c(.1, .8, .8)
+      sigma_weight <- c(.001, .001, .001)
     }
     if (sex == "XO") {
       kappa_prior <- c(.85, .15)
-      mu_prior <- list(c(0, 0), c(-0.5, 2))
+      mu_prior <- list(c(1, 1), c(1, 1.5))
+      kappa_weight = c(.9, .9)
+      mu_weight <- c(.1, .8)
+      sigma_weight <- c(.001, .001)
     }
     if (sex == "OO") {
       kappa_prior <- c(.1)
-      mu_prior <- list(c(0, 0))
+      mu_prior <- list(c(1, 1))
+      kappa_weight <- c(.9)
+      mu_weight <- c(.1, .8)
+      sigma_weight <- c(.001)
     }
     parallel::mclapply(as.list(1:nboot), function(iter, x, bootsize, sex) {
         x_boot <- x %>% 
@@ -37,7 +46,10 @@ boostrap_EM_constraint <- function(x, sex = "XY", threshold = 1, nboot = 100, bo
         res <- x_boot %>% 
               EM_constraint(
                 kappa_prior = kappa_prior,
-                mu_prior = mu_prior
+                mu_prior = mu_prior,
+                kappa_weight = kappa_weight,
+                mu_weight = mu_weight,
+                sigma_weight = sigma_weight
               )
         res$bic <- compute_bic_constraint(x_boot, loglik = res$loglik, sex = sex)
         return(tibble(loglik = res$loglik, BIC = res$BIC))

R/em_constraint.R

@@ -15,12 +15,12 @@
 EM_constraint <- function(
     x,
     threshold = 1,
-    frac = .8,
-    kappa_weight = nrow(x) * .5,
-    mu_weight = nrow(x) * .01,
-    sigma_weight = nrow(x) * .0001,
-    kappa_prior = c(.8, .15, .05),
-    mu_prior = list(c(0, 0), c(-0.5, 2), c(0,-2)),
+    frac = 1,
+    kappa_weight = c(.9, .9, .9),
+    mu_weight = c(.1, .8, .8),
+    sigma_weight = c(.001, .001, .001),
+    kappa_prior = c(.7, .2, .1),
+    mu_prior = list(c(1, 1), c(1, 1.5), c(1, .5)),
     verbose = F
   ) {
   params <- EM_init_params(

R/em_init_params.R

@@ -15,35 +15,37 @@
 #' @noRd
 EM_init_params <- function(
     x,
-    kappa_weight = nrow(x) * .5,
-    mu_weight = nrow(x) * .01,
-    sigma_weight = nrow(x) * .0001,
-    kappa_prior = c(.8, .15, .05),
-    mu_prior = list(c(0, 0), c(-.5, 2), c(0,-2))
+    kappa_weight = c(.9, .9, .9),
+    mu_weight = c(.1, .8, .8),
+    sigma_weight = c(.001, .001, .001),
+    kappa_prior = c(.7, .2, .1),
+    mu_prior = list(c(1, 1), c(1, 1.5), c(1, .5))
     ) {
   x_mean <- colMeans(x)
   params <- list(
-    gamma = kappa_prior * kappa_weight + 1
+    gamma = kappa_prior * nrow(x) * kappa_weight + 1
   )
   params$rho <- list()
   params$nu <- list()
-  for (k in 1:length(params$gamma)) {
-    params$rho[[k]] <- x_mean + mu_prior[[k]]
-    params$nu[[k]] <- rep(mu_weight, ncol(x))
-  }
-  params$kappa <- as.vector(MCMCpack::rdirichlet(1, params$gamma))
-  params$kappa <- params$kappa / sum(params$kappa)
   params$beta <- list()
   params$mu <- list()
   params$sigma <- list()
   params$sigma_prior <- list()
-  # params$alpha <- rep((sigma_weight + ncol(x)) / 2, length(params$gamma))
-  params$alpha <- rep(ncol(x), length(params$gamma))
+  
+  for (k in 1:length(params$gamma)) {
+    params$rho[[k]] <- x_mean * mu_prior[[k]]
+    params$nu[[k]] <- nrow(x) * mu_weight[k]
+  }
+  params$kappa <- as.vector(MCMCpack::rdirichlet(1, params$gamma))
+  params$kappa <- params$kappa / sum(params$kappa)
+  params$alpha <- (nrow(x) * sigma_weight + ncol(x)) / 2
   for (k in 1:length(params$gamma)) {
-    params$beta[[k]] <- sigma_weight / 2 * cov(x)
-    params$sigma_prior[[k]] <- retry::retry(MCMCpack::riwish(
-      params$alpha[k], params$beta[[k]]
-    ), when = "random error", max_tries = 5)
+    params$alpha[k] <- max(c(params$alpha[k], ncol(x)))
+    params$beta[[k]] <- nrow(x) * sigma_weight[k] / 2 * cov(x)
+    # params$sigma_prior[[k]] <- retry::retry(MCMCpack::riwish(
+    #   params$alpha[k], params$beta[[k]]
+    # ), when = "random error", max_tries = 5)
+    params$sigma_prior[[k]] <- cov(x)
     params$sigma[[k]] <- params$sigma_prior[[k]]
     params$mu[[k]] <- mvtnorm::rmvnorm(1, mean = params$rho[[k]], sigma = params$nu[[k]]^-1 * params$sigma[[k]])
   }

R/plot_model.R

@@ -30,8 +30,6 @@ plot_model <- function(x, h, params) {
     chr_name <- c("A")
   }
   p1 <- x %>%
-    exp() %>% 
-    log10() %>% 
     as_tibble() %>% 
     mutate(
         clust_proba = clust_proba,

dev/0-dev_history.Rmd

@@ -21,7 +21,7 @@ fusen::fill_description(
     `Authors@R` = c(
       person("Laurent", "Modolo", email = "laurent.modolo@ens-lyon.fr", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7606-4110"))
     ),
-    Version = "0.0.3",
+    Version = "0.0.4",
     "Date/Publication" = format(Sys.time(), "%Y/%m/%d")
   ),
   overwrite = T

dev/config_fusen.yaml

@@ -2,7 +2,7 @@ flat_full.Rmd:
   path: dev/flat_full.Rmd
   state: active
   R:
-  - R/sim_kmer.R
+  - R/compute_tpm.R
   - R/expand_theta.R
   - R/parse_annotation.R
   tests: []
@@ -11,7 +11,7 @@ flat_full.Rmd:
     flat_file: dev/flat_full.Rmd
     vignette_name: Getting started
     open_vignette: false
-    check: true
+    check: false
     document: true
     overwrite: 'yes'
 flat_full_bayesian.Rmd:
@@ -41,6 +41,6 @@ flat_full_bayesian.Rmd:
     flat_file: dev/flat_full_bayesian.Rmd
     vignette_name: Getting started bayesian version
     open_vignette: false
-    check: true
+    check: false
     document: true
     overwrite: 'yes'

dev/config_not_registered.csv

 "type","path","origin"
-"R","R/compute_tpm.R","Possibly deprecated file. Please check its link with detected flat source: dev/flat_full.Rmd"
 "R","R/kmerclust-package.R","No existing source path found."
+"R","R/sim_kmer.R","Possibly deprecated file. Please check its link with detected flat source: dev/flat_full.Rmd"
 "R","R/utils-pipe.R","No existing source path found."

dev/flat_full_bayesian.Rmd

@@ -206,7 +206,7 @@ EM_init_params <- function(
   params$sigma_prior <- list()
   
   for (k in 1:length(params$gamma)) {
-    params$rho[[k]] <- x_mean
+    params$rho[[k]] <- x_mean * mu_prior[[k]]
     params$nu[[k]] <- nrow(x) * mu_weight[k]
   }
   params$kappa <- as.vector(MCMCpack::rdirichlet(1, params$gamma))
@@ -215,23 +215,15 @@ EM_init_params <- function(
   for (k in 1:length(params$gamma)) {
     params$alpha[k] <- max(c(params$alpha[k], ncol(x)))
     params$beta[[k]] <- nrow(x) * sigma_weight[k] / 2 * cov(x)
-    params$sigma_prior[[k]] <- retry::retry(MCMCpack::riwish(
-      params$alpha[k], params$beta[[k]]
-    ), when = "random error", max_tries = 5)
+    # params$sigma_prior[[k]] <- retry::retry(MCMCpack::riwish(
+    #   params$alpha[k], params$beta[[k]]
+    # ), when = "random error", max_tries = 5)
+    params$sigma_prior[[k]] <- cov(x)
     params$sigma[[k]] <- params$sigma_prior[[k]]
     params$mu[[k]] <- mvtnorm::rmvnorm(1, mean = params$rho[[k]], sigma = params$nu[[k]]^-1 * params$sigma[[k]])
   }
   return(params)
 }
-test_model_XY <- function(data) {
-  data %>%
-      dplyr::select(count_m, count_f) %>%
-      as.matrix() %>% 
-      compute_tpm() %>%
-      EM_constraint()
-}
-model_XY <- test_model_XY(data)
-plot_model(x = model_XY$x, h = model_XY$h, params = model_XY$params)
 ```
 
 ## cluster id
@@ -540,16 +532,25 @@ compute_bic_constraint <- function(x, loglik, sex = "OO") {
 #' res <- boostrap_EM(data, nboot = 10, bootsize = 0.01)
 boostrap_EM_constraint <- function(x, sex = "XY", threshold = 1, nboot = 100, bootsize = 1000, core = 6) {
     if (sex == "XY") {
-      kappa_prior <- c(.8, .15, .05)
-      mu_prior <- list(c(0, 0), c(-0.5, 2), c(0,-2))
+      kappa_prior <- c(.7, .2, .1)
+      mu_prior <- list(c(1, 1), c(1, 1.5), c(1, .5))
+      kappa_weight <- c(.9, .9, .9)
+      mu_weight <- c(.1, .8, .8)
+      sigma_weight <- c(.001, .001, .001)
     }
     if (sex == "XO") {
       kappa_prior <- c(.85, .15)
-      mu_prior <- list(c(0, 0), c(-0.5, 2))
+      mu_prior <- list(c(1, 1), c(1, 1.5))
+      kappa_weight = c(.9, .9)
+      mu_weight <- c(.1, .8)
+      sigma_weight <- c(.001, .001)
     }
     if (sex == "OO") {
       kappa_prior <- c(.1)
-      mu_prior <- list(c(0, 0))
+      mu_prior <- list(c(1, 1))
+      kappa_weight <- c(.9)
+      mu_weight <- c(.1, .8)
+      sigma_weight <- c(.001)
     }
     parallel::mclapply(as.list(1:nboot), function(iter, x, bootsize, sex) {
         x_boot <- x %>% 
@@ -560,7 +561,10 @@ boostrap_EM_constraint <- function(x, sex = "XY", threshold = 1, nboot = 100, bo
         res <- x_boot %>% 
               EM_constraint(
                 kappa_prior = kappa_prior,
-                mu_prior = mu_prior
+                mu_prior = mu_prior,
+                kappa_weight = kappa_weight,
+                mu_weight = mu_weight,
+                sigma_weight = sigma_weight
               )
         res$bic <- compute_bic_constraint(x_boot, loglik = res$loglik, sex = sex)
         return(tibble(loglik = res$loglik, BIC = res$BIC))
@@ -627,7 +631,10 @@ test_model_XO <- function(data) {
       compute_tpm() %>%
       EM_constraint(
         kappa_prior = c(.85, .15),
-        mu_prior = list(c(0, 0), c(-0.5, 2))
+        mu_prior = list(c(1, 1), c(1, 1.5)),
+        kappa_weight = c(.9, .9),
+        mu_weight = c(.1, .8),
+        sigma_weight = c(.001, .001)
       )
 }
 model_XO <- test_model_XO(data)
@@ -644,7 +651,10 @@ test_model_OO <- function(data) {
       compute_tpm() %>%
       EM_constraint(
         kappa_prior = c(1),
-        mu_prior = list(c(0, 0))
+        mu_prior = list(c(1, 1)),
+        kappa_weight = c(.9),
+        mu_weight = c(.1),
+        sigma_weight = c(.001)
       )
 }
 model_OO <- test_model_OO(data)

vignettes/getting-started-bayesian-version.Rmd

@@ -103,7 +103,10 @@ test_model_XO <- function(data) {
       compute_tpm() %>%
       EM_constraint(
         kappa_prior = c(.85, .15),
-        mu_prior = list(c(0, 0), c(-0.5, 2))
+        mu_prior = list(c(1, 1), c(1, 1.5)),
+        kappa_weight = c(.9, .9),
+        mu_weight = c(.1, .8),
+        sigma_weight = c(.001, .001)
       )
 }
 model_XO <- test_model_XO(data)
@@ -120,7 +123,10 @@ test_model_OO <- function(data) {
       compute_tpm() %>%
       EM_constraint(
         kappa_prior = c(1),
-        mu_prior = list(c(0, 0))
+        mu_prior = list(c(1, 1)),
+        kappa_weight = c(.9),
+        mu_weight = c(.1),
+        sigma_weight = c(.001)
       )
 }
 model_OO <- test_model_OO(data)
@@ -134,7 +140,6 @@ comparison <- data %>%
     dplyr::select(count_m, count_f) %>%
     as.matrix() %>% 
     compute_tpm() %>%
-    # log() %>%
     compare_models_constaint(nboot = 10, bootsize = 1000, core = 1)
 comparison %>% 
   ggplot(aes(x = name, y = loglik, fill = name)) +