From e1e76c4690c9a8d61d2262f9a8b3f9ee8960f093 Mon Sep 17 00:00:00 2001
From: Laurent Modolo <laurent.modolo@ens-lyon.fr>
Date: Thu, 9 Nov 2023 15:57:54 +0100
Subject: [PATCH] add kmerclust_clust.R
---
src/bin/kmerclust_clust.R | 34 ++++++++++++++++++++++++++++++++++
1 file changed, 34 insertions(+)
create mode 100644 src/bin/kmerclust_clust.R
diff --git a/src/bin/kmerclust_clust.R b/src/bin/kmerclust_clust.R
new file mode 100644
index 0000000..346e6a8
--- /dev/null
+++ b/src/bin/kmerclust_clust.R
@@ -0,0 +1,34 @@
+library(kmerclust)
+library(tidyverse)
+args <- commandArgs(trailingOnly = TRUE)
+print(args)
+
+load(file = paste0(args[1], ".Rdata"))
+
+params = list(
+ A = list(kappa = 0.4, l1 = mean(count$count_m), l2 = mean(count$count_m), l3 = mean(count$count_m)),
+)
+if (args[2] %in% c("XY", "XO")) {
+ params$X <- list(kappa = 0.3, l1 = 1, l2 = mean(count$count_m), l3 = mean(count$count_m))
+}
+if (args[2] == "XY") {
+ params$Y <- list(kappa = 0.3, l1 = mean(count$count_m), l2 = 1, l3 = mean(count$count_m))
+}
+
+res <- count %>%
+ dplyr::select(count_m, count_f) %>%
+ mutate(
+ count_m = round(count_m),
+ count_f = round(count_f)
+ ) %>%
+ dplyr::filter(count_m + count_f > 0) %>%
+ as.matrix() %>%
+ em_bipoiss_clust(nbatch = 100, max_iter = 1000)
+
+save(res, file = paste0(args[1], "_clust_", args[3], ".Rdata"))
+
+count %>%
+ dplyr::mutate(chromosome = res$chromosome) %>%
+ dplyr::select(name, chromsome) %>%
+ write_csv2(file = paste0(args[1], "_clust_", args[3], ".csv"))
+
--
GitLab