Bug fixes to bring up methylseq analysis Tasks
authorTim Reddy Tim <treddy@hudsonalpha.org>
Mon, 8 Dec 2008 14:40:24 +0000 (14:40 +0000)
committerTim Reddy Tim <treddy@hudsonalpha.org>
Mon, 8 Dec 2008 14:40:24 +0000 (14:40 +0000)
htswanalysis/scripts/ConfigureTasks.pm
htswanalysis/src/Methylseq/Methylseq_Analysis.cpp

index 869af453954d623c046d8eed9f8d27ff31500320..8fe851ea12611df2f123b8d9cf8d753aea6e927e 100755 (executable)
@@ -52,12 +52,13 @@ my $SUMMARIZE_PROJECT="$root_dir/scripts/SummarizeProject2.pm";
 ###
 ### Important: get the xml file form the server.
 ###
-my $projects_file = getProjectsXML($parm);
+#my $projects_file = getProjectsXML($parm);
+my $projects_file = "/Users/Data/Projects.xml";
 
 my $xmldoc = XML::Simple->new();
 
 my $xmlfile = $projects_file;
-my $xml = $xmldoc->XMLin($xmlfile, ForceArray => ['Project','MethylSeq', 'ComparePeakCalls','CompareLibraries','PeakCalling','ProfileReads','qPCR','MotifFinding'], KeepRoot=>1);
+my $xml = $xmldoc->XMLin($xmlfile, ForceArray => ['Project','Methylseq', 'ComparePeakCalls','CompareLibraries','PeakCalling','ProfileReads','qPCR','MotifFinding'], KeepRoot=>1);
 
 my $projects = "";
 
@@ -100,7 +101,7 @@ close(MAKE);
 
 print STDERR "Project makefile compete. Starting to build tasks\n";
 
-`cd $data_dir/Projects && make -j 4 -f TaskMakefile > make_tasks.log 2> make_tasks.err`;
+#`cd $data_dir/Projects && make -j 4 -f TaskMakefile > make_tasks.log 2> make_tasks.err`;
 
 sub writeProject {
   my $project = shift;
@@ -181,12 +182,12 @@ sub WriteMethylseqTasks {
 
       my $outfile = "methylseq_summary";
   
-      my $seqcheck  = "if [ ! -e $data_dir/Libraries/$msp1.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for sequencing.\"; fi;"; 
-         $seqcheck .= "if [ ! -e $data_dir/Libraries/$hpa2.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for sequencing.\"; fi;"; 
-         $seqcheck .= "if [ ! -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;"; 
+      my $seqcheck  = "if [ ! -e $data_dir/Libraries/$msp1.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for sequencing.\"; fi;\n"; 
+         $seqcheck .= "\tif [ ! -e $data_dir/Libraries/$hpa2.txt ]; then $root_dir/scripts/analys_track_main.py updsts $task \"Waiting for sequencing.\"; fi;\n"; 
+         $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"; 
 
       my $cmd  = "$outfile: $data_dir/Libraries/$msp1.txt $data_dir/Libraries/$hpa2.txt\n";
-         $cmd .= "\t$METHYLSEQDIR/Methylseq_Analysis $data_dir/Libraries/$msp1.txt $root_dir/reference_data/".$genome."_methylseq_regions $data_dir/Libraries/$hpa2.txt > \$@\n";
+         $cmd .= "\t$METHYLSEQDIR/Methylseq_Analysis $data_dir/Libraries/$msp1.txt $data_dir/Libraries/$hpa2.txt $root_dir/reference_data/".$genome."_methylseq_regions > \$@\n";
       writeTask($mseq, $root, $outfile, $cmd, $seqcheck);
 
       $tasks .= $task." ";
@@ -195,10 +196,6 @@ sub WriteMethylseqTasks {
   return $tasks;
 }
 
-sub WriteProfileTasks {
-
-} 
-
 sub WriteQPCRTasks {
   my $project = shift;
   my $tasks = "";
index 14a1248d6f1ab3209e4978738bfaf9ea17d59f7c..ce47ea5df704d504f00b4cbc69c7a3d1db108ca0 100755 (executable)
@@ -170,7 +170,7 @@ unsigned int count_regions_assayed(Regions& r, unsigned int threshold = 1) {
 unsigned int count_regions_methylated(Regions& r, unsigned int msp_threshold = 1, unsigned int hpa_threshold = 1) {
   unsigned int n_methylated = 0;
   for(Regions::iterator i = r.begin(); i != r.end(); ++i) {
-    if(i->msp1_count >= msp_threshold && i->hpa2_count >= hpa_threshold) { n_methylated++; }
+    if(i->msp1_count >= msp_threshold && i->hpa2_count < hpa_threshold) { n_methylated++; }
   }
   return n_methylated;
 }
@@ -262,7 +262,7 @@ int main(int argc, char** argv) {
     if(regions[idx].msp1_count >= msp_threshold) {
       assayed << regions[idx].chr << "\t" << regions[idx].start << "\t" << regions[idx].end << endl;
       is_assayed = 1;
-      if(regions[idx].hpa2_count >= hpa_threshold) {
+      if(regions[idx].hpa2_count < hpa_threshold) {
         methylated << regions[idx].chr << "\t" << regions[idx].start << "\t" << regions[idx].end << endl;
         is_methylated = 1;
       }