htswanalysis/scripts/ConfigureTasks.pm

   1 #!/usr/bin/perl -w
   2
   3 use XML::Simple;
   4 use Cwd;
   5 use Config;
   6
   7 use threads;
   8 use threads::shared;
   9 use Thread::Queue;
  10
  11 use strict;
  12 use warnings;
  13
  14 #
  15 # Number of threads to run simultaneously.
  16 #
  17 my $Nthread = 1;
  18 my $CmdQueue = Thread::Queue->new;
  19
  20 my $root_dir = shift;
  21 my $data_dir = shift;
  22 my $parm = shift;
  23
  24 my %libs;
  25
  26 my $BIOP = "$root_dir/bin/BioProspector.mac";
  27 my $QUESTDIR = "$root_dir/bin/QuEST";
  28 my $MACSDIR  = "$root_dir/bin";
  29 my $WINGPEAKSDIR = "$root_dir/bin";
  30 my $WINGPEAKSGENOMEDIR = "$root_dir/reference_data";
  31 my $PROFILEDIR = "$root_dir/bin";
  32 my $GENOMEDIR = "/Volumes/Genomes";
  33 my $QPCRDIR = "$root_dir/bin";
  34 my $QPCRTESTDIR = "$root_dir/reference_data/qPCR_Tests";
  35 my $QPCRBACKGROUND = "$root_dir/reference_data/GenericBackground";
  36
  37 my $SUMMARIZE_PROJECT="$root_dir/scripts/SummarizeProject2.pm";
  38
  39 ###
  40 ### Check directory for project and task folders. Make them if necessary
  41 ### May not have permissions for this. These should be moved to an install script
  42 `if [ ! -e $data_dir/Tasks ]; then mkdir $data_dir/Tasks; fi;`;
  43 `if [ ! -e $data_dir/Projects ]; then mkdir $data_dir/Projects; fi;`;
  44
  45
  46 ###
  47 ### Important: get the xml file form the server.
  48 ###
  49 my $projects_file = getProjectsXML($parm);
  50
  51 my $xmldoc = XML::Simple->new();
  52
  53 my $xmlfile = $projects_file;
  54 my $xml = $xmldoc->XMLin($xmlfile, ForceArray => ['Project','ComparePeakCalls','CompareLibraries','PeakCalling','ProfileReads','qPCR','MotifFinding'], KeepRoot=>1);
  55
  56 my $projects = "";
  57
  58 my @tasks;
  59
  60 if(exists($xml->{Projects}->{Project})) {
  61   for my $project (@{$xml->{Projects}->{Project}}) {
  62     push @tasks, @{writeProject($project)};
  63
  64     my $projectname = $project->{Name};
  65     my $projectId   = $project->{ProjectId};
  66
  67     $projects .= "$projectId ";
  68   }
  69 }
  70
  71 print STDERR "Projects all complete\n";
  72
  73 my %saw; my $unique_tasks = ""; my $unique_task_dirs = "";
  74 for(grep(!$saw{$_}++, @tasks)) {
  75   $unique_tasks .= "$_ ";
  76   $unique_task_dirs .= "../Tasks/$_ ";
  77 }
  78
  79 print STDERR "writing tasks\n";
  80
  81 open(MAKE,">$data_dir/Projects/TaskMakefile");
  82 print MAKE "all: $unique_task_dirs\n\n.PHONY: $unique_task_dirs\n\n";
  83 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"; }
  84 close(MAKE);
  85
  86 print STDERR "Task makefile complete\n";
  87
  88 my $index_list = "";
  89 for(split(/ /,$projects)) { $index_list .= "$_/index.html "; }
  90
  91 open(MAKE,">$data_dir/Projects/ProjectMakefile");
  92 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 ../../LibraryInfo.xml > `basename $@`'."\n";
  93 close(MAKE);
  94
  95 print STDERR "Project makefile compete. Starting to build tasks\n";
  96
  97 `make -j 4 -f TaskMakefile > make_tasks.log 2> make_tasks.err`;
  98
  99 sub writeProject {
 100   my $project = shift;
 101   my $projectname = $project->{Name}; $projectname =~ s/\s/_/g;
 102   my $projectid = $project->{ProjectId};
 103
 104   my $projectdir = $projectid;
 105
 106   my $tasks = "";
 107   $tasks .= WriteQPCRTasks($project);
 108   $tasks .= WriteProfileTasks($project);
 109   $tasks .= WriteCompareLibTasks($project);
 110   $tasks .= WriteMotifFindingTasks($project);
 111   $tasks .= WritePeakCallingTasks($project);
 112
 113   print STDERR "Wrote tasks for project '$projectname' with tasks ($tasks)\n";
 114
 115   my @tasks = split(/ /,$tasks);
 116
 117   if(! (-e $projectdir && -d $projectdir)) {
 118     `mkdir $data_dir/Projects/$projectdir`;
 119   }
 120
 121   print STDERR "Created directory $projectdir\n";
 122
 123   $xmldoc->XMLout($project, OutputFile=>"$data_dir/Projects/$projectid/Project.xml", RootName=>"Project", XMLDecl=>1);
 124
 125   print STDERR "Wrote configureation for project $projectdir\n";
 126
 127   return \@tasks;
 128 }
 129
 130 sub registerTask {
 131   my $taskid = shift;
 132   `$root_dir/scripts/analys_track_main.py updsts $taskid`;
 133   print STDERR "Task $taskid registered.\n";
 134 }
 135
 136 sub writeTask {
 137   my $task = shift;
 138   my $root = shift;
 139   my $outfile = shift;
 140   my $cmd = shift;
 141
 142   my $taskid = $task->{TaskId};
 143   my $taskdir = "$data_dir/Tasks/".$taskid;
 144
 145   if(! (-e $taskdir && -d $taskdir)) {
 146     `mkdir $taskdir`;
 147   }
 148
 149   $xmldoc->XMLout($task, OutputFile=>"$taskdir/Task.xml", RootName=>$root, XMLDecl=>1);
 150
 151   open(MAKEFILE, ">$taskdir/Makefile");
 152   print MAKEFILE "all: .notify $outfile | .start\n\n.PHONY: .notify .start\n\n";
 153   print MAKEFILE ".start:\n\t$root_dir/scripts/analys_track_main.py updsts $taskid \"Processing\"\n\ttouch .start\n\n";
 154
 155   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";
 156   print MAKEFILE "$cmd";
 157   close(MAKEFILE);
 158   registerTask($taskid);
 159 }
 160
 161 sub WriteQPCRTasks {
 162   my $project = shift;
 163   my $tasks = "";
 164   if(exists($project->{qPCR})) {
 165     for my $qpcr (@{$project->{qPCR}}) {
 166       my $task   = $qpcr->{TaskId};
 167       my $lib    = $qpcr->{Library};
 168       my $genome = $qpcr->{Genome};
 169       my $name   = $qpcr->{Name}; $name =~ tr/ \(\)\./____/;
 170
 171       my $background = $QPCRBACKGROUND;
 172       my $testdir =  $QPCRTESTDIR;
 173       my $outfile = "$name.qPCR";
 174
 175       my $seqcheck = "if [ ! -e ~Data/Libraries/$lib.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for sequencing.\"; fi;";
 176       my $cmd = "$outfile: ~Data/Libraries/$lib.txt\n\t$seqcheck\n\t$QPCRDIR/qPCR \$< $background $testdir > \$@\n";
 177       writeTask($qpcr, "qPCR", $outfile, $cmd);
 178
 179       $tasks .= $task." ";
 180     }
 181   }
 182   return $tasks;
 183 }
 184
 185 sub WriteProfileTasks {
 186   my $project = shift;
 187   my $tasks = "";
 188   if(exists($project->{ProfileReads})) {
 189     for my $profile (@{$project->{ProfileReads}}) {
 190       my $task   = $profile->{TaskId};
 191       my $genome = $profile->{Genome};
 192       my $lib =    $profile->{Library};
 193       my $name =   $profile->{Name};
 194       $name =~ tr/ /_/;
 195       $name =~ s/\s\(\)\./____/g;
 196
 197       my $outfile = "$lib.wig.gz $lib.profile.gif";
 198
 199       my $seqcheck = "if [ ! -e ~Data/Libraries/$lib.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for sequencing.\"; fi;";
 200       my $cmds .= "$lib.wig.gz: ~Data/Libraries/$lib.txt\n";
 201       $cmds .= "\t$seqcheck\n";
 202       $cmds .= "\t".$PROFILEDIR.'/profile_reads_wig ~Data/Libraries/'.$lib.'.txt "'.$name.'" "'.$name.'" | gzip > '.$lib.'.wig.gz';
 203       $cmds .= "\n\n";
 204       $cmds .= $lib.'.profile.gif: ~Data/Libraries/'.$lib.'.txt '.$root_dir.'/reference_data/'.$genome.'_tx_start_sites'."\n";
 205       $cmds .= "\t".$PROFILEDIR.'/profile_reads_against_features $^ | '.$root_dir.'/scripts/profile_to_svg.pm | /opt/local/bin/convert - $@'."\n";
 206       $cmds .= "\n";
 207
 208       writeTask($profile, "ProfileReads",  $outfile, $cmds);
 209       $tasks .= $task." ";
 210     }
 211   }
 212   return $tasks;
 213 }
 214
 215 sub WriteCompareLibTasks {
 216   my $project = shift;
 217   my $tasks = "";
 218   if(exists($project->{CompareLibraries})) {
 219     for my $cmp (@{$project->{CompareLibraries}}) {
 220       #Usage: /Users/ENCODE/EXPTRACK/QC/count_reads_in_peaks TF label1 fearure1 reads1 label2 feature2 reads2
 221
 222       my $task = $cmp->{TaskId};
 223       my $tf = $cmp->{TF};
 224       print STDERR "Before: $tf\n"; $tf =~ s/\s/_/g; $tf =~ s/\(\)\./___/g; $tf =~ s/'//;
 225       print STDERR "After: $tf\n";
 226       my $genome = $cmp->{Genome};
 227       my $features = "$root_dir/reference_data/".$genome."_upstream5k_downstream1k";
 228       my $name1 = $cmp->{Library}->[0]->{Library};
 229       my $name2 = $cmp->{Library}->[1]->{Library};
 230
 231       $name1 =~ s/\s\(\)\./____/g;
 232       $name2 =~ s/\s\(\)\./____/g;
 233
 234       my $outfile = $name1."_".$name2.".compare ";
 235
 236       my $seqcheck = "if [ ! -e ~Data/Libraries/$name1.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for $name1 sequencing.\"; fi;";
 237       $seqcheck .= "\n\tif [ ! -e ~Data/Libraries/$name2.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for $name2 sequencing.\"; fi;";
 238
 239       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";
 240
 241       writeTask($cmp, "CompareLibraries", $outfile, $cmd);
 242       $tasks .= $task." ";
 243     }
 244   }
 245   return $tasks;
 246 }
 247
 248
 249 sub WriteMotifFindingTasks {
 250   my $project = shift;
 251   my $tasks = "";
 252   if(exists($project->{MotifFinding})) {
 253     for my $motiffind (@{$project->{MotifFinding}}) {
 254       #my $caller =  $motiffind->{Caller};
 255       #my $name =    $motiffind->{Genome};
 256       #my $set =     $motiffind->{Set};
 257       #my $width =   $motiffind->{Width};
 258       #my $options = $motiffind->{Options};
 259 #
 260       #my $index = -1;
 261       #for my $i (0..scalar(@{$project->{PeakCalling}})-1) {
 262         #if($project->{PeakCalling}->[$i]->{Name} eq $set && $project->{PeakCalling}->[$i]->{Caller} eq $caller) { $index = $i; }
 263       #}
 264 #
 265       #if($index < 0) { print STDERR "Error: Peak calling $set ($caller) not found in configuration.\n"; }
 266       #if($index >= 0) {
 267         ##Usage: BioProspector -i fastafile
 268         #$set =~ s/ /_/g;
 269         #my $outfile .= $set."-".$caller.".w".$width.".biop";
 270   #
 271         #my $input_file;
 272         #if($caller eq "QuEST") { $input_file = "peak_caller.ChIP.out.fasta"; }
 273         #elsif($caller eq "WingPeaks") { $input_file = "$set.peaks.fasta"; }
 274         #else { print STDERR "Unable to compare peaks from $caller\n"; }
 275   #
 276         #my $set_fasta = $caller."_".$set."/".$input_file;
 277     #
 278         #my $cmd = "$outfile: $set_fasta\n";
 279         #$cmd .= "\t$BIOP $options -W $width -i ".'$< -o $@'."\n";
 280   #
 281         #$file_list .= "$outfile ";
 282         #$cmds .= $cmd."\n";
 283       #}
 284     }
 285   }
 286   return $tasks;
 287 }
 288
 289 sub WritePeakCallingTasks {
 290   my $project = shift;
 291   my $tasks = "";
 292   if(exists($project->{PeakCalling})) {
 293     for my $peakcall (@{$project->{PeakCalling}}) {
 294       my $task =   $peakcall->{TaskId};
 295       my $name =   $peakcall->{Name};
 296       my $caller = $peakcall->{Caller};
 297
 298       my $signal = $peakcall->{Signal}->{Library};
 299       my $bg =     $peakcall->{Background}->{Library};
 300       my $genome = $peakcall->{Genome};
 301
 302       my $seqcheck = "if [ ! -e ~Data/Libraries/$signal.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for $signal sequencing.\"; fi;";
 303       $seqcheck .= "\n\tif [ ! -e ~Data/Libraries/$bg.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for $bg sequencing.\"; fi;";
 304
 305       $tasks .= $task." ";
 306
 307       $name =~ s/\s/_/g;
 308
 309       my $outfile;
 310       my $cmd = "";
 311       if($caller eq "QuEST") {
 312         $outfile .= "peak_caller.ChIP.out peak_caller.ChIP.out.bedgraph peak_caller.ChIP.out.fasta";
 313
 314         $cmd   .= "peak_caller.ChIP.out: ~Data/Libraries/$signal.txt ~Data/Libraries/$bg.txt\n";
 315         $cmd   .= "\t$seqcheck\n";
 316         $cmd   .= "\trm -f background_RX_noIP.align.txt\n";
 317         $cmd   .= "\trm -f pseudo_ChIP_RX_noIP.align.txt\n";
 318         $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";
 319         $cmd   .= "\t".$QUESTDIR.'/run_QuEST_with_param_file.pl -p QuEST.batch.pars >> '.$name.'.QuEST.log;'."\n";
 320         $cmd   .= "\t".'rm -rf scores;'."\n\n";
 321
 322         $cmd  .= "peak_caller.ChIP.out.tab: peak_caller.ChIP.out\n";
 323         $cmd .= "\t$root_dir/scripts/tabify_quest.sh \$< > \$@\n\n";
 324
 325         $cmd .= "peak_caller.ChIP.out.bedgraph: peak_caller.ChIP.out\n";
 326         $cmd .= "\t$root_dir/scripts/QuEST_2_BED.pm \$< $name\n";
 327
 328         $cmd .= "peak_caller.ChIP.out.fasta: peak_caller.ChIP.out.tab\n";
 329         $cmd .= "\t".'cat $< | '.$root_dir.'/scripts/extract_peaks.pm '.$root_dir.'/reference_data/hg18_chrom_list.txt > $@'."\n";
 330
 331       } elsif($caller eq "WingPeaks") {
 332         $outfile .= "$name.peaks $name.peaks.fasta ";
 333
 334         $cmd .= "$name.peaks: ~Data/Libraries/$signal.txt ~Data/Libraries/$bg.txt\n";
 335         $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";
 336
 337         $cmd .= "%.peaks.tab: %.peaks\n";
 338         $cmd .= "\t".'cat $< | awk \'{ print $$1"\t"$$2"\t"$$3"\t"$$4"\t"$$7}\' > $@'."\n\n";
 339
 340         $cmd .= "%.peaks.fasta: %.peaks.tab\n";
 341         $cmd .= "\t".'cat $< | '.$root_dir.'/scripts/extract_peaks.pm '.$root_dir.'/reference_data/hg18_chrom_list.txt > $@'."\n\n";
 342       } elsif($caller eq "MACS") {
 343         $outfile .= $name."_peaks.bed ".$name."_peaks.fasta ";
 344
 345         $cmd .= $name."_peaks.bed: ~Data/Libraries/$signal.txt ~Data/Libraries/$bg.txt\n";
 346         $cmd .= "\t$seqcheck\n";
 347         $cmd .= "\tcat ~Data/Libraries/$signal.txt | $root_dir/scripts/align_to_bed.pm > $signal.bed\n";
 348         $cmd .= "\tcat ~Data/Libraries/$bg.txt | $root_dir/scripts/align_to_bed.pm > $bg.bed\n";
 349         $cmd .= "\t$MACSDIR/macs -t $signal.bed -c ./$bg.bed --name=$name --pvalue=1e-10 > $name.log 2> $name.err\n";
 350         $cmd .= "\t".'echo "track name="'.$name.'" description="'.$name.'"GR_EtOH_Rep2_Peak_Calls" > header';
 351         $cmd .= "\t".'cat header '.$name.'_peaks.bed > t; mv t '.$name.'_peaks.bed'."\n";
 352         $cmd .= "\trm -f $signal.bed $bg.bed header\n";
 353         $cmd .= "\t".'exit `grep -c "^CRITICAL" '.$name.'.err`'."\n\n";
 354
 355         $cmd .= "%_peaks.tab: %_peaks.bed\n";
 356         $cmd .= "\t".'cat $< | awk \'{print NR"\t"$$1"\t"$$2"\t"$$3"\t1"}\' > $@'."\n\n";
 357
 358         $cmd .= "%_peaks.fasta: %_peaks.tab\n";
 359         $cmd .= "\t".'cat $< | '.$root_dir.'/scripts/extract_peaks.pm '.$root_dir.'/reference_data/hg18_chrom_list.txt > $@'."\n\n";
 360       }
 361
 362       writeTask($peakcall, "PeakCalling", $outfile, $cmd);
 363     }
 364   }
 365   return $tasks;
 366 }
 367
 368 sub WriteComparePeakCallingTasks {
 369   my $project = shift;
 370   my $tasks = "";
 371   if(exists($project->{ComparePeakCalls})) {
 372     for my $cmppeakcall (@{$project->{ComparePeakCalls}}) {
 373       my $caller1 = $cmppeakcall->{Caller1};
 374       my $caller2 = $cmppeakcall->{Caller2};
 375       my $name = $cmppeakcall->{Genome};
 376       my $set1 = $cmppeakcall->{Set1};
 377       my $set2 = $cmppeakcall->{Set2};
 378
 379       my $index1 = -1;
 380       my $index2 = -1;
 381
 382       #for my $i (0..scalar(@{$xml->{PeakCalling}})-1) {
 383         #if($xml->{PeakCalling}->[$i]->{Name} eq $set1 && $xml->{PeakCalling}->[$i]->{Caller} eq $caller1) { $index1 = $i; }
 384         #if($xml->{PeakCalling}->[$i]->{Name} eq $set2 && $xml->{PeakCalling}->[$i]->{Caller} eq $caller2) { $index2 = $i; }
 385       #}
 386 #
 387       #if($index1 < 0) { print STDERR "Error: Peak calling $set1 ($caller1) not found in configuration.\n"; }
 388       #if($index2 < 0) { print STDERR "Error: Peak calling $set2 ($caller2) not found in configuration.\n"; }
 389       #if($index1 >= 0 && $index2 >= 0) {
 390         ##Usage: /Users/ENCODE/EXPTRACK/QC/count_reads_in_peaks TF label1 fearure1 reads1 label2 feature2 reads2
 391         #$set1 =~ s/ /_/g;
 392         #$set2 =~ s/ /_/g;
 393         #my $outfile .= $set1."-".$caller1."_v_".$set2."-".$caller2.".compare";
 394 #
 395         #my $peak_file_1;
 396         #if($caller1 eq "QuEST") { $peak_file_1 = "peak_caller.ChIP.out.tab"; }
 397         #elsif($caller1 eq "WingPeaks") { $peak_file_1 = "$set1.peaks.tab"; }
 398         #else { print STDERR "Unable to compare peaks from $caller1\n"; }
 399 #
 400         #my $peak_file_2;
 401         #if($caller2 eq "QuEST") { $peak_file_2 = "peak_caller.ChIP.out.tab"; }
 402         #elsif($caller2 eq "WingPeaks") { $peak_file_2 = "$set2.peaks.tab"; }
 403         #else { print STDERR "Unable to compare peaks from $caller2\n"; }
 404 #
 405         #my $set1_feat = $caller1."_".$set1."/".$peak_file_1;
 406         #my $set2_feat = $caller2."_".$set2."/".$peak_file_2;
 407   #
 408         #my $cmd = "$outfile: $set1_feat $set2_feat\n";
 409         #$cmd .= "\t~/EXPTRACK/QC/count_reads_in_peaks NA $set1-$caller1 $set1_feat ~Data/Libraries/".$xml->{PeakCalling}->[$index1]->{Signal}->{Library}.".txt ";
 410         #$cmd .= "$set2-$caller2 $set2_feat ~Data/Libraries/".$xml->{PeakCalling}->[$index2]->{Signal}->{Library}.".txt > \$@\n";
 411   #
 412         #$file_list .= "$outfile ";
 413         #$cmds .= $cmd."\n";
 414       #}
 415     }
 416   }
 417   return $tasks;
 418 }
 419
 420 sub getProjectsXML {
 421   my $option = shift;
 422   my $dir= "$data_dir/Projects";
 423   my $filename = "$dir/Projects.xml";
 424   if(!defined($option)) { $option = ""; }
 425   `$root_dir/scripts/analys_track_main.py getProjects $option $dir`;
 426   `cat $filename | sed -e "s/\&/_and_/" > t; mv t $filename`;
 427   return $filename;
 428 }