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