###
### 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 = "";
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;
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." ";
return $tasks;
}
-sub WriteProfileTasks {
-
-}
-
sub WriteQPCRTasks {
my $project = shift;
my $tasks = "";
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;
}
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;
}