test/cmdlinetest

   1 #!/bin/bash
   2 # Time-stamp: <2024-02-15 12:30:51 Tao Liu>
   3
   4 # integrative subcmds testing
   5
   6 if [ $# -lt 1 ];then
   7     echo "Run all tests for subcommands of MACS3. Need >=1 parameter for a tag name! A unique string combining date, time and MACS3 version is recommended. ./cmdlinetest <TAG> [UPDATE]"
   8     echo ""
   9     echo "If UPDATE is set with anything, aka, this command sees a second parameter, the standard output results will be replaced."
  10     exit
  11 fi
  12
  13 TAG=$1
  14 if [ "$#" = "2" ];then
  15     REPLACE="Y"
  16 else
  17     REPLACE="N"
  18 fi
  19
  20 TEMPDIR=../temp
  21 OUTPUTDIR_PREFIX=${TEMPDIR}/${TAG}
  22
  23 # make temp directory if necessary
  24 if [ ! -d $TEMPDIR ]; then mkdir $TEMPDIR; fi
  25
  26 CHIP=CTCF_SE_ChIP_chr22_50k.bed.gz
  27 CTRL=CTCF_SE_CTRL_chr22_50k.bed.gz
  28
  29 CHIPPE=CTCF_PE_ChIP_chr22_50k.bam
  30 CTRLPE=CTCF_PE_CTRL_chr22_50k.bam
  31
  32 CHIPBEDPE=CTCF_PE_ChIP_chr22_50k.bedpe.gz
  33 CTRLBEDPE=CTCF_PE_CTRL_chr22_50k.bedpe.gz
  34
  35 CHIPCONTIGS50K=contigs50k.bed.gz
  36
  37 CHIPBIGSPEEDTEST=CTCF_12878_5M.bed.gz
  38 CTRLBIGSPEEDTEST=input_12878_5M.bed.gz
  39
  40 CALLVARPEAK=callvar_testing.narrowPeak
  41
  42 ATACSEQBAM=yeast_500k_SRR1822137.bam
  43 ATACSEQBED=yeast_500k_SRR1822137.bedpe.gz
  44
  45 # callpeak
  46 echo "1. callpeak"
  47 echo "1.1 callpeak narrow"
  48
  49 mkdir ${OUTPUTDIR_PREFIX}_run_callpeak_narrow
  50
  51 macs3 callpeak -g 52000000 -t $CHIP -c $CTRL -n run_callpeak_narrow0 -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_narrow --cutoff-analysis &> ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0.log
  52 macs3 callpeak -g 52000000 -t $CHIP -c $CTRL -n run_callpeak_narrow1 -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_narrow --d-min 15 --call-summits &> ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow1.log
  53 macs3 callpeak -g 52000000 -t $CHIP -c $CTRL -n run_callpeak_narrow2 -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_narrow --nomodel --extsize 100 &> ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow2.log
  54 macs3 callpeak -g 52000000 -t $CHIP -c $CTRL -n run_callpeak_narrow3 -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_narrow --nomodel --extsize 100 --shift -50 &> ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow3.log
  55 macs3 callpeak -g 52000000 -t $CHIP -c $CTRL -n run_callpeak_narrow4 -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_narrow --nomodel --nolambda --extsize 100 --shift -50 &> ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow4.log
  56 macs3 callpeak -g 52000000 -t $CHIP -c $CTRL -n run_callpeak_narrow5 -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_narrow --scale-to large &> ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow5.log
  57
  58 echo "1.2 callpeak broad"
  59 mkdir ${OUTPUTDIR_PREFIX}_run_callpeak_broad
  60
  61 macs3 callpeak -g 52000000 -t $CHIP -c $CTRL -n run_callpeak_broad -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_broad --broad &> ${OUTPUTDIR_PREFIX}_run_callpeak_broad/run_callpeak_broad.log
  62
  63 echo "1.3 callpeak on PE narrow/broad"
  64
  65 mkdir  ${OUTPUTDIR_PREFIX}_run_callpeak_pe_narrow
  66
  67 macs3 callpeak -g 52000000 -f BAMPE -t $CHIPPE -c $CTRLPE -n run_callpeak_bampe_narrow -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_pe_narrow --call-summits &> ${OUTPUTDIR_PREFIX}_run_callpeak_pe_narrow/run_callpeak_bampe_narrow.log
  68 macs3 callpeak -g 52000000 -f BEDPE -t $CHIPBEDPE -c $CTRLBEDPE -n run_callpeak_bedpe_narrow -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_pe_narrow --call-summits &> ${OUTPUTDIR_PREFIX}_run_callpeak_pe_narrow/run_callpeak_bedpe_narrow.log
  69 macs3 callpeak -g 52000000 -f BEDPE -t $CHIPBEDPE -n run_callpeak_pe_narrow_onlychip -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_pe_narrow &> ${OUTPUTDIR_PREFIX}_run_callpeak_pe_narrow/run_callpeak_pe_narrow_onlychip.log
  70
  71 mkdir  ${OUTPUTDIR_PREFIX}_run_callpeak_pe_broad
  72
  73 macs3 callpeak -g 52000000 -f BAMPE -t $CHIPPE -c $CTRLPE -n run_callpeak_bampe_broad -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_pe_broad --broad &> ${OUTPUTDIR_PREFIX}_run_callpeak_pe_broad/run_callpeak_bampe_broad.log
  74 macs3 callpeak -g 52000000 -f BEDPE -t $CHIPBEDPE -c $CTRLBEDPE -n run_callpeak_bedpe_broad -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_pe_broad --broad &> ${OUTPUTDIR_PREFIX}_run_callpeak_pe_broad/run_callpeak_bedpe_broad.log
  75
  76 # pileup
  77 echo "2. pileup"
  78
  79 mkdir ${OUTPUTDIR_PREFIX}_run_pileup
  80
  81 macs3 pileup -f BED -i $CHIP --extsize 200 --outdir ${OUTPUTDIR_PREFIX}_run_pileup -o run_pileup_ChIP.bed.bdg &> ${OUTPUTDIR_PREFIX}_run_pileup/run_pileup_ChIP.bed.log
  82 macs3 pileup -f BED -i $CTRL --extsize 200 --outdir ${OUTPUTDIR_PREFIX}_run_pileup -o run_pileup_CTRL.bed.bdg &> ${OUTPUTDIR_PREFIX}_run_pileup/run_pileup_CTRL.bed.log
  83
  84 macs3 pileup -f BAMPE -i $CHIPPE --outdir ${OUTPUTDIR_PREFIX}_run_pileup -o run_pileup_ChIPPE.bampe.bdg &> ${OUTPUTDIR_PREFIX}_run_pileup/run_pileup_ChIPPE.bampe.log
  85 macs3 pileup -f BAMPE -i $CTRLPE --outdir ${OUTPUTDIR_PREFIX}_run_pileup -o run_pileup_CTRLPE.bampe.bdg &> ${OUTPUTDIR_PREFIX}_run_pileup/run_pileup_CTRLPE.bampe.log
  86
  87 macs3 pileup -f BEDPE -i $CHIPBEDPE --outdir ${OUTPUTDIR_PREFIX}_run_pileup -o run_pileup_ChIPPE.bedpe.bdg &> ${OUTPUTDIR_PREFIX}_run_pileup/run_pileup_ChIPPE.bedpe.log
  88 macs3 pileup -f BEDPE -i $CTRLBEDPE --outdir ${OUTPUTDIR_PREFIX}_run_pileup -o run_pileup_CTRLPE.bedpe.bdg &> ${OUTPUTDIR_PREFIX}_run_pileup/run_pileup_CTRLPE.bedpe.log
  89
  90 # filterdup
  91 echo "3. filterdup"
  92
  93 mkdir ${OUTPUTDIR_PREFIX}_run_filterdup
  94
  95 macs3 filterdup -g 52000000 -i $CHIP --outdir ${OUTPUTDIR_PREFIX}_run_filterdup -o run_filterdup_result.bed --dry-run &> ${OUTPUTDIR_PREFIX}_run_filterdup/run_filterdup_d.log
  96 macs3 filterdup -g 52000000 -i $CHIP --outdir ${OUTPUTDIR_PREFIX}_run_filterdup -o run_filterdup_result.bed &> ${OUTPUTDIR_PREFIX}_run_filterdup/run_filterdup.log
  97
  98 macs3 filterdup -g 52000000 -i $CHIPPE -f BAMPE --outdir ${OUTPUTDIR_PREFIX}_run_filterdup -o run_filterdup_result_pe.bedpe --dry-run &> ${OUTPUTDIR_PREFIX}_run_filterdup/run_filterdup_pe_d.log
  99 macs3 filterdup -g 52000000 -i $CHIPPE -f BAMPE --outdir ${OUTPUTDIR_PREFIX}_run_filterdup -o run_filterdup_result_pe.bedpe &> ${OUTPUTDIR_PREFIX}_run_filterdup/run_filterdup_pe.log
 100
 101 # predictd
 102 echo "4. predictd"
 103
 104 mkdir ${OUTPUTDIR_PREFIX}_run_predictd
 105
 106 macs3 predictd -i $CHIPPE -f BAMPE --outdir ${OUTPUTDIR_PREFIX}_run_predictd &> ${OUTPUTDIR_PREFIX}_run_predictd/run_predictd_bampe.log
 107 grep "Average insertion length" ${OUTPUTDIR_PREFIX}_run_predictd/run_predictd_bampe.log | perl -pe 's/^.*\s(\d+)\sbps/$1/' > ${OUTPUTDIR_PREFIX}_run_predictd/run_predictd_bampe.txt
 108
 109 macs3 predictd -i $CHIPBEDPE -f BEDPE --outdir ${OUTPUTDIR_PREFIX}_run_predictd &> ${OUTPUTDIR_PREFIX}_run_predictd/run_predictd_bedpe.log
 110 grep "Average insertion length" ${OUTPUTDIR_PREFIX}_run_predictd/run_predictd_bedpe.log | perl -pe 's/^.*\s(\d+)\sbps/$1/' > ${OUTPUTDIR_PREFIX}_run_predictd/run_predictd_bedpe.txt
 111
 112 macs3 predictd -g 52000000 -i $CHIP --d-min 10 --outdir ${OUTPUTDIR_PREFIX}_run_predictd --rfile run_predictd.R &> ${OUTPUTDIR_PREFIX}_run_predictd/run_predictd.log
 113 grep "predicted fragment length" ${OUTPUTDIR_PREFIX}_run_predictd/run_predictd.log | perl -pe 's/^.*\s(\d+)\sbps/$1/' > ${OUTPUTDIR_PREFIX}_run_predictd/run_predictd.txt
 114 grep "alternative fragment length" ${OUTPUTDIR_PREFIX}_run_predictd/run_predictd.log | perl -pe 's/^.*\s(\S+)\sbps/$1/' >> ${OUTPUTDIR_PREFIX}_run_predictd/run_predictd.txt
 115
 116 # randsample
 117 echo "5. randsample"
 118
 119 mkdir ${OUTPUTDIR_PREFIX}_run_randsample
 120
 121 macs3 randsample -i $CHIP -n 10000 --seed 31415926 --outdir ${OUTPUTDIR_PREFIX}_run_randsample -o run_randsample.bed &> ${OUTPUTDIR_PREFIX}_run_randsample/run_randsample.log
 122
 123 macs3 randsample -f BAMPE -i $CHIPPE -n 10000 --seed 31415926 --outdir ${OUTPUTDIR_PREFIX}_run_randsample -o run_randsample_bampe.bedpe &> ${OUTPUTDIR_PREFIX}_run_randsample/run_randsample_bampe.log
 124
 125 macs3 randsample -f BEDPE -i $CHIPBEDPE -n 10000 --seed 31415926 --outdir ${OUTPUTDIR_PREFIX}_run_randsample -o run_randsample_bedpe.bedpe &> ${OUTPUTDIR_PREFIX}_run_randsample/run_randsample_bedpe.log
 126
 127 # refinepeak
 128 echo "6. refinepeak"
 129
 130 mkdir ${OUTPUTDIR_PREFIX}_run_refinepeak
 131
 132 macs3 refinepeak -b ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_peaks.narrowPeak -i $CHIP --outdir ${OUTPUTDIR_PREFIX}_run_refinepeak --o-prefix run_refinepeak_w_prefix &> ${OUTPUTDIR_PREFIX}_run_refinepeak/run_refinepeak_w_prefix.log
 133
 134 macs3 refinepeak -b ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_peaks.narrowPeak -i $CHIP --outdir ${OUTPUTDIR_PREFIX}_run_refinepeak -o run_refinepeak_w_ofile.bed &> ${OUTPUTDIR_PREFIX}_run_refinepeak/run_refinepeak_w_ofile.log
 135
 136 # bdgcmp
 137 echo "7. bdgcmp"
 138
 139 mkdir ${OUTPUTDIR_PREFIX}_run_bdgcmp
 140
 141 macs3 bdgcmp -t ${OUTPUTDIR_PREFIX}_run_pileup/run_pileup_ChIP.bed.bdg -c ${OUTPUTDIR_PREFIX}_run_pileup/run_pileup_CTRL.bed.bdg  -m ppois FE -p 1 --outdir ${OUTPUTDIR_PREFIX}_run_bdgcmp --o-prefix run_bdgcmp &> ${OUTPUTDIR_PREFIX}_run_bdgcmp/run_bdgcmp.log
 142
 143 # bdgpeakcall
 144 echo "8. bdgpeakcall"
 145
 146 mkdir ${OUTPUTDIR_PREFIX}_run_bdgpeakcall
 147
 148 macs3 bdgpeakcall -i ${OUTPUTDIR_PREFIX}_run_bdgcmp/run_bdgcmp_FE.bdg -c 2 --outdir ${OUTPUTDIR_PREFIX}_run_bdgpeakcall --o-prefix run_bdgpeakcall_w_prefix &> ${OUTPUTDIR_PREFIX}_run_bdgpeakcall/run_bdgpeakcall_w_prefix.log
 149
 150 macs3 bdgpeakcall -i ${OUTPUTDIR_PREFIX}_run_bdgcmp/run_bdgcmp_FE.bdg -c 2 --outdir ${OUTPUTDIR_PREFIX}_run_bdgpeakcall -o run_bdgpeakcall_cutoff.txt --cutoff-analysis &> ${OUTPUTDIR_PREFIX}_run_bdgpeakcall/run_bdgpeakcall_cutoff.log
 151
 152 # bdgbroadcall
 153 echo "9. bdgbroadcall"
 154
 155 mkdir ${OUTPUTDIR_PREFIX}_run_bdgbroadcall
 156
 157 macs3 bdgbroadcall -i ${OUTPUTDIR_PREFIX}_run_bdgcmp/run_bdgcmp_FE.bdg -c 2 -C 1.5 --outdir ${OUTPUTDIR_PREFIX}_run_bdgbroadcall --o-prefix run_bdgbroadcall_w_prefix &> ${OUTPUTDIR_PREFIX}_run_bdgbroadcall/run_bdgbroadcall_w_prefix.log
 158
 159 # bdgdiff
 160 echo "10. bdgdiff"
 161
 162 mkdir ${OUTPUTDIR_PREFIX}_run_callpeak_narrow_revert
 163 mkdir ${OUTPUTDIR_PREFIX}_run_bdgdiff
 164
 165 macs3 callpeak -g 10000000 --nomodel --extsize 250 -c $CHIP -t $CTRL -n run_callpeak_narrow_revert -B --outdir ${OUTPUTDIR_PREFIX}_run_callpeak_narrow_revert &> ${OUTPUTDIR_PREFIX}_run_callpeak_narrow_revert/run_callpeak_narrow_revert.log
 166
 167 macs3 bdgdiff --t1 ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_treat_pileup.bdg --c1 ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_control_lambda.bdg --t2 ${OUTPUTDIR_PREFIX}_run_callpeak_narrow_revert/run_callpeak_narrow_revert_treat_pileup.bdg --c2 ${OUTPUTDIR_PREFIX}_run_callpeak_narrow_revert/run_callpeak_narrow_revert_control_lambda.bdg --o-prefix run_bdgdiff_prefix --outdir ${OUTPUTDIR_PREFIX}_run_bdgdiff &> ${OUTPUTDIR_PREFIX}_run_bdgdiff/run_bdgdiff_w_prefix.log
 168
 169 macs3 bdgdiff --t1 ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_treat_pileup.bdg --c1 ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_control_lambda.bdg --t2 ${OUTPUTDIR_PREFIX}_run_callpeak_narrow_revert/run_callpeak_narrow_revert_treat_pileup.bdg --c2 ${OUTPUTDIR_PREFIX}_run_callpeak_narrow_revert/run_callpeak_narrow_revert_control_lambda.bdg -o cond1.bed cond2.bed common.bed --outdir ${OUTPUTDIR_PREFIX}_run_bdgdiff &> ${OUTPUTDIR_PREFIX}_run_bdgdiff/run_bdgdiff_w_o_file.log
 170
 171 # cmbreps
 172 echo "11. cmbreps"
 173
 174 mkdir ${OUTPUTDIR_PREFIX}_run_cmbreps
 175
 176 macs3 cmbreps -i ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_treat_pileup.bdg ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_control_lambda.bdg ${OUTPUTDIR_PREFIX}_run_bdgcmp/run_bdgcmp_ppois.bdg -m max -o run_cmbreps_max.bdg --outdir ${OUTPUTDIR_PREFIX}_run_cmbreps &> ${OUTPUTDIR_PREFIX}_run_cmbreps/run_cmbreps_max.log
 177 macs3 cmbreps -i ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_treat_pileup.bdg ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_control_lambda.bdg ${OUTPUTDIR_PREFIX}_run_bdgcmp/run_bdgcmp_ppois.bdg -m mean -o run_cmbreps_mean.bdg --outdir ${OUTPUTDIR_PREFIX}_run_cmbreps &> ${OUTPUTDIR_PREFIX}_run_cmbreps/run_cmbreps_mean.log
 178 macs3 cmbreps -i ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_treat_pileup.bdg ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_control_lambda.bdg ${OUTPUTDIR_PREFIX}_run_bdgcmp/run_bdgcmp_ppois.bdg -m fisher -o run_cmbreps_fisher.bdg --outdir ${OUTPUTDIR_PREFIX}_run_cmbreps &> ${OUTPUTDIR_PREFIX}_run_cmbreps/run_cmbreps_fisher.log
 179
 180 # bdgopt
 181 echo "12. bdgopt"
 182
 183 mkdir ${OUTPUTDIR_PREFIX}_run_bdgopt
 184
 185 macs3 bdgopt -m min -p 10 -i ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_treat_pileup.bdg -o run_bdgopt_min.bdg --outdir ${OUTPUTDIR_PREFIX}_run_bdgopt &> ${OUTPUTDIR_PREFIX}_run_bdgopt/run_bdgopt_min.log
 186
 187 macs3 bdgopt -m max -p 2 -i ${OUTPUTDIR_PREFIX}_run_callpeak_narrow/run_callpeak_narrow0_treat_pileup.bdg -o run_bdgopt_max.bdg --outdir ${OUTPUTDIR_PREFIX}_run_bdgopt &> ${OUTPUTDIR_PREFIX}_run_bdgopt/run_bdgopt_max.log
 188
 189 # callvar
 190 echo "13. callvar"
 191
 192 mkdir ${OUTPUTDIR_PREFIX}_run_callvar
 193
 194 macs3 callvar -b ${CALLVARPEAK} -t ${CHIPPE} -c ${CTRLPE} -o ${OUTPUTDIR_PREFIX}_run_callvar/PEsample.vcf &> ${OUTPUTDIR_PREFIX}_run_callvar/run_callvar_PE.log
 195
 196 # test large amount of contigs
 197 echo "14. 50k contigs with buffersize"
 198
 199 mkdir ${OUTPUTDIR_PREFIX}_run_50kcontigs
 200
 201 echo "14.1 callpeak"
 202 macs3 callpeak -g 10000000 -t $CHIPCONTIGS50K -n run_callpeak_50kcontigs --outdir ${OUTPUTDIR_PREFIX}_run_50kcontigs --buffer-size 1000 --nomodel --extsize 200 &> ${OUTPUTDIR_PREFIX}_run_50kcontigs/run_callpeak_50kcontigs.log
 203
 204 echo "14.2 filterdup"
 205 macs3 filterdup -g 10000000 -i $CHIPCONTIGS50K --outdir ${OUTPUTDIR_PREFIX}_run_50kcontigs -o run_filterdup_result.bed --buffer-size 1000 &> ${OUTPUTDIR_PREFIX}_run_50kcontigs/run_filterdup_50kcontigs.log
 206
 207 echo "14.3 pileup"
 208 macs3 pileup -f BED -i $CHIPCONTIGS50K --extsize 200 --outdir ${OUTPUTDIR_PREFIX}_run_50kcontigs -o run_pileup_ChIP.bed.bdg --buffer-size 1000 &> ${OUTPUTDIR_PREFIX}_run_50kcontigs/run_pileup_ChIP.bed.log
 209
 210 echo "14.4 randsample"
 211 macs3 randsample -i $CHIPCONTIGS50K -n 100000 --seed 31415926 --outdir ${OUTPUTDIR_PREFIX}_run_50kcontigs -o run_randsample.bed --buffer-size 1000 &> ${OUTPUTDIR_PREFIX}_run_50kcontigs/run_randsample.log
 212
 213 echo "15. hmmratac"
 214 mkdir ${OUTPUTDIR_PREFIX}_run_hmmratac
 215
 216 macs3 hmmratac -i $ATACSEQBAM -n hmmratac_yeast500k --save-training-data --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac.log
 217
 218 macs3 hmmratac -i $ATACSEQBED -n hmmratac_yeast500k_bedpe -f BEDPE --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_bedpe.log
 219
 220 TRAINING_REGIONS_BED=${OUTPUTDIR_PREFIX}_run_hmmratac/hmmratac_yeast500k_training_regions.bed
 221 HMM_MODEL=${OUTPUTDIR_PREFIX}_run_hmmratac/hmmratac_yeast500k_model.json
 222
 223 echo "15.1 hmmratac load training regions from bedfile"
 224 macs3 hmmratac -i $ATACSEQBAM -n hmmratac_yeast500k_load_training_regions -t ${TRAINING_REGIONS_BED} --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_load_training_regions.log
 225
 226 echo "15.2 hmmratac load hmm model file"
 227 macs3 hmmratac -i $ATACSEQBAM -n hmmratac_yeast500k_load_hmm_model --model ${HMM_MODEL} --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_load_hmm_model.log
 228
 229 echo "16. search for errors or warnings in log files"
 230 flag=0
 231 for i in `ls ${OUTPUTDIR_PREFIX}_run_*/*.log`;do
 232     echo " checking $i..."
 233     egrep -i "error|warning" $i;
 234     if [[ $? == 0 ]];then
 235         echo " ... error/warning found!"
 236         echo "the log file:"
 237         cat $i;
 238         let "flag=1";
 239     else
 240         echo " ... clear!"
 241     fi
 242 done;
 243
 244 if [[ $flag == 1 ]]; then
 245     echo " Error or Warning can be found! Quit the test!"
 246     exit 1;
 247 fi
 248
 249 echo "17. compare with standard outputs"
 250
 251 flag=0
 252 # we do not check the consistency of random subsampling results
 253 subfolders=(50kcontigs bdgbroadcall bdgcmp bdgdiff bdgpeakcall callpeak_broad callpeak_narrow callpeak_narrow_revert callpeak_pe_broad callpeak_pe_narrow cmbreps filterdup pileup refinepeak bdgopt predictd callvar hmmratac)
 254
 255 m=0
 256 n=0
 257 for i in ${subfolders[@]};do
 258     let "m+=1"
 259     let "n=0"
 260     for j in `ls standard_results_${i}/*`; do
 261         k=`basename $j`
 262         let "n+=1"
 263         echo "16.${m}.${n} checking $i $k ..."
 264         fq=${OUTPUTDIR_PREFIX}_run_${i}/${k}
 265         fs=standard_results_${i}/${k}
 266         echo "  checking the files: $fq vs $fs"
 267         if [ "${REPLACE}" = "Y" ]; then
 268             cp $fq $fs
 269             echo " ... replaced!"
 270         else
 271             sort $fq| grep -v ^# > tmp_fq.txt
 272             sort $fs| grep -v ^# > tmp_fs.txt
 273             d=`diff tmp_fq.txt tmp_fs.txt`
 274             if [ -z "$d" ]; then
 275                 # the two files are the same
 276                 echo " ... success! Results are the same!"
 277             else
 278                 # if the two files are not the same, and if the file
 279                 # is a peak file, here we will check if the files are
 280                 # the similar using jaccard index
 281                 if [[ "$fq" == *.bed || "$fq" == *.gappedPeak || "$fq" == *.narrowPeak || "$fq" == *.bed12 ]]; then
 282                     # for peaks, we will use our jaccard index script to check
 283                     ji=$(python ./jaccard.py "$fq" "$fs")
 284                     # Compare the output value
 285                     if (( $(echo "$ji > 0.99" | bc -l) )); then
 286                         echo " ... success! Results are different but Jaccard Index is $ji (>0.99)"
 287                     else
 288                         echo " ... failed! Results are different and Jaccard Index is $ji (<0.99)"
 289                         let "flag=1"
 290                     fi
 291                 else
 292                     echo " ... failed! Results are different (top10 lines of difference):"
 293                     diff tmp_fq.txt tmp_fs.txt | head -10
 294                     let "flag=1";
 295                 fi
 296             fi
 297         fi
 298     done
 299 done
 300 rm -f tmp_fq.txt tmp_fs.txt
 301
 302 # exit with 1 if test fails
 303 if [ $flag -eq 1 ];then
 304     exit $flag;
 305 fi
 306
 307 echo "18. brief speed test"
 308
 309 mkdir ${OUTPUTDIR_PREFIX}_speedtest
 310 macs3 callpeak -t CTCF_12878_5M.bed.gz -c Input_12878_5M.bed.gz -n speedtest -B --outdir ${OUTPUTDIR_PREFIX}_speedtest &> ${OUTPUTDIR_PREFIX}_speedtest/speedtest.log &
 311 # monitor the macs3 callpeak run and output CPU and mem usage
 312 ./prockreport 0.5 $!
 313 wait;
 314
 315 # END of test