15 # Number of threads to run simultaneously.
18 my $CmdQueue = Thread::Queue->new;
24 if(!defined($root_dir) || !defined($data_dir)) {
25 print "Usage: ./ConfigureTasks.pm [root_dir] [data_dir] {option} \n";
31 my $BIOP = "$root_dir/bin/BioProspector.mac";
32 my $QUESTDIR = "$root_dir/bin/QuEST";
33 my $MACSDIR = "$root_dir/bin";
34 my $WINGPEAKSDIR = "$root_dir/bin";
35 my $WINGPEAKSGENOMEDIR = "$root_dir/reference_data";
36 my $PROFILEDIR = "$root_dir/bin";
37 my $GENOMEDIR = "/Volumes/Genomes";
38 my $QPCRDIR = "$root_dir/bin";
39 my $QPCRTESTDIR = "$root_dir/reference_data/qPCR_Tests";
40 my $QPCRBACKGROUND = "$root_dir/reference_data/GenericBackground";
42 my $SUMMARIZE_PROJECT="$root_dir/scripts/SummarizeProject2.pm";
45 ### Check directory for project and task folders. Make them if necessary
46 ### May not have permissions for this. These should be moved to an install script
47 `if [ ! -e $data_dir/Tasks ]; then mkdir $data_dir/Tasks; fi;`;
48 `if [ ! -e $data_dir/Projects ]; then mkdir $data_dir/Projects; fi;`;
52 ### Important: get the xml file form the server.
54 my $projects_file = getProjectsXML($parm);
56 my $xmldoc = XML::Simple->new();
58 my $xmlfile = $projects_file;
59 my $xml = $xmldoc->XMLin($xmlfile, ForceArray => ['Project','ComparePeakCalls','CompareLibraries','PeakCalling','ProfileReads','qPCR','MotifFinding'], KeepRoot=>1);
65 if(exists($xml->{Projects}->{Project})) {
66 for my $project (@{$xml->{Projects}->{Project}}) {
67 push @tasks, @{writeProject($project)};
69 my $projectname = $project->{Name};
70 my $projectId = $project->{ProjectId};
72 $projects .= "$projectId ";
76 print STDERR "Projects all complete\n";
78 my %saw; my $unique_tasks = ""; my $unique_task_dirs = "";
79 for(grep(!$saw{$_}++, @tasks)) {
80 $unique_tasks .= "$_ ";
81 $unique_task_dirs .= "../Tasks/$_ ";
84 print STDERR "writing tasks\n";
86 open(MAKE,">$data_dir/Projects/TaskMakefile");
87 print MAKE "all: $unique_task_dirs\n\n.PHONY: $unique_task_dirs\n\n";
88 for(split(/ /,$unique_tasks)) { print MAKE "../Tasks/$_:\n\t".'if ! $(MAKE) -C $@; then '.$root_dir.'/scripts/analys_track_main.py updsts '.$_.' "Error"; fi;'."\n\n"; }
91 print STDERR "Task makefile complete\n";
94 for(split(/ /,$projects)) { $index_list .= "$_/index.html "; }
96 open(MAKE,">$data_dir/Projects/ProjectMakefile");
97 print MAKE "all: $index_list | .start\n\n.PHONY: $index_list .start\n\n.start:\n\ttouch .start; echo \"Projects updated at `date`\";\n\n$index_list:\n\t".'cd `dirname $@` && '.$SUMMARIZE_PROJECT.' '.$root_dir.' '.$data_dir.'/LibraryInfo.xml > `basename $@`'."\n";
100 print STDERR "Project makefile compete. Starting to build tasks\n";
102 `cd $data_dir/Projects && make -j 4 -f TaskMakefile > make_tasks.log 2> make_tasks.err`;
106 my $projectname = $project->{Name}; $projectname =~ s/\s/_/g;
107 my $projectid = $project->{ProjectId};
109 my $projectdir = $projectid;
112 $tasks .= WriteQPCRTasks($project);
113 $tasks .= WriteProfileTasks($project);
114 $tasks .= WriteCompareLibTasks($project);
115 $tasks .= WriteMotifFindingTasks($project);
116 $tasks .= WritePeakCallingTasks($project);
118 print STDERR "Wrote tasks for project '$projectname' with tasks ($tasks)\n";
120 my @tasks = split(/ /,$tasks);
122 if(! (-e $projectdir && -d $projectdir)) {
123 `mkdir $data_dir/Projects/$projectdir`;
126 print STDERR "Created directory $projectdir\n";
128 $xmldoc->XMLout($project, OutputFile=>"$data_dir/Projects/$projectid/Project.xml", RootName=>"Project", XMLDecl=>1);
130 print STDERR "Wrote configureation for project $projectdir\n";
137 `$root_dir/scripts/analys_track_main.py updsts $taskid`;
138 print STDERR "Task $taskid registered.\n";
147 my $taskid = $task->{TaskId};
148 my $taskdir = "$data_dir/Tasks/".$taskid;
150 if(! (-e $taskdir && -d $taskdir)) {
154 $xmldoc->XMLout($task, OutputFile=>"$taskdir/Task.xml", RootName=>$root, XMLDecl=>1);
156 open(MAKEFILE, ">$taskdir/Makefile");
157 print MAKEFILE "all: .notify $outfile | .start\n\n.PHONY: .notify .start\n\n";
158 print MAKEFILE ".start:\n\t$root_dir/scripts/analys_track_main.py updsts $taskid \"Processing\"\n\ttouch .start\n\n";
160 print MAKEFILE ".notify: | .start $outfile .start\n\techo \"Complete\"\n\t$root_dir/scripts/analys_track_main.py updsts $taskid \"Complete\"\n\ttouch .notify\n\n";
161 print MAKEFILE "$cmd";
163 registerTask($taskid);
169 if(exists($project->{qPCR})) {
170 for my $qpcr (@{$project->{qPCR}}) {
171 my $task = $qpcr->{TaskId};
172 my $lib = $qpcr->{Library};
173 my $genome = $qpcr->{Genome};
174 my $name = $qpcr->{Name}; $name =~ tr/ \(\)\./____/;
176 my $background = $QPCRBACKGROUND;
177 my $testdir = $QPCRTESTDIR;
178 my $outfile = "$name.qPCR";
180 my $seqcheck = "if [ ! -e ~Data/Libraries/$lib.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for sequencing.\"; fi;";
181 my $cmd = "$outfile: ~Data/Libraries/$lib.txt\n\t$seqcheck\n\t$QPCRDIR/qPCR \$< $background $testdir > \$@\n";
182 writeTask($qpcr, "qPCR", $outfile, $cmd);
190 sub WriteProfileTasks {
193 if(exists($project->{ProfileReads})) {
194 for my $profile (@{$project->{ProfileReads}}) {
195 my $task = $profile->{TaskId};
196 my $genome = $profile->{Genome};
197 my $lib = $profile->{Library};
198 my $name = $profile->{Name};
200 $name =~ s/\s\(\)\./____/g;
202 my $outfile = "$lib.wig.gz $lib.profile.gif";
204 my $seqcheck = "if [ ! -e ~Data/Libraries/$lib.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for sequencing.\"; fi;";
205 my $cmds .= "$lib.wig.gz: ~Data/Libraries/$lib.txt\n";
206 $cmds .= "\t$seqcheck\n";
207 $cmds .= "\t".$PROFILEDIR.'/profile_reads_wig ~Data/Libraries/'.$lib.'.txt "'.$name.'" "'.$name.'" | gzip > '.$lib.'.wig.gz';
209 $cmds .= $lib.'.profile.gif: ~Data/Libraries/'.$lib.'.txt '.$root_dir.'/reference_data/'.$genome.'_tx_start_sites'."\n";
210 $cmds .= "\t".$PROFILEDIR.'/profile_reads_against_features $^ | '.$root_dir.'/scripts/profile_to_svg.pm | /opt/local/bin/convert - $@'."\n";
213 writeTask($profile, "ProfileReads", $outfile, $cmds);
220 sub WriteCompareLibTasks {
223 if(exists($project->{CompareLibraries})) {
224 for my $cmp (@{$project->{CompareLibraries}}) {
225 #Usage: /Users/ENCODE/EXPTRACK/QC/count_reads_in_peaks TF label1 fearure1 reads1 label2 feature2 reads2
227 my $task = $cmp->{TaskId};
229 print STDERR "Before: $tf\n"; $tf =~ s/\s/_/g; $tf =~ s/\(\)\./___/g; $tf =~ s/'//;
230 print STDERR "After: $tf\n";
231 my $genome = $cmp->{Genome};
232 my $features = "$root_dir/reference_data/".$genome."_upstream5k_downstream1k";
233 my $name1 = $cmp->{Library}->[0]->{Library};
234 my $name2 = $cmp->{Library}->[1]->{Library};
236 $name1 =~ s/\s\(\)\./____/g;
237 $name2 =~ s/\s\(\)\./____/g;
239 my $outfile = $name1."_".$name2.".compare ";
241 my $seqcheck = "if [ ! -e ~Data/Libraries/$name1.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for $name1 sequencing.\"; fi;";
242 $seqcheck .= "\n\tif [ ! -e ~Data/Libraries/$name2.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for $name2 sequencing.\"; fi;";
244 my $cmd = "$outfile: ~Data/Libraries/$name1.txt ~Data/Libraries/$name2.txt\n\t$seqcheck\n\t$root_dir/bin/count_reads_in_peaks $tf $name1 $features ~Data/Libraries/$name1.txt $name2 $features ~Data/Libraries/$name2.txt > \$@\n";
246 writeTask($cmp, "CompareLibraries", $outfile, $cmd);
254 sub WriteMotifFindingTasks {
257 if(exists($project->{MotifFinding})) {
258 for my $motiffind (@{$project->{MotifFinding}}) {
259 #my $caller = $motiffind->{Caller};
260 #my $name = $motiffind->{Genome};
261 #my $set = $motiffind->{Set};
262 #my $width = $motiffind->{Width};
263 #my $options = $motiffind->{Options};
266 #for my $i (0..scalar(@{$project->{PeakCalling}})-1) {
267 #if($project->{PeakCalling}->[$i]->{Name} eq $set && $project->{PeakCalling}->[$i]->{Caller} eq $caller) { $index = $i; }
270 #if($index < 0) { print STDERR "Error: Peak calling $set ($caller) not found in configuration.\n"; }
272 ##Usage: BioProspector -i fastafile
274 #my $outfile .= $set."-".$caller.".w".$width.".biop";
277 #if($caller eq "QuEST") { $input_file = "peak_caller.ChIP.out.fasta"; }
278 #elsif($caller eq "WingPeaks") { $input_file = "$set.peaks.fasta"; }
279 #else { print STDERR "Unable to compare peaks from $caller\n"; }
281 #my $set_fasta = $caller."_".$set."/".$input_file;
283 #my $cmd = "$outfile: $set_fasta\n";
284 #$cmd .= "\t$BIOP $options -W $width -i ".'$< -o $@'."\n";
286 #$file_list .= "$outfile ";
294 sub WritePeakCallingTasks {
297 if(exists($project->{PeakCalling})) {
298 for my $peakcall (@{$project->{PeakCalling}}) {
299 my $task = $peakcall->{TaskId};
300 my $name = $peakcall->{Name};
301 my $caller = $peakcall->{Caller};
303 my $signal = $peakcall->{Signal}->{Library};
304 my $bg = $peakcall->{Background}->{Library};
305 my $genome = $peakcall->{Genome};
307 my $seqcheck = "if [ ! -e ~Data/Libraries/$signal.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for $signal sequencing.\"; fi;";
308 $seqcheck .= "\n\tif [ ! -e ~Data/Libraries/$bg.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for $bg sequencing.\"; fi;";
316 if($caller eq "QuEST") {
317 $outfile .= "peak_caller.ChIP.out peak_caller.ChIP.out.bedgraph peak_caller.ChIP.out.fasta";
319 $cmd .= "peak_caller.ChIP.out: ~Data/Libraries/$signal.txt ~Data/Libraries/$bg.txt\n";
320 $cmd .= "\t$seqcheck\n";
321 $cmd .= "\trm -f background_RX_noIP.align.txt\n";
322 $cmd .= "\trm -f pseudo_ChIP_RX_noIP.align.txt\n";
323 $cmd .= "\t".$QUESTDIR.'/generate_QuEST_parameters.pl -rp '.$GENOMEDIR.'/QuEST_'.$genome.' -solexa_align_ChIP ~Data/Libraries/'.$signal.'.txt -solexa_align_RX_noIP ~Data/Libraries/'.$bg.'.txt -ap '.$data_dir.'/Tasks/'.$task.' -silent > '.$name.'.QuEST.log;'."\n";
324 $cmd .= "\t".$QUESTDIR.'/run_QuEST_with_param_file.pl -p QuEST.batch.pars >> '.$name.'.QuEST.log;'."\n";
325 $cmd .= "\t".'rm -rf scores;'."\n\n";
327 $cmd .= "peak_caller.ChIP.out.tab: peak_caller.ChIP.out\n";
328 $cmd .= "\t$root_dir/scripts/tabify_quest.sh \$< > \$@\n\n";
330 $cmd .= "peak_caller.ChIP.out.bedgraph: peak_caller.ChIP.out\n";
331 $cmd .= "\t$root_dir/scripts/QuEST_2_BED.pm \$< $name\n";
333 $cmd .= "peak_caller.ChIP.out.fasta: peak_caller.ChIP.out.tab\n";
334 $cmd .= "\t".'cat $< | '.$root_dir.'/scripts/extract_peaks.pm '.$root_dir.'/reference_data/hg18_chrom_list.txt > $@'."\n";
336 } elsif($caller eq "WingPeaks") {
337 $outfile .= "$name.peaks $name.peaks.fasta ";
339 $cmd .= "$name.peaks: ~Data/Libraries/$signal.txt ~Data/Libraries/$bg.txt\n";
340 $cmd .= "\t".$WINGPEAKSDIR.'/ChIPSeq_PeakCaller_ENCODE -gn '.$WINGPEAKSGENOMEDIR.'/'.$genome.'_chrlist.cod -it '.$name.' -in 1 -if ~Data/Libraries/'.$signal.'.txt -ct Background -cn 1 -cf ~Data/Libraries/'.$bg.'.txt -ot '.$name.' > '.$name.'.log;'."\n";
342 $cmd .= "%.peaks.tab: %.peaks\n";
343 $cmd .= "\t".'cat $< | awk \'{ print $$1"\t"$$2"\t"$$3"\t"$$4"\t"$$7}\' > $@'."\n\n";
345 $cmd .= "%.peaks.fasta: %.peaks.tab\n";
346 $cmd .= "\t".'cat $< | '.$root_dir.'/scripts/extract_peaks.pm '.$root_dir.'/reference_data/hg18_chrom_list.txt > $@'."\n\n";
347 } elsif($caller eq "MACS") {
348 $outfile .= $name."_peaks.bed ".$name."_peaks.fasta ";
350 $cmd .= $name."_peaks.bed: ~Data/Libraries/$signal.txt ~Data/Libraries/$bg.txt\n";
351 $cmd .= "\t$seqcheck\n";
352 $cmd .= "\tcat ~Data/Libraries/$signal.txt | $root_dir/scripts/align_to_bed.pm > $signal.bed\n";
353 $cmd .= "\tcat ~Data/Libraries/$bg.txt | $root_dir/scripts/align_to_bed.pm > $bg.bed\n";
354 $cmd .= "\t$MACSDIR/macs -t $signal.bed -c ./$bg.bed --name=$name --pvalue=1e-10 > $name.log 2> $name.err\n";
355 $cmd .= "\t".'echo "track name="'.$name.'" description="'.$name.'"GR_EtOH_Rep2_Peak_Calls" > header';
356 $cmd .= "\t".'cat header '.$name.'_peaks.bed > t; mv t '.$name.'_peaks.bed'."\n";
357 $cmd .= "\trm -f $signal.bed $bg.bed header\n";
358 $cmd .= "\t".'exit `grep -c "^CRITICAL" '.$name.'.err`'."\n\n";
360 $cmd .= "%_peaks.tab: %_peaks.bed\n";
361 $cmd .= "\t".'cat $< | awk \'{print NR"\t"$$1"\t"$$2"\t"$$3"\t1"}\' > $@'."\n\n";
363 $cmd .= "%_peaks.fasta: %_peaks.tab\n";
364 $cmd .= "\t".'cat $< | '.$root_dir.'/scripts/extract_peaks.pm '.$root_dir.'/reference_data/hg18_chrom_list.txt > $@'."\n\n";
367 writeTask($peakcall, "PeakCalling", $outfile, $cmd);
373 sub WriteComparePeakCallingTasks {
376 if(exists($project->{ComparePeakCalls})) {
377 for my $cmppeakcall (@{$project->{ComparePeakCalls}}) {
378 my $caller1 = $cmppeakcall->{Caller1};
379 my $caller2 = $cmppeakcall->{Caller2};
380 my $name = $cmppeakcall->{Genome};
381 my $set1 = $cmppeakcall->{Set1};
382 my $set2 = $cmppeakcall->{Set2};
387 #for my $i (0..scalar(@{$xml->{PeakCalling}})-1) {
388 #if($xml->{PeakCalling}->[$i]->{Name} eq $set1 && $xml->{PeakCalling}->[$i]->{Caller} eq $caller1) { $index1 = $i; }
389 #if($xml->{PeakCalling}->[$i]->{Name} eq $set2 && $xml->{PeakCalling}->[$i]->{Caller} eq $caller2) { $index2 = $i; }
392 #if($index1 < 0) { print STDERR "Error: Peak calling $set1 ($caller1) not found in configuration.\n"; }
393 #if($index2 < 0) { print STDERR "Error: Peak calling $set2 ($caller2) not found in configuration.\n"; }
394 #if($index1 >= 0 && $index2 >= 0) {
395 ##Usage: /Users/ENCODE/EXPTRACK/QC/count_reads_in_peaks TF label1 fearure1 reads1 label2 feature2 reads2
398 #my $outfile .= $set1."-".$caller1."_v_".$set2."-".$caller2.".compare";
401 #if($caller1 eq "QuEST") { $peak_file_1 = "peak_caller.ChIP.out.tab"; }
402 #elsif($caller1 eq "WingPeaks") { $peak_file_1 = "$set1.peaks.tab"; }
403 #else { print STDERR "Unable to compare peaks from $caller1\n"; }
406 #if($caller2 eq "QuEST") { $peak_file_2 = "peak_caller.ChIP.out.tab"; }
407 #elsif($caller2 eq "WingPeaks") { $peak_file_2 = "$set2.peaks.tab"; }
408 #else { print STDERR "Unable to compare peaks from $caller2\n"; }
410 #my $set1_feat = $caller1."_".$set1."/".$peak_file_1;
411 #my $set2_feat = $caller2."_".$set2."/".$peak_file_2;
413 #my $cmd = "$outfile: $set1_feat $set2_feat\n";
414 #$cmd .= "\t~/EXPTRACK/QC/count_reads_in_peaks NA $set1-$caller1 $set1_feat ~Data/Libraries/".$xml->{PeakCalling}->[$index1]->{Signal}->{Library}.".txt ";
415 #$cmd .= "$set2-$caller2 $set2_feat ~Data/Libraries/".$xml->{PeakCalling}->[$index2]->{Signal}->{Library}.".txt > \$@\n";
417 #$file_list .= "$outfile ";
427 my $dir= "$data_dir/Projects";
428 my $filename = "$dir/Projects.xml";
429 if(!defined($option)) { $option = ""; }
430 `$root_dir/scripts/analys_track_main.py getProjects $option $dir`;
431 `cat $filename | sed -e "s/\&/_and_/" > t; mv t $filename`;