Extentions2 ks (#1015)

* seacr lenient * fix * fix * test dag * test dag * test dag * test dag * test dag * log * seacr mode * memory * memory * relaxed * osx --------- Co-authored-by: [email protected] <[email protected]>
maxplanck-ie · Jun 6, 2024 · e91c58c · e91c58c
1 parent 7d94a6f
commit e91c58c
Show file tree

Hide file tree

Showing 14 changed files with 144 additions and 33 deletions.
diff --git a/.ci_stuff/test_dag.sh b/.ci_stuff/test_dag.sh
@@ -221,18 +221,18 @@ if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 628 ]; then exit 1 ; fi
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --bigWigType log2ratio .ci_stuff/organism.yaml .ci_stuff/ChIP.sample_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 562 ]; then exit 1 ; fi
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --peakCaller SEACR .ci_stuff/organism.yaml .ci_stuff/ChIP.sample_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 815 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 927 ]; then exit 1 ; fi
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --peakCaller SEACR --useSpikeInForNorm .ci_stuff/spikein_organism.yaml .ci_stuff/ChIP.sample_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1259 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1371 ]; then exit 1 ; fi
 #noInput
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml .ci_stuff/ChIP.sample_noControl_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 403 ]; then exit 1 ; fi
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --peakCaller Genrich .ci_stuff/organism.yaml .ci_stuff/ChIP.sample_noControl_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 349 ]; then exit 1 ; fi
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --peakCaller SEACR .ci_stuff/organism.yaml .ci_stuff/ChIP.sample_noControl_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 525 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 637 ]; then exit 1 ; fi
 WC=`ChIP-seq -d BAM_input --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --peakCaller SEACR --useSpikeInForNorm .ci_stuff/spikein_organism.yaml .ci_stuff/ChIP.sample_noControl_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 752 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 864 ]; then exit 1 ; fi
 # fromBAM
 WC=`ChIP-seq -d outdir --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --fromBAM BAM_input/filtered_bam/ .ci_stuff/organism.yaml .ci_stuff/ChIP.sample_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1187 ]; then exit 1 ; fi
@@ -289,7 +289,7 @@ WC=`ChIP-seq -d outdir --useSpikeInForNorm --sampleSheet .ci_stuff/test_sampleSh
 if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1045 ]; then exit 1 ; fi
 #multiComp and noInput and fromBam
 WC=`ChIP-seq -d outdir --sampleSheet .ci_stuff/test_sampleSheet_multiComp.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --peakCaller SEACR --fromBAM BAM_input/filtered_bam/ .ci_stuff/spikein_organism.yaml .ci_stuff/ChIP.sample_noControl_config.yaml | tee >(cat 1>&2) | grep -v "conda installation\|Conda environment" | sed '/^\s*$/d' | wc -l`
-if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1078 ]; then exit 1 ; fi
+if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1190 ]; then exit 1 ; fi
 
 
 # mRNA-seq

diff --git a/.github/workflows/osx.yml b/.github/workflows/osx.yml
@@ -52,5 +52,6 @@ jobs:
       - name: createEnvsOSX
         run: |
           micromamba activate snakePipes_CI
+          conda config --add subdirs osx-64
           snakePipes createEnvs --autodetectCondaEnvDir --only ${{matrix.envs}}
      
diff --git a/snakePipes/common_functions.py b/snakePipes/common_functions.py
@@ -23,6 +23,7 @@ def set_env_yamls():
     return {'CONDA_SHARED_ENV': 'envs/shared.yaml',
             'CONDA_CREATE_INDEX_ENV': 'envs/createIndices.yaml',
             'CONDA_RNASEQ_ENV': 'envs/rna_seq.yaml',
+            'CONDA_SALMON_ENV': 'envs/salmon.yaml',
             'CONDA_SLEUTH_ENV': 'envs/sleuth.yaml',
             'CONDA_RMATS_ENV': 'envs/rMats.yaml',
             'CONDA_scRNASEQ_ENV': 'envs/sc_rna_seq.yaml',

diff --git a/snakePipes/shared/rules/ChIP_peak_calling.snakefile b/snakePipes/shared/rules/ChIP_peak_calling.snakefile
@@ -248,14 +248,14 @@ rule prep_bedgraph:
         bedtools genomecov -bg -i filtered_bedgraph/{params.sample}.fragments.bed -g {params.genome} > {output}
         """
 
-rule SEACR_peaks:
+rule SEACR_peaks_stringent:
     input:
         chip = "filtered_bedgraph/{chip_sample}.fragments.bedgraph",
         control = lambda wildcards: "filtered_bedgraph/"+get_control(wildcards.chip_sample)+".fragments.bedgraph" if get_control(wildcards.chip_sample)
                  else []
     output:
         "SEACR/{chip_sample}.filtered.stringent.bed"
-    log: "SEACR/logs/{chip_sample}.log"
+    log: "SEACR/logs/{chip_sample}_stringent.log"
     params:
         fdr = lambda wildcards,input: fdr if not input.control else "",
         prefix = os.path.join(outdir,"SEACR/{chip_sample}.filtered"),
@@ -265,7 +265,24 @@ rule SEACR_peaks:
         bash {params.script} {input.chip} {input.control} {params.fdr} "norm" "stringent" {params.prefix} 2>{log}
         """
 
-rule SEACR_peak_qc:
+rule SEACR_peaks_relaxed:
+    input:
+        chip = "filtered_bedgraph/{chip_sample}.fragments.bedgraph",
+        control = lambda wildcards: "filtered_bedgraph/"+get_control(wildcards.chip_sample)+".fragments.bedgraph" if get_control(wildcards.chip_sample)
+                 else []
+    output:
+        "SEACR/{chip_sample}.filtered.relaxed.bed"
+    log: "SEACR/logs/{chip_sample}_relaxed.log"
+    params:
+        fdr = lambda wildcards,input: fdr if not input.control else "",
+        prefix = os.path.join(outdir,"SEACR/{chip_sample}.filtered"),
+        script=os.path.join(maindir, "shared","tools/SEACR-1.3/SEACR_1.3.sh")
+    conda: CONDA_SEACR_ENV
+    shell: """
+        bash {params.script} {input.chip} {input.control} {params.fdr} "norm" "relaxed" {params.prefix} 2>{log}
+        """
+
+rule SEACR_peak_stringent_qc:
     input:
         bam = "filtered_bam/{sample}.filtered.bam",
         peaks = "SEACR/{sample}.filtered.stringent.bed"
@@ -274,7 +291,41 @@ rule SEACR_peak_qc:
     params:
         genome_index = genome_index
     benchmark:
-        "SEACR/.benchmark/SEACR_peak_qc.{sample}.filtered.benchmark"
+        "SEACR/.benchmark/SEACR_peak_stringent_qc.{sample}.filtered.benchmark"
+    conda: CONDA_SHARED_ENV
+    shell: """
+        # get the number of peaks
+        peak_count=`wc -l < {input.peaks}`
+        
+        # get the number of mapped reads
+        mapped_reads=`samtools view -c -F 4 {input.bam}`
+        
+        #calculate the number of alignments overlapping the peaks
+        # exclude reads flagged as unmapped (unmapped reads will be reported when using -L)
+        reads_in_peaks=`samtools view -c -F 4 -L {input.peaks} {input.bam}`
+        
+        # calculate Fraction of Reads In Peaks
+        frip=`bc -l <<< "$reads_in_peaks/$mapped_reads"`
+        # compute peak genome coverage
+        peak_len=`awk '{{total+=$3-$2}}END{{print total}}' {input.peaks}`
+        genome_size=`awk '{{total+=$3-$2}}END{{print total}}' {params.genome_index}`
+        genomecov=`bc -l <<< "$peak_len/$genome_size"`
+        
+        # write peak-based QC metrics to output file
+        printf "peak_count\tFRiP\tpeak_genome_coverage\n%d\t%5.3f\t%6.4f\n" $peak_count $frip $genomecov > {output.qc}
+        """
+
+
+rule SEACR_peak_relaxed_qc:
+    input:
+        bam = "filtered_bam/{sample}.filtered.bam",
+        peaks = "SEACR/{sample}.filtered.relaxed.bed"
+    output:
+        qc = "SEACR/{sample}.filtered.relaxed_peaks.qc.txt"
+    params:
+        genome_index = genome_index
+    benchmark:
+        "SEACR/.benchmark/SEACR_peak_relaxed_qc.{sample}.filtered.benchmark"
     conda: CONDA_SHARED_ENV
     shell: """
         # get the number of peaks

diff --git a/snakePipes/shared/rules/ChIP_peak_calling_spikein.snakefile b/snakePipes/shared/rules/ChIP_peak_calling_spikein.snakefile
@@ -248,14 +248,14 @@ rule prep_bedgraph:
         bigWigToBedGraph {input} {output}
         """
 
-rule SEACR_peaks:
+rule SEACR_peaks_stringent:
     input:
         chip = "filtered_bedgraph/{chip_sample}_host.fragments.bedgraph",
         control = lambda wildcards: "filtered_bedgraph/"+get_control(wildcards.chip_sample)+"_host.fragments.bedgraph" if get_control(wildcards.chip_sample)
                  else []
     output:
         "SEACR/{chip_sample}_host.stringent.bed"
-    log: "SEACR/logs/{chip_sample}.log"
+    log: "SEACR/logs/{chip_sample}_stringent.log"
     params:
         fdr = lambda wildcards,input: fdr if not input.control else "",
         prefix = os.path.join(outdir,"SEACR/{chip_sample}_host"),
@@ -265,8 +265,25 @@ rule SEACR_peaks:
         bash {params.script} {input.chip} {input.control} {params.fdr} "non" "stringent" {params.prefix} 2>{log}
         """
 
+rule SEACR_peaks_relaxed:
+    input:
+        chip = "filtered_bedgraph/{chip_sample}_host.fragments.bedgraph",
+        control = lambda wildcards: "filtered_bedgraph/"+get_control(wildcards.chip_sample)+"_host.fragments.bedgraph" if get_control(wildcards.chip_sample)
+                 else []
+    output:
+        "SEACR/{chip_sample}_host.relaxed.bed"
+    log: "SEACR/logs/{chip_sample}_relaxed.log"
+    params:
+        fdr = lambda wildcards,input: fdr if not input.control else "",
+        prefix = os.path.join(outdir,"SEACR/{chip_sample}_host"),
+        script=os.path.join(maindir, "shared","tools/SEACR-1.3/SEACR_1.3.sh")
+    conda: CONDA_SEACR_ENV
+    shell: """
+        bash {params.script} {input.chip} {input.control} {params.fdr} "non" "relaxed" {params.prefix} 2>{log}
+        """
+
 
-rule SEACR_peak_qc:
+rule SEACR_peak_stringent_qc:
     input:
         bam = "split_bam/{sample}_host.bam",
         peaks = "SEACR/{sample}_host.stringent.bed"
@@ -275,7 +292,7 @@ rule SEACR_peak_qc:
     params:
         genome_index = genome_index
     benchmark:
-        "SEACR/.benchmark/SEACR_peak_qc.{sample}_host.benchmark"
+        "SEACR/.benchmark/SEACR_peak_qc.{sample}_host_stringend.benchmark"
     conda: CONDA_SHARED_ENV
     shell: """
         # get the number of peaks
@@ -299,4 +316,35 @@ rule SEACR_peak_qc:
         printf "peak_count\tFRiP\tpeak_genome_coverage\n%d\t%5.3f\t%6.4f\n" $peak_count $frip $genomecov > {output.qc}
         """
 
-
+rule SEACR_peak_relaxed_qc:
+    input:
+        bam = "split_bam/{sample}_host.bam",
+        peaks = "SEACR/{sample}_host.relaxed.bed"
+    output:
+        qc = "SEACR/{sample}_host.relaxed_peaks.qc.txt"
+    params:
+        genome_index = genome_index
+    benchmark:
+        "SEACR/.benchmark/SEACR_peak_qc.{sample}_host_relaxed.benchmark"
+    conda: CONDA_SHARED_ENV
+    shell: """
+        # get the number of peaks
+        peak_count=`wc -l < {input.peaks}`
+        
+        #get the number of mapped reads
+        mapped_reads=`samtools view -c -F 4 {input.bam}`
+        
+        #calculate the number of alignments overlapping the peaks
+        # exclude reads flagged as unmapped (unmapped reads will be reported when using -L)
+        reads_in_peaks=`samtools view -c -F 4 -L {input.peaks} {input.bam}`
+        
+        # calculate Fraction of Reads In Peaks
+        frip=`bc -l <<< "$reads_in_peaks/$mapped_reads"`
+        # compute peak genome coverage
+        peak_len=`awk '{{total+=$3-$2}}END{{print total}}' {input.peaks}`
+        genome_size=`awk '{{total+=$3-$2}}END{{print total}}' {params.genome_index}`
+        genomecov=`bc -l <<< "$peak_len/$genome_size"`
+        
+        # write peak-based QC metrics to output file
+        printf "peak_count\tFRiP\tpeak_genome_coverage\n%d\t%5.3f\t%6.4f\n" $peak_count $frip $genomecov > {output.qc}
+        """
diff --git a/snakePipes/shared/rules/Salmon.snakefile b/snakePipes/shared/rules/Salmon.snakefile
@@ -80,7 +80,7 @@ else:
             gtf = genes_gtf,
             lib_type = getSalmon_libtype(pairedEnd, libraryType)
         threads: 8
-        conda: CONDA_RNASEQ_ENV
+        conda: CONDA_SALMON_ENV
         shell: """
             salmon quant -p {threads} --softclipOverhangs --validateMappings --numBootstraps 50 -i {params.index} -l {params.lib_type} -r {input.fastq} -o {params.outdir} 2> {log}
             """

diff --git a/snakePipes/shared/rules/Salmon_allelic.snakefile b/snakePipes/shared/rules/Salmon_allelic.snakefile
@@ -79,7 +79,7 @@ samtools fastq -1 {output.r1} -2 {output.r2} -0 /dev/null -s /dev/null -n  \\
             gtf = genes_gtf,
             lib_type = getSalmon_libtype(pairedEnd, libraryType)
         threads: 8
-        conda: CONDA_RNASEQ_ENV
+        conda: CONDA_SALMON_ENV
         shell: """
             salmon quant -p {threads} --softclipOverhangs --validateMappings --numBootstraps 50 -i {params.index} -l {params.lib_type} -1 {input.r1} -2 {input.r2} -o {params.outdir} 2> {log}
             """
@@ -114,7 +114,7 @@ samtools fastq -1 /dev/null -2 /dev/null -0 /dev/null -s {output.r1} -n  \\
             gtf = genes_gtf,
             lib_type = getSalmon_libtype(pairedEnd, libraryType)
         threads: 8
-        conda: CONDA_RNASEQ_ENV
+        conda: CONDA_SALMON_ENV
         shell: """
             salmon quant -p {threads} --softclipOverhangs --validateMappings --numBootstraps 50 -i {params.index} -l {params.lib_type} -r {input.fastq} -o {params.outdir} 2> {log}
             """

diff --git a/snakePipes/shared/rules/createIndices.snakefile b/snakePipes/shared/rules/createIndices.snakefile
@@ -411,7 +411,7 @@ rule Salmon_index_joint_fa:
         err = "SalmonIndex_RNAVelocity/logs/SalmonIndex.err",
         out = "SalmonIndex_RNAVelocity/logs/SalmonIndex.out"
     threads: lambda wildcards: 16 if 16<max_thread else max_thread
-    conda: CONDA_RNASEQ_ENV
+    conda: CONDA_SALMON_ENV
     shell:"""
         cat {input.joint_fasta} {input.genome_fasta} > {output.seq_fa}
         salmon index -p {threads} -t {output.seq_fa} -d {input.decoys} -i SalmonIndex_RNAVelocity {params.salmonIndexOptions} > {log.out} 2> {log.err}

diff --git a/snakePipes/shared/rules/envs/rna_seq.yaml b/snakePipes/shared/rules/envs/rna_seq.yaml
@@ -1,4 +1,4 @@
-name: snakepipes_RNAseq_environment_2.0
+name: snakepipes_RNAseq_environment_3.0
 channels:
  - conda-forge
  - bioconda
@@ -9,8 +9,8 @@ dependencies:
  - subread = 2.0.1
  - hisat2 = 2.2.1
  - star = 2.7.10b
- - salmon = 1.10.1
- - r-base = 4.1.3
+# - salmon = 1.10.1
+ - r-base = 4.3.0
  - r-wasabi
  - r-dplyr
  - r-ggplot2
@@ -20,7 +20,7 @@ dependencies:
  - pandoc
  - r-dt
  - r-rcolorbrewer
- - bioconductor-deseq2
+ - bioconductor-deseq2 = 1.42.0
  - bioconductor-tximport
  - r-knitcitations
  - bioconductor-apeglm

diff --git a/snakePipes/shared/rules/envs/salmon.yaml b/snakePipes/shared/rules/envs/salmon.yaml
@@ -0,0 +1,6 @@
+name: snakepipes_salmon_environment_1.0
+channels:
+ - conda-forge
+ - bioconda
+dependencies:
+ - salmon = 1.10.1
diff --git a/snakePipes/shared/rules/scRNAseq_Alevin.snakefile b/snakePipes/shared/rules/scRNAseq_Alevin.snakefile
@@ -39,7 +39,7 @@ rule SalmonAlevin:
             out =  "Alevin/logs/alevin.{sample}.out",
             err = "Alevin/logs/alevin.{sample}.err"
         #Use RNAseq env because Salmon already installed there (no need for duplication).
-        conda: CONDA_RNASEQ_ENV
+        conda: CONDA_SALMON_ENV
         threads: 8
         shell:"""
             salmon alevin -l {params.libtype} -1 {input.R1} -2 {input.R2} {params.protocol} -i {params.index} -p {threads} -o {params.outdir} --tgMap {input.tgMap} --dumpFeatures --dumpMtx --numCellBootstraps 100 > {log.out} 2> {log.err}
@@ -124,7 +124,7 @@ rule AlevinForVelocity:
             out =  "AlevinForVelocity/logs/alevin.{sample}.out",
             err = "AlevinForVelocity/logs/alevin.{sample}.err"
         #Use RNAseq env because Salmon already installed there (no need for duplication).
-        conda: CONDA_RNASEQ_ENV
+        conda: CONDA_SALMON_ENV
         threads: 8
         shell:"""
             salmon alevin -l {params.libtype} -1 {input.R1} -2 {input.R2} {params.protocol} -i {params.velo_index} -p {threads} -o {params.outdir} --tgMap {params.tgMap} --dumpFeatures --dumpMtx --numCellBootstraps 100 > {log.out} 2> {log.err}

diff --git a/snakePipes/workflows/ChIP-seq/Snakefile b/snakePipes/workflows/ChIP-seq/Snakefile
@@ -209,11 +209,11 @@ def run_macs2():
 def run_seacr():
     if (peakCaller == "SEACR"):
         if useSpikeInForNorm:
-            file_list = expand("SEACR/{chip_sample}_host.stringent.bed", chip_sample = chip_samples)
-            file_list.append(expand("SEACR/{chip_sample}_host.stringent_peaks.qc.txt",chip_sample=chip_samples))
+            file_list = expand("SEACR/{chip_sample}_host.{mode}.bed", chip_sample = chip_samples,mode=["stringent","relaxed"])
+            file_list.append(expand("SEACR/{chip_sample}_host.{mode}_peaks.qc.txt",chip_sample=chip_samples,mode=["stringent","relaxed"]))
         else:
-            file_list = expand("SEACR/{chip_sample}.filtered.stringent.bed", chip_sample = chip_samples)
-            file_list.append(expand("SEACR/{chip_sample}.filtered.stringent_peaks.qc.txt",chip_sample=chip_samples))
+            file_list = expand("SEACR/{chip_sample}.filtered.{mode}.bed", chip_sample = chip_samples,mode=["stringent","relaxed"])
+            file_list.append(expand("SEACR/{chip_sample}.filtered.{mode}_peaks.qc.txt",chip_sample=chip_samples,mode=["stringent","relaxed"]))
         return (file_list)
     return ([])
 

diff --git a/snakePipes/workflows/ChIP-seq/cluster.yaml b/snakePipes/workflows/ChIP-seq/cluster.yaml
@@ -8,3 +8,7 @@ Genrich_peaks:
     memory: 20G
 namesort_bams:
     memory: 6G
+SEACR_peaks_stringent:
+    memory: 10G
+SEACR_peaks_lenient:
+    memory: 20G
diff --git a/tests/test_jobcounts.py b/tests/test_jobcounts.py
@@ -589,7 +589,7 @@ def test_seacr(self, ifs):
         print(' '.join([str(i) for i in ci]))
         _p = sp.run(ci, capture_output=True, text=True)
         assert _p.returncode == 0
-        assert parseSpOut(_p) == 77
+        assert parseSpOut(_p) == 89
     def test_seacr_spikein(self, ifs):
         ci = [
             "ChIP-seq",
@@ -608,7 +608,7 @@ def test_seacr_spikein(self, ifs):
         print(' '.join([str(i) for i in ci]))
         _p = sp.run(ci, capture_output=True, text=True)
         assert _p.returncode == 0
-        assert parseSpOut(_p) == 118
+        assert parseSpOut(_p) == 130
     def test_SE(self, ifs):
         ci = [
             "ChIP-seq",
@@ -695,7 +695,7 @@ def test_seacr_noInput(self, ifs):
         print(' '.join([str(i) for i in ci]))
         _p = sp.run(ci, capture_output=True, text=True)
         assert _p.returncode == 0
-        assert parseSpOut(_p) == 50
+        assert parseSpOut(_p) == 62
     def test_seacr_spikein_noInput(self, ifs):
         ci = [
             "ChIP-seq",
@@ -714,7 +714,7 @@ def test_seacr_spikein_noInput(self, ifs):
         print(' '.join([str(i) for i in ci]))
         _p = sp.run(ci, capture_output=True, text=True)
         assert _p.returncode == 0
-        assert parseSpOut(_p) == 71
+        assert parseSpOut(_p) == 83
     def test_frombam(self, ifs):
         ci = [
             "ChIP-seq",
@@ -1121,7 +1121,7 @@ def test_multicomp_fromBam_noInput_SEACR(self, ifs):
         print(' '.join([str(i) for i in ci]))
         _p = sp.run(ci, capture_output=True, text=True)
         assert _p.returncode == 0
-        assert parseSpOut(_p) == 108
+        assert parseSpOut(_p) == 120
 
 class TestmRNAseq:
     def test_default(self, ifs):