search

kundajelab / atac_dnase_pipelines

ATAC-seq and DNase-seq processing pipeline
BSD 3-Clause "New" or "Revised" License
161 stars 81 forks source link

Unknown Error in IDR Step #53

Closed smuchnik closed 7 years ago

smuchnik commented 7 years ago

Hello, I am attempting to run your ATAC-Seq pipeline, just on some sample ENCODE data at the moment (thank you for sharing it!). I'm consistently running into an error at the IDR step that I'm not sure how to interpret and would appreciate any insight you have into what is going wrong.

Thank you!

Full log: https://gist.github.com/smuchnik/f496b3a880fed9a06cbf51e6904f4af2

Error message:

Task failed: Program & line : '/home/sm2576/project/ATAC/atac_dnase_pipelines/modules/callpeak_idr.bds', line 74 Task Name : 'idr2 rep1-pr' Task ID : 'atac.bds.20170605_132057_849/task.callpeak_idr.idr2_rep1_pr.line_74.id_43' Task PID : '8589' Task hint : 'idr --samples /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/pseudo_reps/rep1/pr1/Test_Read1.trim.PE2SE.nodup.pr1.tn5.pf.pval0.1.500K.narrowPeak.gz /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/pseudo_reps/rep1/pr2/Test_Read1.trim.PE2SE.nodup.pr2.tn5.pf' Task resources : 'cpus: -1 mem: -1.0 B wall-timeout: 8640000' State : 'ERROR' Dependency state : 'ERROR' Retries available : '1' Input files : '[/ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/pseudo_reps/rep1/pr1/Test_Read1.trim.PE2SE.nodup.pr1.tn5.pf.pval0.1.500K.narrowPeak.gz, /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/pseudo_reps/rep1/pr2/Test_Read1.trim.PE2SE.nodup.pr2.tn5.pf.pval0.1.500K.narrowPeak.gz, /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/rep1/Test_Read1.trim.PE2SE.nodup.tn5.pf.pval0.1.500K.narrowPeak.gz]' Output files : '[/ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.filt.narrowPeak.gz, /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.unthresholded-peaks.txt.png, /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.log.txt, /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.unthresholded-peaks.txt.gz, /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.filt.12-col.bed.gz]' Script file : '/ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/atac.bds.20170605_132057_849/task.callpeak_idr.idr2_rep1_pr.line_74.id_43.sh' Exit status : '1' Program : 

    # SYS command. line 76

     if [[ -f $(which conda) && $(conda env list | grep bds_atac_py3 | wc -l) != "0" ]]; then source activate bds_atac_py3; sleep 5; fi;  export PATH=/home/sm2576/project/ATAC/atac_dnase_pipelines/.:/home/sm2576/project/ATAC/atac_dnase_pipelines/modules:/home/sm2576/project/ATAC/atac_dnase_pipelines/utils:${PATH}:/bin:/usr/bin:/usr/local/bin:${HOME}/.bds; set -o pipefail; STARTTIME=$(date +%s)

    # SYS command. line 78

     idr --samples /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/pseudo_reps/rep1/pr1/Test_Read1.trim.PE2SE.nodup.pr1.tn5.pf.pval0.1.500K.narrowPeak.gz /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/pseudo_reps/rep1/pr2/Test_Read1.trim.PE2SE.nodup.pr2.tn5.pf.pval0.1.500K.narrowPeak.gz --peak-list /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/rep1/Test_Read1.trim.PE2SE.nodup.tn5.pf.pval0.1.500K.narrowPeak.gz --input-file-type narrowPeak \
                --output-file /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.unthresholded-peaks.txt --rank p.value --soft-idr-threshold 0.1 \
                --plot --use-best-multisummit-IDR --log-output-file /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.log.txt

    # SYS command. line 82

     idr_thresh_transformed=$(awk -v p=0.1 'BEGIN{print -log(p)/log(10)}')

    # SYS command. line 85

     awk 'BEGIN{OFS="\t"} $12>='"${idr_thresh_transformed}"' {if ($2<0) $2=0; print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12,"0"}' /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.unthresholded-peaks.txt \
                | sort | uniq | sort -k7n,7n | gzip -nc > /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.13-col.bed.gz

    # SYS command. line 88

     zcat /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.13-col.bed.gz | awk 'BEGIN{OFS="\t"} {print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10}' | gzip -nc > /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.narrowPeak.gz

    # SYS command. line 89

     zcat /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.13-col.bed.gz | awk 'BEGIN{OFS="\t"} {print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12}' | gzip -nc > /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.12-col.bed.gz

    # SYS command. line 91

     bedtools intersect -v -a <(zcat -f /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.13-col.bed.gz) -b <(zcat -f /ycga-gpfs/project/fas/sestan/sm2576/ATAC/genome_data/hg38/hg38.blacklist.bed.gz) | grep -P 'chr[\dXY]+[ \t]' | awk 'BEGIN{OFS="\t"} {if ($5>1000) $5=1000; print $0}' | gzip -nc > /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.filt.13-col.bed.gz

    # SYS command. line 92

     zcat /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.filt.13-col.bed.gz | awk 'BEGIN{OFS="\t"} {print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10}' | gzip -nc > /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.filt.narrowPeak.gz

    # SYS command. line 93

     zcat /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.filt.13-col.bed.gz | awk 'BEGIN{OFS="\t"} {print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12}' | gzip -nc > /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.filt.12-col.bed.gz

    # SYS command. line 95

     gzip -f /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.unthresholded-peaks.txt

    # SYS command. line 96

     rm -f /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.13-col.bed.gz /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.IDR0.1.filt.13-col.bed.gz

    # SYS command. line 98

     TASKTIME=$[$(date +%s)-${STARTTIME}]; echo "Task has finished (${TASKTIME} seconds)."

Fatal error: /home/sm2576/project/ATAC/atac_dnase_pipelines/atac.bds, line 1277, pos 2. Task/s failed.

Creating checkpoint file: Config or command line option disabled checkpoint file creation, nothing done.

leepc12 commented 7 years ago

Can I see files in /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1?

$ ls -l /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1
smuchnik commented 7 years ago

Sure, that folder contains the following:

total 256 -rw-r--r-- 1 sm2576 sestan 20 Jun 5 13:55 screen_test5_rep1-pr.IDR0.1.12-col.bed.gz -rw-r--r-- 1 sm2576 sestan 20 Jun 5 13:55 screen_test5_rep1-pr.IDR0.1.13-col.bed.gz -rw-r--r-- 1 sm2576 sestan 20 Jun 5 13:55 screen_test5_rep1-pr.IDR0.1.filt.13-col.bed.gz -rw-r--r-- 1 sm2576 sestan 20 Jun 5 13:55 screen_test5_rep1-pr.IDR0.1.narrowPeak.gz -rw-r--r-- 1 sm2576 sestan 12642 Jun 5 13:55 screen_test5_rep1-pr.unthresholded-peaks.txt -rw-r--r-- 1 sm2576 sestan 88898 Jun 5 13:55 screen_test5_rep1-pr.unthresholded-peaks.txt.noalternatesummitpeaks.png

leepc12 commented 7 years ago

It looks like none of the peaks passed IDR threshold of 0.1. Can you check the following? Can you also post or email me /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.unthresholded-peaks.txt?

$ idr_thresh_transformed=$(awk -v p=0.1 'BEGIN{print -log(p)/log(10)}')

$ awk 'BEGIN{OFS="\t"} $12>='"${idr_thresh_transformed}"' {if ($2<0) $2=0; print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12,"0"}' \
    /ycga-gpfs/project/fas/sestan/sm2576/ATAC/screen_test5/out/peak/macs2/idr/pseudo_reps/rep1/screen_test5_rep1-pr.unthresholded-peaks.txt \
    | sort | uniq | sort -k7n,7n | wc -l
smuchnik commented 7 years ago

Ah I see, that makes sense. Thank you!