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