htswanalysis/src/complexity_count.cpp

   1 #include <iostream>
   2 #include <fstream>
   3 #include <vector>
   4 #include <map>
   5 #include <queue>
   6 #include <math.h>
   7 #include <string>
   8
   9 #include <gsl/gsl_statistics.h>
  10
  11 #define WINDOW 25
  12
  13 using namespace std;
  14
  15
  16 void int_to_chromstr(int chrom, char* chr);
  17 int chromstr_to_int(string chrom);
  18
  19 void split (const string& text, const string& separators, vector<string>& words);
  20
  21 class Hit {
  22   public:
  23     int chr;
  24     int start;
  25     int end;
  26     bool strand;
  27
  28     string name;
  29
  30     unsigned int count[2];
  31
  32     Hit(string name, int chr, int start, int end, bool strand) { this->name = name; this->chr = chr; this->start = start; this->end = end; this->strand = strand; count[0] = 0; count[1] = 0; }
  33 };
  34
  35 ostream &operator<<( ostream &out, const Hit &h ) {
  36   out << h.name << "\tchr" << h.chr << "\t" << h.start << "\t" << h.end << "\t" << h.strand << "\t" << h.count[0] << "\t" << h.count[1];
  37   return out;
  38 }
  39
  40 typedef vector<Hit> Hits;
  41
  42 bool compare_hits(const Hit &a, const Hit &b) {
  43   if(a.chr == b.chr) {
  44     return a.end < b.start;
  45   } else {
  46     return a.chr < b.chr;
  47   }
  48 }
  49
  50 void read_align_file(char* filename, Hits& features) {
  51   string delim = " \n";
  52   string location_delim = ":";
  53   char strand_str[2]; strand_str[1] = '\0';
  54   ifstream seqs(filename);
  55   string name("");
  56   while(seqs.peek() != EOF) {
  57     char line[2048];
  58     seqs.getline(line,2048,'\n');
  59     string line_str(line);
  60     vector<string> fields;
  61     split(line_str, delim, fields);
  62     if(fields.size() <= 3) { continue; }
  63     vector<string> location; split(fields[3], location_delim, location);
  64     int chr = chromstr_to_int(location[0]);
  65     if(chr == -1) { continue; }
  66     int pos = atoi(location[1].c_str());
  67     bool strand = ((fields[4].c_str())[0] == 'F')?0:1;
  68
  69     if(strand == 0) {
  70       Hit hit(name,chr,pos,pos,strand);
  71       features.push_back(hit);
  72     } else {
  73       Hit hit(name,chr,pos+25,pos+25,strand);
  74       features.push_back(hit);
  75     }
  76   }
  77   seqs.close();
  78
  79   //sort the data so we can run through it once
  80   std::stable_sort(features.begin(),features.end(),compare_hits);
  81 }
  82
  83 void print_features(Hits& hits, string filename, unsigned int size_1, unsigned int size_2) {
  84   ofstream out(filename.c_str());
  85   for(Hits::iterator i = hits.begin(); i != hits.end(); ++i) {
  86     Hit h = (*i);
  87     out << h.name << "\tchr" << h.chr << "\t" << h.start << "\t" << h.end << "\t" << h.strand << "\t" << h.count[0]/(double)(size_1/1000000.0) << "\t" << h.count[1]/(double)(size_2/1000000.0) << "\t" << (double)log(((double)h.count[0]/(double)size_1)/(double)(h.count[1]/(double)size_2))/log(2) << endl;
  88   }
  89   out.close();
  90 }
  91
  92 unsigned int count_complexity(Hits& data) {
  93   int chrom = -1;
  94
  95   unsigned int n_island = 0;
  96   unsigned int island_size = 0;
  97   for(Hits::iterator i = data.begin()+1; i != data.end()-1 && i != data.end(); ++i) {
  98     if(chrom == -1 || i->chr != chrom) {
  99       chrom = i->chr;
 100     } else if( (i+1)->chr != chrom ) { continue; }
 101     int prev = (i-1)->start;
 102     int count = 1;
 103     while( (i+1) != data.end() && (i+1)->chr == i->chr && (i+1)->start == i->start) { i++; count++; }
 104     int next = (i+1)->start;
 105
 106     if(i->chr != (i+1)->chr || (i-1)->chr != i->chr) { continue; }
 107
 108     if(prev < i->start - 25 && next > i->start + 25 && count >= 5) { n_island++; island_size += count; }
 109   }
 110   cerr << n_island << " Islands\n";
 111   return island_size;
 112 }
 113
 114 int main(int argc, char** argv) {
 115   if(argc < 3) { cerr << "Usage: " << argv[0] << " label reads\n"; exit(1); }
 116
 117   char label1[128]; strcpy(label1,argv[1]);
 118   Hits align1; read_align_file(argv[2], align1);
 119
 120   int complexity;
 121   if(align1.size() == 0) complexity = 0; else complexity = count_complexity(align1);
 122
 123   cout << label1 << "\t" << complexity << "\t" << align1.size()/1000000 << "M\t" << (double)complexity/(double)align1.size() << endl;
 124 }
 125
 126 int chromstr_to_int(string chromstr) {
 127   const char* chrom = chromstr.c_str();
 128   if(chrom[0] != 'c' || chrom[1] != 'h' || chrom[2] != 'r') { return -1; }
 129   if(chrom[3] == 'X') { return 23; } else
 130   if(chrom[3] == 'Y') { return 24; } else
 131   if(chrom[3] == 'M') { return 25; } else
 132   if(chrom[3] >= '0' && chrom[3] <= '9') { return atoi(chrom+3); } else
 133   { return -1; }
 134 }
 135
 136 void int_to_chromstr(int chrom, char* chr) {
 137   strcpy(chr,"chr");
 138   switch(chrom) {
 139     case 23: chr[3]='X'; break;
 140     case 24: chr[3]='Y'; break;
 141     case 25: chr[3]='M'; break;
 142     default: if(chrom < 10) { chr[3]=(char)(48+chrom); chr[4] = '\0'; }
 143              else { chr[3]=(char)(48+(int)(chrom/10));
 144                     chr[4]=(char)(48+(int)(chrom - (int)(chr[3]-48)));
 145                     chr[5]='\0';
 146                   }
 147   }
 148 }
 149
 150
 151 void split (const string& text, const string& separators, vector<string>& words) {
 152
 153     size_t n     = text.length ();
 154     size_t start = text.find_first_not_of (separators);
 155
 156     while (start < n) {
 157         size_t stop = text.find_first_of (separators, start);
 158         if (stop > n) stop = n;
 159         words.push_back (text.substr (start, stop-start));
 160         start = text.find_first_not_of (separators, stop+1);
 161     }
 162 }