htswanalysis/src/profile_reads_against_features.cpp

   1 #include <iostream>
   2 #include <fstream>
   3 #include <vector>
   4 #include <queue>
   5 #include <math.h>
   6 #include <string>
   7
   8 #define WINDOW 25
   9
  10 using namespace std;
  11
  12 vector<int> profile(2000,0);
  13
  14 void split (const string& text, const string& separators, vector<string>& words);
  15
  16 class Hit {
  17   public:
  18     string chr;
  19     unsigned int pos;
  20     bool strand;
  21
  22     Hit(string chr, unsigned int pos, bool strand) {
  23       this->chr = chr; this->pos = pos; this->strand = strand;
  24     }
  25 };
  26
  27 typedef vector<Hit> Hits;
  28
  29 bool compare_hits(const Hit &a, const Hit &b)
  30 {
  31   if(a.chr == b.chr) {
  32     return a.pos < b.pos;
  33   } else {
  34     return a.chr < b.chr;
  35   }
  36 }
  37
  38 int main(int argc, char** argv) {
  39   if(argc < 2) { cerr << "Usage: " << argv[0] << " [aligned reads filename] [file with list of feature starts]\n"; exit(0); }
  40   char filename[128]; strcpy(filename,argv[1]);
  41   char feature_filename[128]; strcpy(feature_filename,argv[2]);
  42   Hits data;
  43
  44   Hits features;
  45
  46   string delim("\t");
  47   string location_delim(":");
  48
  49   ifstream feat(feature_filename);
  50   size_t N = 0;
  51   while(feat.peek() != EOF) {
  52     char line[1024];
  53     feat.getline(line,1024,'\n');
  54     N++;
  55     string line_str(line);
  56     vector<string> fields;
  57     split(line_str, delim, fields);
  58     if(fields.size() != 3) { cerr << "Error: wrong number of fields in feature list (line " << N << " has " << fields.size() << " fields)\n"; }
  59
  60     string chr = fields[0];
  61
  62     int pos = atoi(fields[1].c_str());
  63     bool strand = (bool)(atoi(fields[2].c_str()));
  64
  65     Hit feat(chr,pos,strand);
  66     features.push_back(feat);
  67   }
  68
  69   cerr << "Found " << features.size() << " features\n";
  70
  71   //sort the features so we can run through it once
  72   std::sort(features.begin(),features.end(),compare_hits);
  73
  74   delim = " \n";
  75   location_delim = ":";
  76   char strand_str[2]; strand_str[1] = '\0';
  77   ifstream seqs(filename);
  78   unsigned int linenum = 0;
  79   while(seqs.peek() != EOF) {
  80     cerr << ++linenum << endl;
  81     char line[2048];
  82     seqs.getline(line,2048,'\n');
  83     string line_str(line);
  84     vector<string> fields;
  85     split(line_str, delim, fields);
  86     if(fields.size() == 3) { continue; }
  87
  88     vector<string> location; split(fields[3], location_delim, location);
  89     string chr = location[0];
  90
  91     int pos = atoi(location[1].c_str());
  92     bool strand = ((fields[4].c_str())[0] == 'F')?0:1;
  93
  94     Hit hit(chr,pos,strand);
  95     data.push_back(hit);
  96   }
  97   seqs.close();
  98
  99   cerr << "Found " << data.size() << " reads\n";
 100
 101   //sort the data so we can run through it once
 102   std::sort(data.begin(),data.end(),compare_hits);
 103
 104   string chrom = "";
 105   unsigned int feat_idx = 0;
 106   int num_hits = 0;
 107   for(Hits::iterator i = data.begin(); i != data.end(); ++i) {
 108     num_hits++;
 109     if(chrom == "" || i->chr != chrom) {
 110       chrom = i->chr;
 111       feat_idx = 0;
 112       while(feat_idx < features.size() && features[feat_idx].chr != chrom) { feat_idx++; }
 113       if(feat_idx == features.size()) { break; }
 114       cerr << chrom.c_str() << " feat_idx: " << feat_idx << endl;
 115     }
 116
 117     int dist_to_feature = i->pos - features[feat_idx].pos;
 118     //if we have passed the last feature, fast forward to the next
 119     while( feat_idx < features.size() && ((features[feat_idx].strand == 0 && dist_to_feature > 1000) || (features[feat_idx].strand == 1 && dist_to_feature > 1000))) {
 120       if(features[feat_idx].chr != i->chr ) { goto end_loop; }
 121       feat_idx++;
 122       dist_to_feature = i->pos - features[feat_idx].pos;
 123     }
 124
 125     if(features[feat_idx].strand == 1) { dist_to_feature *= -1; }
 126
 127     if(dist_to_feature > -1000 && dist_to_feature < 1000) {
 128       profile[dist_to_feature + 1000]++;
 129     }
 130   end_loop:
 131   ;
 132   }
 133
 134   for(unsigned int i = 0; i < profile.size(); i++) {
 135     cout << (int)i - (int)1000 << "\t" << (double)profile[i] / (double)data.size() << endl;;
 136   }
 137 }
 138
 139 void split (const string& text, const string& separators, vector<string>& words) {
 140
 141     size_t n     = text.length ();
 142     size_t start = text.find_first_not_of (separators);
 143
 144     while (start < n) {
 145         size_t stop = text.find_first_of (separators, start);
 146         if (stop > n) stop = n;
 147         words.push_back (text.substr (start, stop-start));
 148         start = text.find_first_not_of (separators, stop+1);
 149     }
 150 }