add kmerclust_clust.R

src/bin/kmerclust_clust.R

+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"))
+