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   while(seqs.peek() != EOF) {
  79     char line[2048];
  80     seqs.getline(line,2048,'\n');
  81     string line_str(line);
  82     vector<string> fields;
  83     split(line_str, delim, fields);
  84     if(fields.size() == 3) { continue; }
  85
  86     vector<string> location; split(fields[3], location_delim, location);
  87     string chr = location[0];
  88
  89     int pos = atoi(location[1].c_str());
  90     bool strand = ((fields[4].c_str())[0] == 'F')?0:1;
  91
  92     Hit hit(chr,pos,strand);
  93     data.push_back(hit);
  94   }
  95   seqs.close();
  96
  97   cerr << "Found " << data.size() << " reads\n";
  98
  99   //sort the data so we can run through it once
 100   std::sort(data.begin(),data.end(),compare_hits);
 101
 102   string chrom = "";
 103   unsigned int feat_idx = 0;
 104   int num_hits = 0;
 105   for(Hits::iterator i = data.begin(); i != data.end(); ++i) {
 106     num_hits++;
 107     if(chrom == "" || i->chr != chrom) {
 108       chrom = i->chr;
 109       feat_idx = 0;
 110       while(feat_idx < features.size() && features[feat_idx].chr != chrom) { feat_idx++; }
 111       if(feat_idx == features.size()) { break; }
 112       cerr << chrom.c_str() << " feat_idx: " << feat_idx << endl;
 113     }
 114
 115     int dist_to_feature = i->pos - features[feat_idx].pos;
 116     //if we have passed the last feature, fast forward to the next
 117     while( feat_idx < features.size() && ((features[feat_idx].strand == 0 && dist_to_feature > 1000) || (features[feat_idx].strand == 1 && dist_to_feature > 1000))) {
 118       if(features[feat_idx].chr != i->chr ) { goto end_loop; }
 119       feat_idx++;
 120       dist_to_feature = i->pos - features[feat_idx].pos;
 121     }
 122
 123     if(features[feat_idx].strand == 1) { dist_to_feature *= -1; }
 124
 125     if(dist_to_feature > -1000 && dist_to_feature < 1000) {
 126       profile[dist_to_feature + 1000]++;
 127     }
 128   end_loop:
 129   ;
 130   }
 131
 132   for(unsigned int i = 0; i < profile.size(); i++) {
 133     cout << (int)i - (int)1000 << "\t" << (double)profile[i] / (double)data.size() << endl;;
 134   }
 135 }
 136
 137 void split (const string& text, const string& separators, vector<string>& words) {
 138
 139     size_t n     = text.length ();
 140     size_t start = text.find_first_not_of (separators);
 141
 142     while (start < n) {
 143         size_t stop = text.find_first_of (separators, start);
 144         if (stop > n) stop = n;
 145         words.push_back (text.substr (start, stop-start));
 146         start = text.find_first_not_of (separators, stop+1);
 147     }
 148 }