In tangchao7498/DecodeR: Barcode demultiplexing of native RNAs nanopore sequencing

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, fig.align = 'center')

library(PorexploreR)
library(data.table)
library(rhdf5)
library(dplyr)
library(ggplot2)
library(smoother)
library(parallel)
library(spatstat)
library(changepoint)
library(segclust2d)
library(scales)
library(pbapply)

MyChangePoint <- function(sig, MinLength = 10, ChangePoints = 68, StateStat = "Mean") {
  if(is.null(StateStat) | is.na(StateStat)) {
    stop("StateStat must be one of Mean or Median")
  }

  if(length(StateStat) != 1) {
    stop("StateStat must be one of Mean or Median")
  }

  if(!is.element(StateStat, c("Mean", "Median"))) {
    stop("StateStat must be one of Mean or Median")
  }

  cp0 <- suppressWarnings(changepoint::cpt.meanvar(data = sig, 
                                                   Q = ChangePoints, 
                                                   penalty = "Manual", 
                                                   method = "BinSeg", 
                                                   class = FALSE, 
                                                   minseglen = MinLength, 
                                                   param.estimates = FALSE, 
                                                   pen.value = 0.0001)) - 0.5
  bins <- cut(seq_along(sig), c(0, cp0, length(sig)), include.lowest = T, labels = FALSE)

  if(StateStat == "Mean") {
    bin_sig <- as.numeric(by(sig, bins, mean))
  } else {
    bin_sig <- as.numeric(by(sig, bins, median))
  }
  return(bin_sig)
}

dir_bs <- "/mnt/raid5/Personal/tangchao/project/Nanopore/BarcodeDecomplex/analysis/03.BarcodeProcessing/01.NormalBarcodeSignal/20210703"
files <- list.files(path = dir_bs, pattern = ".signal", full.names = TRUE)

pblapply(FUN = function(x) {
  load(x)
  barcode <- mclapply(barcode, function(x) MyChangePoint(sig = x, ChangePoints = 98, MinLength = 10, StateStat = "Mean"), mc.cores = 4)
  save(barcode, file = gsub("01.NormalBarcodeSignal", "02.BarcodeSignalBinExp", x))
  rm("barcode")
}, files)