---
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()