htswanalysis/src/profile_reads_wig.cpp

   1 #include <iostream>
   2 #include <fstream>
   3 #include <vector>
   4 #include <queue>
   5 #include <math.h>
   6
   7 #define WINDOW 25
   8
   9 using namespace std;
  10
  11 class Read {
  12   public:
  13     string chr;
  14     int pos;
  15     bool strand;
  16
  17     Read(string chr, int pos, bool strand) { this->chr = chr; this->pos = pos; this->strand = strand; }
  18     Read(const Read& r) { this->chr = r.chr; this->pos = r.pos; this->strand = r.strand; }
  19     Read& operator=(const Read& r) {this->chr = r.chr; this->pos = r.pos; this->strand = r.strand; return *this;}
  20 };
  21
  22 typedef vector<Read> Reads;
  23
  24 void print_wig_header(char* name, char* desc) {
  25   //cout << "browser position chr1" << endl;
  26   cout << "browser hide all" << endl;
  27   cout << "browser full altGraph" << endl;
  28   cout << "track type=wig name=" << name << " description=" << desc << " visibility=full autoScale=on color=0,0,0 priority=10" << endl;
  29 }
  30
  31 bool compare_reads(const Read &a, const Read &b) { if(a.chr == b.chr) { return a.pos < b.pos; } else { return a.chr < b.chr; } }
  32
  33 void split (const string& text, const string& separators, vector<string>& words);
  34 void read_align_file(char* filename, Reads& features);
  35
  36 int main(int argc, char** argv) {
  37   if(argc != 4) { cerr << "Usage: " << argv[0] << " Aligned_reads Wiggle_Track_Name Wiggle_Track_Description" << endl; return(0); }
  38   char filename[1024]; strcpy(filename,argv[1]);
  39   Reads data;
  40
  41   char track_name[1024]; strcpy(track_name, argv[2]);
  42   char track_desc[1024]; strcpy(track_desc, argv[3]);
  43
  44   print_wig_header(track_name, track_desc);
  45   read_align_file(filename, data);
  46
  47   unsigned int n_aligned_reads = data.size();
  48   cerr << n_aligned_reads/1000000.00 << "M reads\n";
  49
  50   queue<Read> Q;
  51   for(Reads::iterator i = data.begin(); i != data.end();) {
  52     if(Q.empty() || Q.front().chr != (*i).chr) {
  53       while(!Q.empty()) { Q.pop(); }
  54       cout << "variableStep chrom=" << i->chr.c_str() << endl;
  55     }
  56
  57     string chr_string = i->chr;
  58     int posnum = i->pos;
  59     while(i != data.end() && i->chr == chr_string && i->pos == posnum) {
  60       Q.push(*i); ++i;
  61     }
  62
  63     while(!Q.empty() && Q.front().pos + WINDOW < Q.back().pos) {
  64       Q.pop();
  65     }
  66
  67     if(Q.size() > 1) {
  68       cout << Q.back().pos << " " << Q.size() / (double)(n_aligned_reads/1000000.0) << endl;
  69     }
  70   }
  71 }
  72
  73 void read_align_file(char* filename, Reads& features) {
  74   string delim(" \n");
  75   string location_delim(":");
  76   char strand_str[2]; strand_str[1] = '\0';
  77   ifstream seqs(filename);
  78   string name("");
  79   while(seqs.peek() != EOF) {
  80     char line[2048];
  81     seqs.getline(line,2048,'\n');
  82
  83     string line_str(line);
  84     vector<string> fields;
  85     split(line_str, delim, fields);
  86     if(fields.size() != 7) { continue; }
  87
  88     vector<string> location; split(fields[3], location_delim, location);
  89     string chr = location[0];
  90     if(chr == "NA") { continue; }
  91     int pos = atoi(location[1].c_str());
  92     bool strand = ((fields[4].c_str())[0] == 'F')?0:1;
  93
  94     if(strand == 0) {
  95       Read read(chr,pos,false);
  96       features.push_back(read);
  97     } else {
  98       Read read(chr,pos+25,true);
  99       features.push_back(read);
 100     }
 101   }
 102   seqs.close();
 103
 104   //sort the data so we can run through it once
 105   std::sort(features.begin(),features.end(),compare_reads);
 106   cerr << "Reads sorted\n";
 107 }
 108
 109 void split (const string& text, const string& separators, vector<string>& words) {
 110
 111     size_t n     = text.length ();
 112     size_t start = text.find_first_not_of (separators);
 113
 114     while (start < n) {
 115         size_t stop = text.find_first_of (separators, start);
 116         if (stop > n) stop = n;
 117         words.push_back (text.substr (start, stop-start));
 118         start = text.find_first_not_of (separators, stop+1);
 119     }
 120 }