---
# System prepended metadata

title: 3.1. BWCorresponding.R
tags: [Source]

---

---
tags: Source
---
# 3.1. BWCorresponding.R

"""
bw vs feature count plot

Based on the corresponding data.

TODO:
- check xdata$sample_group in shell mode

"""

library(xcms)
library(CMSITools)
library(tidyverse)
library(stringi)
library(StatTools)
register(bpstart(MulticoreParam(12))) 
source('/tmp/src/getFilesMAX.R')


args <- commandArgs(trailingOnly = TRUE)
filename_in <- args[1]
filepath_in <- args[2]
filename_out <- args[3]
filepath_out <- args[4]



nfeatures <- c()
for (bw in c(seq(0.4,2.6,by=0.2), seq(3,5,by=0.5), seq(5,10,by=2))){
  # RT aligned data
  xdata <- readRDS(file=paste(filepath_in, filename_in, "_bw", bw, "_algn.rds", sep = ""))
  
  pdp <- PeakDensityParam(sampleGroups = xdata$sample_group,
                          minFraction = 0.95, bw = bw, binSize = 0.924)
  xdata <- groupChromPeaks(xdata, param = pdp)
  #saveRDS(xdata, file=paste(filepath_out, filename_out, "_bw", bw, "_corr.rds", sep = ""))
  f <- featureValues(xdata)
  nf <- nrow(f)
  nfeatures <- c(nfeatures, nf)
  print(paste("bw=", bw, " nf=", nf))
  rm(xdata)
}

bw <- c(seq(0.4,2.6,by=0.2), seq(3,5,by=0.5), seq(5,10,by=2))
pdf("bw_nfeature.pdf")
plot(bw,nfeatures, main="NEG, number of features vs bw", ylab= "Number of feature")
dev.off()
