Added wiggle fiole generation

diff --git a/htswanalysis/src/Makefile b/htswanalysis/src/Makefile
index 6563c50f771f8314c3ab37abd074981b2a7374cf..b52b347759bdfb0b8498be50218c72f75ef4eaa9 100644 (file)
--- a/htswanalysis/src/Makefile
+++ b/htswanalysis/src/Makefile
@@ -1,7 +1,7 @@
 CPPFLAGS=-g -Wall -O3 -I/opt/local/include
 LDFLAGS=-lgsl -lgslcblas -lm -L/opt/local/lib
 
-TARGETS=complexity_count qPCR profile_reads_against_features
+TARGETS=complexity_count qPCR profile_reads_against_features profile_reads_wig
 
 all: $(TARGETS)
 

diff --git a/htswanalysis/src/profile_reads_wig.cpp b/htswanalysis/src/profile_reads_wig.cpp
new file mode 100644 (file)
index 0000000..58be9dd
--- /dev/null
+++ b/htswanalysis/src/profile_reads_wig.cpp
@@ -0,0 +1,117 @@
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <queue>
+#include <math.h>
+
+#define WINDOW 25
+
+using namespace std;
+
+typedef pair<unsigned int, unsigned int> Hit;
+typedef vector<Hit> Hits;
+
+void print_wig_header(char* name, char* desc) {
+  //cout << "browser position chr1" << endl;
+  cout << "browser hide all" << endl;
+  cout << "browser full altGraph" << endl;
+  cout << "track type=wig name=" << name << " description=" << desc << " visibility=full autoScale=on color=0,0,0 priority=10" << endl;
+}
+
+bool compare_hits(const Hit &a, const Hit &b)
+{
+  if(a.first == b.first) {
+    return a.second < b.second;
+  } else {
+    return a.first < b.first;
+  }
+}
+
+int main(int argc, char** argv) {
+  if(argc != 4) { cerr << "Usage: " << argv[0] << " Aligned_reads Wiggle_Track_Name Wiggle_Track_Description" << endl; return(0); }
+  char filename[1024]; strcpy(filename,argv[1]);
+  Hits data;
+
+  char track_name[1024]; strcpy(track_name, argv[2]);
+  char track_desc[1024]; strcpy(track_desc, argv[3]);
+
+  print_wig_header(track_name, track_desc);
+  
+  ifstream seqs(filename);
+  while(seqs.peek() != EOF) {
+    char chr_str[32];
+    char pos_str[32];
+    char n_match_str[32];
+    char chr_start[4];
+
+    seqs.ignore(128,' ');
+    seqs.ignore(128,' ');
+    seqs.getline(n_match_str,32,' ');
+    n_match_str[seqs.gcount()] = '\0';
+    unsigned int nmatch = atoi(n_match_str);
+    if(nmatch != 1) { seqs.ignore(1024,'\n'); continue; }
+
+    seqs.get(chr_start,4);
+    chr_start[3] = '\0';
+    if(chr_start[0] != 'c' || chr_start[1] != 'h' || chr_start[2] != 'r') { seqs.ignore(2048,'\n'); continue; }
+    seqs.getline(chr_str,32,':');
+    chr_str[seqs.gcount()] = '\0';
+    seqs.getline(pos_str,32,' ');
+    pos_str[seqs.gcount()] = '\0';
+    seqs.ignore(1024,'\n');
+
+    unsigned int chr;
+
+    if(chr_str[0] == 'X') { chr = 23; } else
+    if(chr_str[0] == 'Y') { chr = 24; } else
+    if(chr_str[0] == 'M' || chr_str[0] == 'm') { chr = 25; } else
+    if(chr_str[0] == 'R') { continue; } else
+    { chr = atoi(chr_str); }
+
+    unsigned int pos = atoi(pos_str);
+
+    Hit hit(chr,pos); 
+    data.push_back(hit);
+  }
+  seqs.close(); 
+
+  std::sort(data.begin(),data.end(),compare_hits);
+
+  unsigned int n_aligned_reads = data.size();
+  cerr << n_aligned_reads/1000000.00 << "M reads\n";
+ 
+  queue<Hit> Q;
+  for(Hits::iterator i = data.begin(); i != data.end();) {
+    if(Q.empty() || Q.front().first != (*i).first) {
+      while(!Q.empty()) { Q.pop(); }
+      char chr[3];
+      switch(i->first) {
+        case 23: chr[0]='X'; chr[1]='\0'; break;
+        case 24: chr[0]='Y'; chr[1]='\0'; break;
+        case 25: chr[0]='M'; chr[1]='\0'; break;
+        case 26: chr[0]='R'; chr[1]='\0'; break;
+        default: chr[0]=(char)(48+(int)(i->first/10));
+                 chr[1]=(char)(48+(int)(i->first - 10*(int)(i->first/10)));
+                 if(chr[0]=='0') { chr[0]=chr[1]; chr[1]='\0'; }
+      }
+      chr[2] = '\0';
+      cout << "variableStep chrom=chr" << chr << endl;
+    }
+
+    unsigned int chrnum = i->first;
+    unsigned int posnum = i->second;
+    while(i != data.end() && i->first == chrnum && i->second == posnum) {
+      Q.push(*i); ++i; 
+    }
+
+    while(!Q.empty() && Q.front().second + WINDOW < Q.back().second) { 
+      Q.pop(); 
+    }
+
+    if(Q.size() > 1) {
+      //cout << Q.back().second << " " << log(Q.size()) << endl;
+      cout << Q.back().second << " " << Q.size() / (double)(n_aligned_reads/1000000.0) << endl;
+    }
+  }
+}
+