diff --git a/src/R/annotatePeaks.R b/src/R/annotatePeaks.R
index abbafbb7a82662ca9ad86f38dfc4e4f6b73110aa..692fdc279128621bef2aa142ac1a838a357c2e8f 100755
--- a/src/R/annotatePeaks.R
+++ b/src/R/annotatePeaks.R
@@ -99,11 +99,11 @@ for (peakType in unique(names(par.l$input_peaks))){
                 x = par.l$output[list.name][[1]], value = T)
     
     # The first can happen when the file is empty
-    if (length(readLines(file, n=1)) == 0 || readLines(file, n=1) == "# No consensus peaks found" ) {
+    if (length(readLines(file, n=1)) == 0 || readLines(file, n=1) == "# No consensus peaks found" || length(readLines(file, n=10)) < 10 ) {
         
         pdf(outIMG)
         plot(c(0, 1), c(0, 1), ann = FALSE, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
-        message = paste0("No consensus peaks found")
+        message = paste0("No or not enough (<10) consensus peaks found")
         text(x = 0.5, y = 0.5, message, cex = 1.6, col = "red")
         dev.off()
         
@@ -129,8 +129,28 @@ for (peakType in unique(names(par.l$input_peaks))){
     
     
     pdf(outIMG)
-    plotAnnoPie(peaks.annotated, main = "Peak Annotation")
-    plot(plotDistToTSS(peaks.annotated))
+    
+    res = tryCatch({ 
+        plotAnnoPie(peaks.annotated, main = "Peak Annotation")
+    }, error = function(cond) {
+        message(paste("Could not plot DistToTSS plot"))
+        # Choose a return value in case of error
+        return(NA)
+    }
+    )
+    
+    
+    
+    res = tryCatch({ 
+        plot(plotDistToTSS(peaks.annotated))
+    }, error = function(cond) {
+        message(paste("Could not plot DistToTSS plot"))
+        # Choose a return value in case of error
+        return(NA)
+    }
+    )
+    
+    
     dev.off()
     
     if (par.l$tabulateOutput){
diff --git a/src/Snakefile b/src/Snakefile
index a1f9483bab0a514d65571bb2764aa0c0d82b9992..571296fd97d97c4874bd6608077de1271ba560dd 100755
--- a/src/Snakefile
+++ b/src/Snakefile
@@ -1764,7 +1764,7 @@ def getGenomeTypeMacs2(assemblyVersion):
 macs2_stringent_outputName     = PEAK_STRINGENT_dir  + '/{basename}' + '.stringent'
 #macs2_stringent_outputName     = expand("{dir}/{sample}.{{specifics}}.stringent", sample = [allSamplesUnique, allIndividualsUnique] if mergeReplicates else allSamplesUnique, dir = PEAK_STRINGENT_dir)
 
-def getControlSample(filePrefix):
+def getControlSample(filePrefix, addCommandLinePrefix):
 #
     # For ATACSeq, no control sample
     if dataType != "ChIPseq":
@@ -1785,18 +1785,22 @@ def getControlSample(filePrefix):
     controlSample = samplesData["control_sampleName"][samplesData[sampleTableCol]==sample_cur].iloc[0]
 
     #print(controlSample)
+    prefix = ""
+    if addCommandLinePrefix:
+        prefix = "-c "
 
     if pandas.isnull(controlSample):
         ctrl_param = ""
     else:
-        ctrl_param = "-c " + FINAL_OUTPUT_dir + "/" + controlSample + ".final.bam"
+        ctrl_param = prefix + FINAL_OUTPUT_dir + "/" + controlSample + ".final.bam"
     return ctrl_param
 
 
 # Runs for all bam files in the final output folder, also the replicated one
 rule macs2_stringent:
     input:
-        bam = FINAL_OUTPUT_dir + "/{basename}.bam"
+        bam   = FINAL_OUTPUT_dir + "/{basename}.bam",
+        index = lambda wildcards: [getControlSample(wildcards.basename, False)]
     output:
         peaks_bedT   = temp(macs2_stringent_outputName + '_peaks.narrowPeak'),
         peaks_bed    = temp(macs2_stringent_outputName + '.narrowPeak.gz'),
@@ -1819,7 +1823,7 @@ rule macs2_stringent:
         outputDir       = PEAK_STRINGENT_dir,
         pairedEndString = macs2_pairedEndMode,
         keepDuplicates  = "--keep-dup all",
-        controlFile = lambda wildcards: getControlSample(wildcards.basename)
+        controlFile = lambda wildcards: getControlSample(wildcards.basename, True)
     # group: "peaks"
     shell:
         """PYTHONPATH="" &&
@@ -1865,7 +1869,7 @@ use rule macs2_stringent as macs2_nonStringent with:
         outputDir       = PEAK_NONSTRINGENT_dir,
         keepDuplicates  = "--keep-dup all",
         pairedEndString = macs2_pairedEndMode,
-        controlFile = lambda wildcards: getControlSample(wildcards.basename)
+        controlFile = lambda wildcards: getControlSample(wildcards.basename, True)
 
 
 # Define the general output file name above the rule. Wildcards are resolved in the resulting string when evaluating the rule
@@ -1873,7 +1877,8 @@ macs2_Encode_outputName     = PEAK_ENCODE_dir + '/{basename}' + '.Encode'
 
 rule macs2_Encode:
     input:
-        bam = expand('{dir}/{{basename}}.bam', dir = FINAL_OUTPUT_dir)
+        bam = expand('{dir}/{{basename}}.bam', dir = FINAL_OUTPUT_dir),
+        index = lambda wildcards: [getControlSample(wildcards.basename, False)]
     output:
         broadPeakfileT  = temp(macs2_Encode_outputName + '_peaks.broadPeak'),
         gappedPeakfileT = temp(macs2_Encode_outputName + '_peaks.gappedPeak'),
@@ -1905,7 +1910,7 @@ rule macs2_Encode:
         keepDuplicates         = "--keep-dup all",
         bdg1                   = macs2_Encode_outputName + '_control_lambda.bdg',
         bdg2                   = macs2_Encode_outputName + '_treat_pileup.bdg',
-        controlFile            = lambda wildcards: getControlSample(wildcards.basename)
+        controlFile            = lambda wildcards: getControlSample(wildcards.basename, True)
 
     shell:
         # 1. First produce broad and gapped peaks, then narrow ones