Adaptations for edge cases and analyses with different genomes, removed the...

Adaptations for edge cases and analyses with different genomes, removed the node constraint in the cluster profile

example/input/cluster.json

@@ -4,7 +4,6 @@
     "group": "zaugg",
     "nCPUs": "{threads}",
     "memory": 20000,
-    "nodes": "avx2",
     "maxTime": "10:00:00",
     "qos": "normal",
     "name": "atac.{rule}.",

example/templates/cluster.json

@@ -4,7 +4,6 @@
     "group": "zaugg",
     "nCPUs": "{threads}",
     "memory": 20000,
-    "nodes": "avx2",
     "maxTime": "10:00:00",
     "qos": "normal",
     "name": "atac.{rule}.",

src/R/aut_stats.R

@@ -167,7 +167,7 @@ for (sampleNameCur in par.l$sampleNames) {
       # 3. Alignment
       #Ivan bug remove * from .s.bam$
 	  file = createFileList(dir_alignment, paste0(sampleNameCur, par.l$fileStatsPattern), verbose = FALSE)
-	  file = file[which(grepl(paste0("/", sampleNameCur,".s.bam"), file))]
+	  file = file[which(grepl(paste0("/", sampleNameCur,".s.bam"), file, fixed = TRUE))]
 	  
 	  if (length(file) != 1) {
 	      errorMessage = paste0("Expecting one file in the directory ", dir_alignment, " that matches ", paste0(sampleNameCur, par.l$fileStatsPattern), ", but found ", length(file), ":", paste0(file, collapse = ","))
@@ -185,7 +185,7 @@ for (sampleNameCur in par.l$sampleNames) {
       
       # 4. Postalign chrM
       file = createFileList(dir_rmChrm, filePatternSample, verbose = FALSE)
-      file = file[which(grepl(paste0("/", sampleNameCur,".", par.l$fileStatsPattern), file))]
+      file = file[which(grepl(paste0("/", sampleNameCur, ".cleaned"), file, fixed = TRUE))]
       if (length(file) != 1) {
           errorMessage = paste0("Expecting one file in the directory ", dir_rmChrm, " that matches ", filePatternSample, ", but found ", length(file), ":", paste0(file, collapse = ","))
           flog.error(errorMessage)
@@ -198,7 +198,7 @@ for (sampleNameCur in par.l$sampleNames) {
       
       # 5. Postalign rmDUP
       file = createFileList(dir_rmDUP, filePatternSample, verbose = FALSE)
-      file = file[which(grepl(paste0("/", sampleNameCur,".", par.l$fileStatsPattern), file))]
+      file = file[which(grepl(paste0("/", sampleNameCur,".cleaned"), file, fixed = TRUE))]
       
       if (length(file) != 1) {
           errorMessage = paste0("Expecting one file in the directory ", dir_rmDUP, " that matches ", filePatternSample, ", but found ", length(file), ":", paste0(file, collapse = ","))
@@ -213,7 +213,7 @@ for (sampleNameCur in par.l$sampleNames) {
     
       # 6. Postalign MAPQ
       file = createFileList(dir_MAPQ, filePatternSample, verbose = FALSE)
-      file = file[which(grepl(paste0("/", sampleNameCur,".", par.l$fileStatsPattern), file))]
+      file = file[which(grepl(paste0("/", sampleNameCur,".cleaned"), file, fixed = TRUE))]
       if (length(file) != 1) {
           errorMessage = paste0("Expecting one file in the directory ", dir_MAPQ, " that matches ", filePatternSample, ", but found ", length(file), ":", paste0(file, collapse = ","))
           flog.error(errorMessage)
@@ -233,7 +233,7 @@ for (sampleNameCur in par.l$sampleNames) {
       
       # Library complexity numbers
       file = createFileList(dir_rmDUP, paste0(sampleNameCur,".", par.l$fileLibraryStatsPattern), verbose = FALSE)
-      file = file[which(grepl(paste0("/", sampleNameCur,".", par.l$fileLibraryStatsPattern), file))]
+      file = file[which(grepl(paste0("/", sampleNameCur,".cleaned"), file, fixed = TRUE))]
       if (length(file) != 1) {
           errorMessage = paste0("Expecting one file in the directory ", dir_rmDUP, " that matches ", paste0(sampleNameCur,".", par.l$fileLibraryStatsPattern), ", but found ", length(file), ":", paste0(file, collapse = ","))
           flog.error(errorMessage)

src/R/functions.R

@@ -182,7 +182,7 @@ checkAndLoadPackages <- function(packages, verbose = TRUE) {
 }
 
 
-.checkAndInstallMissingPackages <- function(packages.vec, verbose = TRUE) {
+.checkAndInstallMissingPackages <- function(packages.vec, verbose = FALSE) {
     
     # if (verbose) cat("Trying to automatically install missing packages. If this fails, install them manually...\n")
     
@@ -350,7 +350,7 @@ createFileList <- function(directory, pattern, recursive = FALSE, ignoreCase = F
 
 
 #' @import checkmate
-generateDefaultReadFlags <- function(pairedEndReads = TRUE, verbose = TRUE) {
+generateDefaultReadFlags <- function(pairedEndReads = TRUE, verbose = FALSE) {
     
     assertFlag(pairedEndReads)
     assertFlag(verbose)
@@ -525,7 +525,7 @@ checkEnrichmentOverTSS <- function(annotationFile, BAMFile, thresholdWithinTSS =
     gr = c(gr.list[[1]], gr.list[[2]])
     gr <- sort(sortSeqlevels(gr))
     
-    flog.info(paste0(" Found", length(gr), " TSS"))
+    flog.info(paste0(" Found ", length(gr), " TSS"))
     
     # If TSS are too close together, the calculation might be biased. Count only reads that uniquely map to one TSS 
     # Exclude TSS that are to close together, delete TSs so that all TSS have at least a distance of 2*thresholdWithinTSS

src/Snakefile

@@ -429,6 +429,9 @@ if not (bamCoverage_normOption == "RPKM") and not (bamCoverage_normOption == "CP
 ###################
 ###################
 
+if annotatePeaks and not config["align"]["assemblyVersion"] in ('hg19', 'hg38', 'mm9', 'mm10'):
+    raise AssertionError("Cannot annotate peaks when the genome assembly is neither human nor mouse. Please set annotatePeaks: false in the config file")
+
 
 peakTypeEncode = ["Encode.broad", "Encode.gapped", "Encode.narrow"]
 peakTypesAll   = ["stringent", "nonStringent"] + peakTypeEncode
@@ -1757,7 +1760,9 @@ def getGenomeTypeMacs2(assemblyVersion):
     elif assemblyVersion == "hg19" or assemblyVersion == "hg38":
         genomeType = "hs"
     else:
-        raise NotImplementedError("Genome assembly version " + assemblyVersion + " not yet implemented for the -g parameter in MACS2.")
+        # Use number instead if not directly supported by MACS2
+        genomeType = str(effectiveGenomeSize)
+        # raise NotImplementedError("Genome assembly version " + assemblyVersion + " not yet implemented for the -g parameter in MACS2.")
 
     return genomeType