|
1
|
1 ########################################################################
|
|
|
2 # JAMMv1.0.7rev1 is a peak finder for joint analysis of NGS replicates.
|
|
|
3 # Copyright (C) 2014-2015 Mahmoud Ibrahim
|
|
|
4 #
|
|
|
5 # This program is free software: you can redistribute it and/or modify
|
|
|
6 # it under the terms of the GNU General Public License as published by
|
|
|
7 # the Free Software Foundation, either version 3 of the License, or
|
|
|
8 # (at your option) any later version.
|
|
|
9 #
|
|
|
10 # This program is distributed in the hope that it will be useful,
|
|
|
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
13 # GNU General Public License for more details.
|
|
|
14 #
|
|
|
15 # You should have received a copy of the GNU General Public License
|
|
|
16 # along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
|
17 #
|
|
|
18 # Contact: mahmoud.ibrahim@mdc-berlin.de
|
|
|
19 ########################################################################
|
|
0
|
20
|
|
|
21
|
|
|
22 ##Finding out the path
|
|
|
23 sPath="`dirname \"$0\"`"
|
|
|
24 sPath="`( cd \"$sPath\" && pwd )`"
|
|
|
25
|
|
|
26
|
|
|
27
|
|
|
28 usage()
|
|
|
29 {
|
|
|
30 cat << EOF
|
|
1
|
31 Welcome to JAMM v1.0.7rev1 (GNU GPLv3). Copyright (C) 2014-2015 Mahmoud Ibrahim.
|
|
0
|
32
|
|
|
33 This program comes with ABSOLUTELY NO WARRANTY; for details visit http://www.gnu.org/licenses/gpl.html. This is free software, and you are welcome to redistribute it under certain conditions; visit http://www.gnu.org/licenses/gpl.html for details.
|
|
|
34
|
|
|
35 OPTIONS:
|
|
|
36 -s directory containing Sample files (required)
|
|
|
37 -g Genome size file (required)
|
|
|
38 -o Output directory (required)
|
|
|
39 -c directory containing input or Control files
|
|
1
|
40 -f Fragment length(s) (default: estimated)
|
|
|
41 -r Resolution, peak or region or window (default: peak)
|
|
0
|
42 -m Mode, normal or narrow (default: normal)
|
|
1
|
43 -i clustering Initialization window selection, deterministic or stochastic (default: deterministic)
|
|
|
44 -b Bin Size (default: estimated)
|
|
0
|
45 -w minimum Window size (default: 2 --- Note: this means minimum_window_size = bin_size x the_value_of_-w)
|
|
1
|
46 -e window Enrichment cutoff, auto or any numeric value (default: 1 --- Set this to "auto" to estimate the window enrichment cutoff)
|
|
|
47 -d keep PCR Dupicates in single-end mode, y or n (default: n --- if -t is "paired", this option has no effect)
|
|
|
48 -t Type, single or paired (default: single, requires BED files. paired requires BEDPE files)
|
|
0
|
49 -p Number of processors used by R scripts (default: 1)
|
|
1
|
50
|
|
0
|
51 EOF
|
|
|
52 }
|
|
|
53
|
|
|
54
|
|
|
55 # =========================
|
|
|
56 # Process Input parameters
|
|
|
57 # =========================
|
|
|
58
|
|
|
59
|
|
|
60 #Defaults -- Change those if you want
|
|
|
61 mode="normal"
|
|
|
62 resol="peak"
|
|
|
63 cores="1"
|
|
|
64 window="2"
|
|
|
65 type="single"
|
|
|
66 windowe="1"
|
|
1
|
67 initModel="deterministic"
|
|
|
68 uniq="n"
|
|
0
|
69
|
|
|
70 #Defaults -- Do not change
|
|
|
71 sdir=""
|
|
|
72 gsize=""
|
|
|
73 out=""
|
|
|
74 binsize="ns"
|
|
|
75 fraglen="ns"
|
|
|
76 ran=$RANDOM
|
|
|
77 wdir=$(mktemp -d)
|
|
1
|
78 export LANG=C #locale defaults
|
|
|
79 export LC_ALL=C #locale defaults
|
|
0
|
80
|
|
1
|
81 while getopts "s:g:o:c:m:r:f:p:w:b:t:e:i:d:" OPTION
|
|
0
|
82 do
|
|
|
83 case $OPTION in
|
|
|
84 s) sdir=$OPTARG
|
|
|
85 ;;
|
|
|
86 g) gsize=$OPTARG
|
|
|
87 ;;
|
|
|
88 o) out=$OPTARG
|
|
|
89 ;;
|
|
|
90 c) bdir=$OPTARG
|
|
|
91 ;;
|
|
|
92 m) mode=$OPTARG
|
|
|
93 ;;
|
|
|
94 r) resol=$OPTARG
|
|
|
95 ;;
|
|
|
96 f) fraglen=$OPTARG
|
|
|
97 ;;
|
|
|
98 p) cores=$OPTARG
|
|
|
99 ;;
|
|
|
100 w) window=$OPTARG
|
|
|
101 ;;
|
|
|
102 b) binsize=$OPTARG
|
|
|
103 ;;
|
|
|
104 t) type=$OPTARG
|
|
|
105 ;;
|
|
|
106 e) windowe=$OPTARG
|
|
|
107 ;;
|
|
1
|
108 i) initModel=$OPTARG
|
|
|
109 ;;
|
|
|
110 d) uniq=$OPTARG
|
|
|
111 ;;
|
|
0
|
112 ?)
|
|
|
113 usage
|
|
|
114 exit
|
|
|
115 ;;
|
|
|
116 esac
|
|
|
117 done
|
|
|
118 if [ "$mode" == "normal" ]; then
|
|
|
119 clustno="2"
|
|
|
120 fi
|
|
|
121 if [ "$mode" == "narrow" ]; then
|
|
|
122 clustno="3"
|
|
|
123 fi
|
|
|
124
|
|
|
125 if [[ -z $sdir ]] || [[ -z $gsize ]] || [[ -z $out ]]
|
|
|
126 then
|
|
|
127 usage
|
|
|
128 exit 1
|
|
|
129 fi
|
|
1
|
130 if [[ -d "$out/peaks" ]]; then
|
|
|
131 printf "\n\nOutput directory $out/peaks already exists. I can't override existing results!\n\n"
|
|
|
132 exit 0
|
|
|
133 fi
|
|
|
134 if [ $fraglen == "ns" ]; then
|
|
|
135 if [[ -d "$out/xcorr" ]]; then
|
|
|
136 printf "\n\nOutput directory $out/xcorr already exists. I can't override existing results!\n\n"
|
|
|
137 exit 0
|
|
|
138 fi
|
|
|
139 fi
|
|
0
|
140 #=======================> DONE!
|
|
|
141
|
|
|
142
|
|
|
143
|
|
1
|
144
|
|
0
|
145 # =============================
|
|
|
146 # Step One: Initial Processing
|
|
|
147 # =============================
|
|
1
|
148 printf "\n\n============================================\nStarted JAMM Pipeline v1.0.7rev1...Hang on!\n============================================\n\n"
|
|
0
|
149
|
|
|
150 if [ ! -d "$wdir" ]; then
|
|
|
151 mkdir $wdir #make working directory
|
|
|
152 fi
|
|
|
153 if [ ! -d "$out" ]; then
|
|
|
154 mkdir $out #make output directory
|
|
|
155 fi
|
|
|
156 mkdir $wdir/bkgd.$ran/ #directory to store background files
|
|
|
157 mkdir $wdir/sizes.$ran/ #chromosomes and sizes
|
|
|
158 mkdir $wdir/samples.$ran/ #store sample files
|
|
|
159
|
|
|
160 dupnum=$(ls -1 $sdir | wc -l) #count how many sample files
|
|
|
161
|
|
|
162
|
|
|
163 #separate chromosome sizes
|
|
|
164 printf "Loading genome size file..."
|
|
|
165 ext="$wdir/sizes.$ran/"
|
|
|
166 awk -v ext="$ext" '{ print >> ext"/size." $1 ".bed" }' $gsize
|
|
|
167 printf "Done!\n"
|
|
|
168
|
|
|
169
|
|
|
170 printf "Processing sample files..."
|
|
|
171 #load each chromosome from each sample file
|
|
|
172 for i in $sdir/*.bed; do
|
|
|
173 samplefile=$(basename $i)
|
|
|
174 for f in $wdir/sizes.$ran/*; do
|
|
|
175 sizefile=$(basename $f)
|
|
|
176 chr=$(echo $sizefile | awk -F"." '{print $2}' | awk -F"." '{print $1}');
|
|
|
177 awk -v chr="$chr" -v ext="$wdir/samples.$ran/" -v samplefile="$samplefile" -F"\t" '$1 == chr { print $2"\t"$6 >> ext"sample."chr"."samplefile }' "$i"
|
|
|
178 done
|
|
|
179 done
|
|
|
180 printf "Done!\n"
|
|
|
181
|
|
|
182
|
|
|
183 if [ ! -z $bdir ]; then
|
|
|
184 #concatenate all background files into one file
|
|
|
185 printf "Processing control files..."
|
|
|
186 cat $bdir/*.bed > $wdir/bkgd.$ran/ctrl.bed
|
|
|
187
|
|
|
188 for f in $wdir/sizes.$ran/*; do
|
|
|
189 sizefile=$(basename $f)
|
|
|
190 chr=$(echo $sizefile | awk -F"." '{print $2}' | awk -F"." '{print $1}');
|
|
|
191 awk -v chr="$chr" -v ext="$wdir/bkgd.$ran/" -F"\t" '$1 == chr { print $2"\t"$6 >> ext"bkgd."chr".ctrl.bed" }' "$wdir/bkgd.$ran/ctrl.bed"
|
|
|
192 done
|
|
|
193
|
|
|
194 printf "Done!\n"
|
|
|
195 fi
|
|
|
196
|
|
|
197 #determine average read lengths
|
|
1
|
198 printf "Getting average read lengths...\n"
|
|
|
199 readL=""
|
|
|
200 if [ ! -z $bdir ]; then
|
|
|
201 readC=$(awk '{a=$3-$2;print a;}' "$wdir/bkgd.$ran/ctrl.bed" | perl -lane '$a+=$_;END{print $a/$.}' | awk '{a=$1+0.5;print a;}' | cut -d"." -f1)
|
|
|
202 printf "Control: $readC\n"
|
|
|
203 fi
|
|
|
204 readL=""
|
|
|
205 for s in $sdir/*.bed; do #and for each sample file
|
|
|
206 file=$(basename $s)
|
|
|
207 samplefile=$(echo $file | awk -F"." '{print $1}');
|
|
|
208 read=$(awk '{a=$3-$2;print a;}' "$s" | perl -lane '$a+=$_;END{print $a/$.}' | awk '{a=$1+0.5;print a;}' | cut -d"." -f1)
|
|
|
209 printf "$samplefile: $read\n"
|
|
|
210 readL="$readL,$read"
|
|
|
211 done
|
|
|
212 readL=${readL#","}
|
|
0
|
213 #=======================> DONE!
|
|
|
214
|
|
|
215
|
|
|
216
|
|
|
217 # =============================
|
|
|
218 # Step Two: Fragment Length
|
|
|
219 # =============================
|
|
|
220 #single-end
|
|
|
221 if [ $type == "single" ]; then
|
|
|
222
|
|
|
223 if [ $fraglen == "ns" ]; then
|
|
|
224
|
|
|
225 ##Counting Where Reads Start and Calculating Cross Correlation
|
|
|
226 mkdir $wdir/stats.$ran/ #store count files
|
|
|
227 mkdir $out/xcorr #final xcorr results
|
|
|
228
|
|
|
229
|
|
|
230 printf "Calculating Fragment Length(s)...\n"
|
|
|
231 for f in $wdir/sizes.$ran/*; do #for each chromosome
|
|
|
232 samplelist=""
|
|
|
233 readlist=""
|
|
|
234 sizefile=$(basename $f)
|
|
|
235 chr=$(echo $sizefile | awk -F"." '{print $2}' | awk -F"." '{print $1}');
|
|
|
236
|
|
|
237 #list of sample bed files and read lengths
|
|
|
238 for s in $wdir/samples.$ran/*.bed; do #and for each sample file
|
|
|
239 samplefile=$(basename $s)
|
|
|
240 chr2=$(echo $samplefile | awk -F"." '{print $2}');
|
|
|
241 if [ $chr == $chr2 ] #belonging to this chromosome
|
|
|
242 then
|
|
|
243 samplelist="$samplelist,$wdir/samples.$ran/$samplefile"
|
|
|
244 fi
|
|
|
245 done
|
|
|
246 readlist="$readL"
|
|
|
247
|
|
|
248 #list of control bed files and read lengths
|
|
|
249 if [ ! -z $bdir ]; then
|
|
|
250 for s in $wdir/bkgd.$ran/*.bed; do #and for each sample file
|
|
|
251 samplefile=$(basename $s)
|
|
|
252 chr2=$(echo $samplefile | awk -F"." '{print $2}');
|
|
|
253 if [ $chr == $chr2 ] #belonging to this chromosome
|
|
|
254 then
|
|
|
255 samplelist="$samplelist,$wdir/bkgd.$ran/$samplefile"
|
|
|
256 readlist="$readL,$readC"
|
|
|
257 fi
|
|
|
258 done
|
|
|
259 fi
|
|
|
260
|
|
|
261 #remove leading comma
|
|
|
262 samplelist=${samplelist#","}
|
|
|
263
|
|
|
264 #call R script for xcorr calculation
|
|
|
265 Rscript "$sPath/xcorr.r" -ibed="$samplelist" -s="$wdir/sizes.$ran/size.$chr.bed" -rl="$readlist" -d="$wdir/stats.$ran" -p="$cores"
|
|
|
266 done
|
|
|
267
|
|
|
268 #report xcorr results (samples)
|
|
|
269 for f in $sdir/*.bed; do
|
|
|
270 file=$(basename $f)
|
|
|
271 samplefile=$(echo $file | awk -F"." '{print $1}');
|
|
|
272 mkdir "$out/xcorr/$samplefile" #final xcorr results
|
|
|
273 if [ -f "$wdir/stats.$ran/xc.$samplefile.tab" ]; then
|
|
|
274 cp $wdir/stats.$ran/xc.$samplefile.tab $out/xcorr/$samplefile/shifts.txt
|
|
|
275 fi
|
|
|
276 Rscript "$sPath/xcorrhelper.r" -infile="$out/xcorr/$samplefile/shifts.txt" -out="$out/xcorr/$samplefile"
|
|
|
277 done
|
|
|
278 #report xcorr results (control)
|
|
|
279 if [ ! -z $bdir ]; then
|
|
|
280 mkdir "$out/xcorr/ctrl" #final xcorr results
|
|
|
281 if [ -f "$wdir/stats.$ran/xc.ctrl.tab" ]; then
|
|
|
282 cp $wdir/stats.$ran/xc.ctrl.tab $out/xcorr/ctrl/shifts.txt
|
|
|
283 fi
|
|
|
284 Rscript "$sPath/xcorrhelper.r" -infile="$out/xcorr/ctrl/shifts.txt" -out="$out/xcorr/ctrl"
|
|
|
285 fi
|
|
|
286
|
|
|
287 fi
|
|
|
288 fi
|
|
|
289
|
|
|
290 #paired-end
|
|
|
291 if [ $type == "paired" ]; then
|
|
|
292 printf "Getting Average Fragment Length(s)...\n"
|
|
|
293 mkdir $out/xcorr #final xcorr results
|
|
|
294
|
|
|
295 for f in $sdir/*.bed; do
|
|
|
296 file=$(basename $f)
|
|
|
297 samplefile=$(echo $file | awk -F"." '{print $1}');
|
|
|
298 mkdir "$out/xcorr/$samplefile"
|
|
|
299 frag=$(awk '{a=$6-$2;print a;}' $f | perl -lane '$a+=$_;END{print $a/$.}' | awk '{a=$1+0.5;print a;}' | cut -d"." -f1)
|
|
|
300 echo "Average_from_paired $frag" > $out/xcorr/$samplefile/shifts.txt
|
|
|
301 Rscript "$sPath/xcorrhelper.r" -infile="$out/xcorr/$samplefile/shifts.txt" -out="$out/xcorr/$samplefile"
|
|
|
302 done
|
|
|
303 fi
|
|
|
304 #=======================> DONE!
|
|
|
305
|
|
|
306
|
|
|
307
|
|
|
308 # =================================
|
|
|
309 # Step Three: Calculating Bin Size
|
|
|
310 # =================================
|
|
|
311 if [ $binsize == "ns" ]; then
|
|
|
312 printf "Getting Bin Size: "
|
|
|
313
|
|
|
314 chr=$(sort -nr -k2 $gsize | head -n 1 | awk -F"\t" '{print $1}');
|
|
|
315 samplelist=""
|
|
|
316 frag=""
|
|
|
317 if [ $fraglen != "ns" ]; then
|
|
|
318 frag=$fraglen
|
|
|
319 k=1
|
|
|
320 fi
|
|
|
321
|
|
|
322 #list of sample bed files and read lengths
|
|
|
323 for s in $wdir/samples.$ran/*.bed; do
|
|
|
324 samplefile=$(basename $s)
|
|
|
325 chr2=$(echo $samplefile | awk -F"." '{print $2}');
|
|
|
326 if [ $chr == $chr2 ]
|
|
|
327 then
|
|
|
328 samplelist="$samplelist,$wdir/samples.$ran/$samplefile"
|
|
|
329 samplename=$(echo $samplefile | awk -F"." '{ print $3 }')
|
|
|
330 samplefilename=$(echo $samplefile | cut -d'.' -f 3-)
|
|
|
331 if [ $fraglen == "ns" ]; then
|
|
|
332 shift=$(awk -F":" '$1 == "Fragment Length" { print $2 }' "$out/xcorr/$samplename/xcorrsummary.txt")
|
|
|
333 frag="$frag,$shift"
|
|
|
334 fi
|
|
|
335 fi
|
|
|
336 done
|
|
|
337 #remove leading comma
|
|
|
338 samplelist=${samplelist#","}
|
|
|
339 frag=${frag#","}
|
|
|
340 Rscript "$sPath/bincalculator.r" -ibed="$samplelist" -s="$gsize" -rl="$readL" -d="$wdir" -p="$cores" -f="$frag" -type="$type"
|
|
|
341 fi
|
|
|
342 if [ $binsize != "ns" ]; then
|
|
|
343 printf "You set a Bin Size: $binsize \n"
|
|
|
344 fi
|
|
|
345 #=======================> DONE!
|
|
|
346
|
|
|
347
|
|
|
348
|
|
|
349 # ===========================
|
|
|
350 # Step Four: Calling Peaks
|
|
|
351 # ===========================
|
|
|
352 mkdir $wdir/peaks.$ran/ #store count files
|
|
|
353 mkdir $out/peaks #store peak files
|
|
|
354
|
|
|
355 printf "Calling Peaks...(mode: $mode, resolution: $resol)\n"
|
|
|
356
|
|
|
357
|
|
|
358 #single-end reads
|
|
|
359 if [ $type == "single" ]; then
|
|
|
360
|
|
|
361 if [ $binsize == "ns" ]; then
|
|
|
362 binsize=$(cat "$wdir/binsize.txt")
|
|
|
363 fi
|
|
|
364
|
|
|
365 counting=1;
|
|
|
366 for f in $wdir/sizes.$ran/*; do #for each chromosome
|
|
|
367 samplelist=""
|
|
|
368 frag=""
|
|
|
369 k=1
|
|
|
370 if [ $fraglen != "ns" ]; then
|
|
|
371 frag=$fraglen
|
|
|
372 fi
|
|
|
373
|
|
|
374
|
|
|
375 sizefile=$(basename $f)
|
|
|
376 chr=$(echo $sizefile | awk -F"." '{print $2}' | awk -F"." '{print $1}');
|
|
|
377
|
|
|
378 printf "Chromosome $chr: "
|
|
|
379
|
|
|
380 #list of sample bed files and fragment lengths
|
|
|
381 for s in $wdir/samples.$ran/*.bed; do #and for each sample file
|
|
|
382 samplefile=$(basename $s)
|
|
|
383 chr2=$(echo $samplefile | awk -F"." '{print $2}');
|
|
|
384 if [ $chr == $chr2 ] #belonging to this chromosome
|
|
|
385 then
|
|
|
386 samplelist="$samplelist,$wdir/samples.$ran/ext.$samplefile"
|
|
|
387 samplename=$(echo $samplefile | awk -F"." '{ print $3 }')
|
|
|
388 samplefilename=$(echo $samplefile | cut -d'.' -f 3-)
|
|
|
389 if [ $fraglen == "ns" ]; then
|
|
|
390 shift=$(awk -F":" '$1 == "Fragment Length" { print $2 }' "$out/xcorr/$samplename/xcorrsummary.txt")
|
|
|
391 frag="$frag,$shift"
|
|
|
392 read=$(echo $readL | cut -f "$k" -d ",")
|
|
|
393 k=$(($k+1))
|
|
|
394 fi
|
|
|
395 if [ $fraglen != "ns" ]; then
|
|
|
396 shift=$(echo $frag | cut -f "$k" -d ",")
|
|
|
397 read=$(echo $readL | cut -f "$k" -d ",")
|
|
|
398 k=$(($k+1))
|
|
|
399 fi
|
|
1
|
400 if [ $uniq == "y" ]; then
|
|
|
401 perl "$sPath/readshifter.pl" "$wdir/samples.$ran/$samplefile" $shift $read > "$wdir/samples.$ran/ext.$samplefile"
|
|
|
402 fi
|
|
|
403 if [ $uniq == "n" ]; then
|
|
|
404 perl "$sPath/readshifter.pl" "$wdir/samples.$ran/$samplefile" $shift $read | uniq > "$wdir/samples.$ran/ext.$samplefile"
|
|
|
405 fi
|
|
0
|
406 fi
|
|
|
407 done
|
|
|
408
|
|
|
409 #control file
|
|
|
410 bkgdfile="None"
|
|
|
411 if [ ! -z $bdir ]; then
|
|
|
412 if [ $fraglen == "ns" ]; then
|
|
|
413 bshift=$(awk -F":" '$1 == "Fragment Length" { print $2 }' "$out/xcorr/ctrl/xcorrsummary.txt")
|
|
|
414 frag="$frag,$bshift"
|
|
|
415 fi
|
|
|
416 if [ $fraglen != "ns" ]; then
|
|
|
417 l=$(($dupnum+1))
|
|
|
418 bshift=$(echo $frag | cut -f "$l" -d ",")
|
|
|
419 fi
|
|
1
|
420
|
|
|
421 if [ $uniq == "y" ]; then
|
|
|
422 perl "$sPath/readshifter.pl" "$wdir/bkgd.$ran/bkgd.$chr.ctrl.bed" $bshift $readC > "$wdir/bkgd.$ran/ext.bkgd.$chr.ctrl.bed"
|
|
|
423 fi
|
|
|
424 if [ $uniq == "n" ]; then
|
|
|
425 perl "$sPath/readshifter.pl" "$wdir/bkgd.$ran/bkgd.$chr.ctrl.bed" $bshift $readC | uniq > "$wdir/bkgd.$ran/ext.bkgd.$chr.ctrl.bed"
|
|
|
426 fi
|
|
0
|
427 bkgdfile="$wdir/bkgd.$ran/ext.bkgd.$chr.ctrl.bed"
|
|
|
428 fi
|
|
|
429
|
|
|
430 #remove leading comma
|
|
|
431 samplelist=${samplelist#","}
|
|
|
432 frag=${frag#","}
|
|
|
433
|
|
|
434 #call the peak calling R script
|
|
1
|
435 Rscript "$sPath/peakfinder.r" -sfile="$f" -bednames="$samplelist" -frag="$frag" -bkgd=$bkgdfile -out="$wdir/peaks.$ran/" -clustnummer="$clustno" -resolution="$resol" -window="$window" -p="$cores" -bin="$binsize" -type="$type" -chrcount="$counting" -initModel="$initModel" -windowe="$windowe"
|
|
0
|
436 counting=$(($counting+1));
|
|
|
437 cp "$wdir/peaks.$ran/$chr.peaks.bed" "$out/peaks/$chr.peaks.bed"
|
|
|
438 rm "$wdir/peaks.$ran/$chr.peaks.bed"
|
|
|
439 done
|
|
|
440 counting=1;
|
|
|
441 fi
|
|
|
442
|
|
|
443
|
|
|
444 #paired-end reads
|
|
|
445 if [ $type == "paired" ]; then
|
|
|
446
|
|
|
447 if [ $binsize == "ns" ]; then
|
|
|
448 binsize=$(cat "$wdir/binsize.txt")
|
|
|
449 fi
|
|
|
450
|
|
|
451
|
|
|
452 counting=1;
|
|
|
453 for f in $wdir/sizes.$ran/*; do #for each chromosome
|
|
|
454 samplelist=""
|
|
|
455
|
|
|
456 sizefile=$(basename $f)
|
|
|
457 chr=$(echo $sizefile | awk -F"." '{print $2}' | awk -F"." '{print $1}');
|
|
|
458
|
|
|
459 printf "Chromosome $chr: "
|
|
|
460
|
|
|
461 #list of sample bed files and fragment lengths
|
|
|
462 for s in $wdir/samples.$ran/*.bed; do #and for each sample file
|
|
|
463 samplefile=$(basename $s)
|
|
|
464 chr2=$(echo $samplefile | awk -F"." '{print $2}');
|
|
|
465 if [ $chr == $chr2 ] #belonging to this chromosome
|
|
|
466 then
|
|
|
467 samplelist="$samplelist,$wdir/samples.$ran/$samplefile"
|
|
|
468 samplename=$(echo $samplefile | awk -F"." '{ print $3 }')
|
|
|
469 samplefilename=$(echo $samplefile | cut -d'.' -f 3-)
|
|
|
470 x="$sdir/$samplefilename"
|
|
|
471 fi
|
|
|
472 done
|
|
|
473
|
|
|
474 #control file
|
|
|
475 bkgdfile="None"
|
|
|
476 if [ ! -z $bdir ]; then
|
|
|
477 bkgdfile="$wdir/bkgd.$ran/bkgd.$chr.ctrl.bed"
|
|
|
478 fi
|
|
|
479
|
|
|
480 #remove leading comma
|
|
|
481 samplelist=${samplelist#","}
|
|
|
482 frag=${frag#","}
|
|
|
483
|
|
|
484 #call the peak calling R script
|
|
1
|
485 Rscript "$sPath/peakfinder.r" -sfile=$f -bednames=$samplelist -frag="NA" -bkgd=$bkgdfile -out="$wdir/peaks.$ran/" -clustnummer="$clustno" -resolution="$resol" -window="$window" -p="$cores" -bin="$binsize" -type="$type" -chrcount="$counting" -initModel="$initModel" -windowe="$windowe"
|
|
0
|
486 counting=$(($counting+1));
|
|
|
487 cp "$wdir/peaks.$ran/$chr.peaks.bed" "$out/peaks/$chr.peaks.bed"
|
|
|
488 rm "$wdir/peaks.$ran/$chr.peaks.bed"
|
|
|
489 done
|
|
|
490 counting=1;
|
|
|
491 fi
|
|
|
492
|
|
|
493
|
|
|
494 cp $wdir/peaks.$ran/min.peaksize $out/peaks/min.peaksize
|
|
|
495
|
|
|
496
|
|
|
497 #concatenate, sort and filter
|
|
1
|
498 cat $out/peaks/*.bed > $out/peaks/all.narrowPeak
|
|
|
499 if [[ -s $out/peaks/all.narrowPeak ]]; then
|
|
|
500 Rscript "$sPath/peakhelper.r" -filelist="$out/peaks/all.narrowPeak"
|
|
|
501 perl "$sPath/peakfilter.pl" $out/peaks/all.narrowPeak | sort -nr -k7 > $out/peaks/filtered.peaks.narrowPeak
|
|
|
502 cut -f1-10 $out/peaks/all.narrowPeak | awk -F"\t" -v j=0 '$7 > j' | sort -nr -k7 > $out/peaks/all.peaks.narrowPeak
|
|
|
503 fi
|
|
0
|
504 rm $out/peaks/all.narrowPeak
|
|
|
505 rm $out/peaks/*.bed
|
|
|
506 rm $out/peaks/min.peaksize
|
|
|
507 #=======================> DONE!
|
|
|
508
|
|
|
509
|
|
|
510
|
|
|
511 rm -rf $wdir
|
|
|
512
|
|
|
513
|
|
|
514
|
|
|
515 printf "\n\n========================================\nWe're done...Congratulations!\n========================================\n\n"
|
