maxplanck-ie / snakepipes

Customizable workflows based on snakemake and python for the analysis of NGS data
http://snakepipes.readthedocs.io
381 stars 85 forks source link

add peak qc for SEACR as for MACS2 #998

Closed katsikora closed 5 months ago

sunta3iouxos commented 6 months ago

maybe this is related: I go the following error that never had using the pipeline and doing peak calling. As a note I am running the developmental version

Error in rule MACS2_peak_qc:
    jobid: 33
    input: filtered_bam/A006200387_220458_S15.filtered.bam, MACS2/A006200387_220458_S15.filtered.BAM_peaks.xls
    output: MACS2/A006200387_220458_S15.filtered.BAM_peaks.qc.txt
    conda-env: /home/tgeorgom/CHN01/bam/None/264d54d86c09dd2633f7597f5b77ab09_
    shell:

        # get the number of peaks
        peak_count=`wc -l < MACS2/A006200387_220458_S15.filtered.BAM_peaks.narrowPeak`

        # get the number of mapped reads
        mapped_reads=`samtools view -c -F 4 filtered_bam/A006200387_220458_S15.filtered.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 MACS2/A006200387_220458_S15.filtered.BAM_peaks.narrowPeak filtered_bam/A006200387_220458_S15.filtered.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}' MACS2/A006200387_220458_S15.filtered.BAM_peaks.narrowPeak`
        genome_size=`awk '{total+=$3-$2}END{print total}' /home/tgeorgom/GRCm38_gencode_release19/genome_fasta/genome.fa.fai`
        genomecov=`bc -l <<< "$peak_len/$genome_size"`

        # write peak-based QC metrics to output file
        printf "peak_count      FRiP    peak_genome_coverage
%d      %5.3f   %6.4f
" $peak_count $frip $genomecov > MACS2/A006200387_220458_S15.filtered.BAM_peaks.qc.txt

        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)

Exiting because a job execution failed. Look above for error message
Complete log: .snakemake/log/2024-04-02T163604.681904.snakemake.log