Merge pull request #39 from UMCUGenetics/change_to_5ppm

Change to 5ppm
UMCUGenetics · Dec 8, 2020 · 296b292 · 296b292
2 parents 56e1619 + be5f897
commit 296b292
Show file tree

Hide file tree

Showing 7 changed files with 18 additions and 14 deletions.
diff --git a/pipeline/run.sh b/pipeline/run.sh
@@ -18,7 +18,7 @@ restart=0
 indir=""
 outdir=""
 rscript="/hpc/local/CentOS7/dbg_mz/R_libs/3.6.2/bin/Rscript" #temp
-ppm=2 #temp
+ppm=5 #temp
 
 # Show usage information
 function show_help() {
@@ -352,7 +352,7 @@ for file in ${outdir}/7-specpks_all_rest/${scanmode}_* ; do
   # 8-peakGrouping.rest
   echo "#!/bin/sh
 
-  /hpc/local/CentOS7/common/lang/R/3.2.2/bin/Rscript ${scripts}/8-peakGrouping.rest.R \$file ${outdir} ${scanmode} ${resol}
+  /hpc/local/CentOS7/common/lang/R/3.2.2/bin/Rscript ${scripts}/8-peakGrouping.rest.R \$file ${outdir} ${scanmode} ${resol} ${ppm}
   " > ${outdir}/jobs/8-peakGrouping.rest/${scanmode}_\${input}.sh
   cur_id=\$(sbatch --job-name=8-peakGroupingRest_\${input}_${scanmode}_${name} --time=${job_8_time} --mem=${job_8_mem} --output=${outdir}/logs/8-peakGrouping.rest/${scanmode}_\${input}.o --error=${outdir}/logs/8-peakGrouping.rest/${scanmode}_\${input}.e ${global_sbatch_parameters} ${outdir}/jobs/8-peakGrouping.rest/${scanmode}_\${input}.sh)
   job_ids+="\${cur_id}:"
@@ -374,7 +374,7 @@ for file in ${outdir}/8-grouping_rest/${scanmode}_* ; do
   # 9-runFillMissing.R part 1
   echo "#!/bin/sh
 
-  /hpc/local/CentOS7/common/lang/R/3.2.2/bin/Rscript ${scripts}/9-runFillMissing.R \$file ${outdir} ${scanmode} ${thresh} ${resol} ${scripts}
+  /hpc/local/CentOS7/common/lang/R/3.2.2/bin/Rscript ${scripts}/9-runFillMissing.R \$file ${outdir} ${scanmode} ${thresh} ${resol} ${scripts} ${ppm}
   " > ${outdir}/jobs/9-runFillMissing/rest_${scanmode}_\${input}.sh
   cur_id=\$(sbatch --job-name=9-runFillMissing_1_\${input}_${scanmode}_${name} --time=${job_9a_time} --mem=${job_9a_mem} --output=${outdir}/logs/9-runFillMissing/rest_${scanmode}_\${input}.o --error=${outdir}/logs/9-runFillMissing/rest_${scanmode}_\${input}.e ${global_sbatch_parameters} ${outdir}/jobs/9-runFillMissing/rest_${scanmode}_\${input}.sh)
   job_ids+="\${cur_id}:"
@@ -396,7 +396,7 @@ job_ids=\${job_ids::-1}
 # 10-collectSamplesFilled
 echo "#!/bin/sh
 
-/hpc/local/CentOS7/common/lang/R/3.2.2/bin/Rscript ${scripts}/10-collectSamplesFilled.R ${outdir} ${scanmode} ${normalization} ${scripts} ${z_score}
+/hpc/local/CentOS7/common/lang/R/3.2.2/bin/Rscript ${scripts}/10-collectSamplesFilled.R ${outdir} ${scanmode} ${normalization} ${scripts} ${z_score}  ${ppm}
 " > ${outdir}/jobs/10-collectSamplesFilled/${scanmode}.sh
 col_id=\$(sbatch --job-name=10-collectSamplesFilled_${scanmode}_${name} --time=${job_10_time} --mem=${job_10_mem} --dependency=afterany:\${job_ids} --output=${outdir}/logs/10-collectSamplesFilled/${scanmode}.o --error=${outdir}/logs/10-collectSamplesFilled/${scanmode}.e ${global_sbatch_parameters} ${outdir}/jobs/10-collectSamplesFilled/${scanmode}.sh)
 

diff --git a/pipeline/scripts/10-collectSamplesFilled.R b/pipeline/scripts/10-collectSamplesFilled.R
@@ -14,6 +14,7 @@ scanmode <- cmd_args[2]
 normalization <- cmd_args[3]
 scripts <- cmd_args[4]
 z_score <- as.numeric(cmd_args[5])
+ppm <- as.numeric(cmd_args[6])
 
 #outdir <- "/Users/nunen/Documents/Metab/processed/test_old"
 #scanmode <- "negative"
@@ -72,6 +73,7 @@ outlist.not.ident = outlist.tot[-index,]
 
 if (z_score == 1) {
   outlist.ident$ppmdev=as.numeric(outlist.ident$ppmdev)
+  outlist.ident <- outlist.ident[which(outlist.ident["ppmdev"] >= -ppm & outlist.ident["ppmdev"] <= ppm),]
 }
 # NAs in theormz_noise <======================================================================= uitzoeken!!!
 outlist.ident$theormz_noise[which(is.na(outlist.ident$theormz_noise))] = 0

diff --git a/pipeline/scripts/7-collectSamplesGroupedHMDB.R b/pipeline/scripts/7-collectSamplesGroupedHMDB.R
@@ -66,8 +66,8 @@ if (n>=sub & (floor(n/sub)-1)>=2){
 
       n_moved = 0
 
-      # Calculate 3ppm and replace border, avoid cut within peakgroup!
-      while ((as.numeric(outlist_i[1,"mzmed.pkt"]) - as.numeric(outlist_i_min_1[min_1_last,"mzmed.pkt"]))*1e+06/as.numeric(outlist_i[1,"mzmed.pkt"]) < 3) {
+      # Calculate ppm and replace border, avoid cut within peakgroup!
+      while ((as.numeric(outlist_i[1,"mzmed.pkt"]) - as.numeric(outlist_i_min_1[min_1_last,"mzmed.pkt"]))*1e+06/as.numeric(outlist_i[1,"mzmed.pkt"]) < ppm) {
         outlist_i_min_1 = rbind(outlist_i_min_1, outlist_i[1,])
         outlist_i = outlist_i[-1,]
         n_moved = n_moved + 1
@@ -92,8 +92,8 @@ outlist_i = outlist[c(start:end),]
 n_moved = 0
 
 if(!is.null(outlist_i_min_1)){
-  # Calculate 4ppm and replace border, avoid cut within peakgroup!
-  while ((as.numeric(outlist_i[1,"mzmed.pkt"]) - as.numeric(outlist_i_min_1[min_1_last,"mzmed.pkt"]))*1e+06/as.numeric(outlist_i[1,"mzmed.pkt"]) < 2*ppm) {
+  # Calculate ppm and replace border, avoid cut within peakgroup!
+  while ((as.numeric(outlist_i[1,"mzmed.pkt"]) - as.numeric(outlist_i_min_1[min_1_last,"mzmed.pkt"]))*1e+06/as.numeric(outlist_i[1,"mzmed.pkt"]) < ppm) {
     outlist_i_min_1 = rbind(outlist_i_min_1, outlist_i[1,])
     outlist_i = outlist_i[-1,]
     n_moved = n_moved + 1

diff --git a/pipeline/scripts/8-peakGrouping.rest.R b/pipeline/scripts/8-peakGrouping.rest.R
@@ -13,13 +13,14 @@ fileIn <- cmd_args[1]
 outdir <- cmd_args[2]
 scanmode <- cmd_args[3]
 resol <- as.numeric(cmd_args[4])
+ppm <- as.numeric(cmd_args[5])
 
 # create output folder
 dir.create(paste(outdir, "8-grouping_rest", sep = "/"),showWarnings = F)
 
 options(digits=16)
 
-groupingRest <- function(outdir, fileIn, scanmode, ppm=2) {
+groupingRest <- function(outdir, fileIn, scanmode, ppm) {
   # fileIn="./results/specpks_all_rest/negative_outlist_i_min_1.1.RData"
   # scanmode="negative"
   # outdir="./results"
@@ -95,4 +96,4 @@ groupingRest <- function(outdir, fileIn, scanmode, ppm=2) {
 
 cat("File to group: ", fileIn)
 
-groupingRest(outdir, fileIn, scanmode)
+groupingRest(outdir, fileIn, scanmode, ppm=ppm)
diff --git a/pipeline/scripts/9-runFillMissing.R b/pipeline/scripts/9-runFillMissing.R
@@ -15,6 +15,7 @@ scanmode <- cmd_args[3]
 thresh  <- as.numeric(cmd_args[4])
 resol <- as.numeric(cmd_args[5])
 scripts <- cmd_args[6]
+ppm <- as.numeric(cmd_args[7])
 
 # file="./results/grouping_rest/negative_1.RData"
 # file="./results/grouping_hmdb/1_negative.RData"
@@ -30,4 +31,4 @@ scripts <- cmd_args[6]
 source(paste(scripts, "AddOnFunctions/sourceDir.R", sep="/"))
 sourceDir(paste(scripts, "AddOnFunctions", sep="/"))
 
-replaceZeros(file,scanmode,resol,outdir,thresh,scripts)
+replaceZeros(file,scanmode,resol,outdir,thresh,scripts,ppm)
diff --git a/pipeline/scripts/AddOnFunctions/ident.hires.noise.HPC.R b/pipeline/scripts/AddOnFunctions/ident.hires.noise.HPC.R
@@ -1,5 +1,5 @@
 # modified identify function to also look for adducts and their isotopes
-ident.hires.noise.HPC <- function(peaklist, allAdducts, scanmode="Negative", look4=c("Cl", "Ac"), identlist=NULL, resol=280000, slope=0, incpt=0, ppm.fixed=10, ppm.iso.fixed=1) {
+ident.hires.noise.HPC <- function(peaklist, allAdducts, scanmode="Negative", look4=c("Cl", "Ac"), identlist=NULL, resol=280000, slope=0, incpt=0, ppm.fixed, ppm.iso.fixed) {
 
   #     peaklist=outlist_Neg
   #     scanmode="Negative"

diff --git a/pipeline/scripts/AddOnFunctions/replaceZeros.R b/pipeline/scripts/AddOnFunctions/replaceZeros.R
@@ -1,4 +1,4 @@
-replaceZeros <- function(file,scanmode,resol,outdir,thresh,scriptDir){
+replaceZeros <- function(file,scanmode,resol,outdir,thresh,scriptDir,ppm){
   # file="./results/grouping_rest/negative_1.RData"
   # scanmode= "negative"
   # scriptDir="./scripts"
@@ -91,7 +91,7 @@ replaceZeros <- function(file,scanmode,resol,outdir,thresh,scriptDir){
   noise.MZ <- noise.MZ[(noise.MZ[ , label] != 0), 1:4]
 
   # Replace "Negative" by "negative" in ident.hires.noise
-  final.outlist.idpat2 = ident.hires.noise.HPC(outpgrlist, allAdducts, scanmode=label2, noise.MZ, look4=look4.add2, resol=resol, slope=0, incpt=0, ppm.fixed=2, ppm.iso.fixed=2)
+  final.outlist.idpat2 = ident.hires.noise.HPC(outpgrlist, allAdducts, scanmode=label2, noise.MZ, look4=look4.add2, resol=resol, slope=0, incpt=0, ppm.fixed=ppm, ppm.iso.fixed=ppm)
   # message(paste(sum(final.outlist.idpat2[ , "assi"] != ""), "assigned noise peaks"))
   tmp <- final.outlist.idpat2[ , c("assi", "theormz")]
   colnames(tmp) <- c("assi_noise",  "theormz_noise")