(no commit message)

diff --git a/htswanalysis/src/Methylseq/Methylseq_Analysis.cpp b/htswanalysis/src/Methylseq/Methylseq_Analysis.cpp
index ce47ea5df704d504f00b4cbc69c7a3d1db108ca0..2581bd5d533ba48c3fc988c59107b8b2a4be302f 100755 (executable)
--- a/htswanalysis/src/Methylseq/Methylseq_Analysis.cpp
+++ b/htswanalysis/src/Methylseq/Methylseq_Analysis.cpp
@@ -222,11 +222,12 @@ unsigned int count_regions_methylated(Regions& r, unsigned int msp_threshold, un
 unsigned int count_errors(Regions& r, unsigned int msp_threshold, unsigned int hpa_threshold);
 
 int main(int argc, char** argv) {
-  if(argc < 4) { cerr << "Usage: " << argv[0] << " msp1_reads hpa2_reads  regions\n"; exit(1); }
+  if(argc < 5) { cerr << "Usage: " << argv[0] << " msp1_reads hpa2_reads  regions name\n"; exit(1); }
 
   char msp1_filename[1024]; strcpy(msp1_filename,argv[1]);
   char hpa2_filename[1024]; strcpy(hpa2_filename,argv[2]);
   char region_filename[1024]; strcpy(region_filename,argv[3]);
+  char name[1024]; strcpy(name,argv[4]);
 
   Reads msp1; read_align_file(msp1_filename, msp1);
   Reads hpa2; read_align_file(hpa2_filename, hpa2);
@@ -239,33 +240,36 @@ int main(int argc, char** argv) {
   pair<unsigned int,unsigned int> total = count_region_tags(regions);
   cout << total.first << " MSP1 tags, " << total.second << " HPA2 tags" << endl;
 
-  unsigned int n_assayed = count_regions_assayed(regions);
+  unsigned int n_assayed = count_regions_assayed(regions,2);
   cout << n_assayed << " of " << regions.size() << " regions assayed (" << (double)n_assayed/(double)regions.size() << ")\n";
 
-  unsigned int n_methylated = count_regions_methylated(regions);
+  unsigned int n_methylated = count_regions_methylated(regions,2);
   cout << n_methylated << " of " << n_assayed << " regions methylated (" << (double)n_methylated/(double)n_assayed << ")\n";
 
-  unsigned int n_error = count_errors(regions);
+  unsigned int n_error = count_errors(regions,2);
   cout << n_error << " of " << regions.size() << " regions with error (" << (double)n_error/(double)regions.size() << ")\n";
 
 
-  ofstream assayed("Assayed.bed");
-  ofstream methylated("Methylated.bed");
+  ofstream methylseqbed("MethylseqCalls.bed");
   ofstream methylcalls("MethylseqCalls.tab");
 
-  unsigned int msp_threshold = 1;
-  unsigned int hpa_threshold = 1;
+  methylseqbed << "track name=\"" << name << "\" description=\"" << name << "\" visibility=1 itemRgb=\"On\"" << endl;
+ 
+
+  unsigned int msp_threshold = 2;
+  unsigned int hpa_threshold = 2;
 
   for(unsigned int idx = 0; idx < regions.size(); idx++) {
     bool is_assayed = 0;
     bool is_methylated = 0;
     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) {
-        methylated << regions[idx].chr << "\t" << regions[idx].start << "\t" << regions[idx].end << endl;
+        methylseqbed << regions[idx].chr << "\t" << regions[idx].start << "\t" << regions[idx].end << "\thpa2\t" << regions[idx].hpa2_count << "\t+\t" << regions[idx].start << "\t" << regions[idx].end << "\t255,140,0" << endl;
         is_methylated = 1;
-      } 
+      } else { 
+        methylseqbed << regions[idx].chr << "\t" << regions[idx].start << "\t" << regions[idx].end << "\tmsp1\t" << regions[idx].msp1_count << "\t+\t" << regions[idx].start << "\t" << regions[idx].end << "\t0,0,205" << endl;
+      }
     }
     methylcalls << regions[idx].name << "\t" << regions[idx].chr << "\t" << regions[idx].start << "\t" << regions[idx].end << "\t" << regions[idx].msp1_count << "\t" << regions[idx].hpa2_count << "\t" << is_assayed << "\t" << is_methylated <<  endl;
   }