htswanalysis/src/GetReadsInSnps/getReadsInSnps.cpp

   1 #include <sys/types.h>
   2 #include <dirent.h>
   3 #include <errno.h>
   4 #include <iostream>
   5 #include <fstream>
   6 #include <vector>
   7 #include <map>
   8 #include <queue>
   9 #include <math.h>
  10 #include <string>
  11
  12 #include <gsl/gsl_statistics.h>
  13
  14 #define WINDOW 25
  15
  16 using namespace std;
  17
  18 void split (const string& text, const string& separators, vector<string>& words);
  19 char *strrevcomp(const string& input);
  20
  21 class Read {
  22   public:
  23     string chr;
  24     unsigned int pos;
  25     string seq;
  26
  27     Read(string chr, unsigned int pos, string seq) { this->chr = chr; this->pos = pos; this->seq = seq; }
  28     Read(const Read& r) { this->chr = r.chr; this->pos = r.pos; this->seq = r.seq; }
  29     Read& operator=(const Read& r) {this->chr = r.chr; this->pos = r.pos; this->seq = r.seq; return *this;}
  30 };
  31
  32 class SNP {
  33   public:
  34
  35     string name;
  36     string chr;
  37     unsigned int pos;
  38     char reference_base;
  39     unsigned int A;
  40     unsigned int C;
  41     unsigned int G;
  42     unsigned int T;
  43     unsigned int N;
  44     unsigned int total;
  45
  46     SNP(string name, string chr, unsigned int pos, char reference_base) {
  47       this->name = name;
  48       this->chr = chr;
  49       this->pos = pos;
  50       this->A = 0; this->C = 0; this->G = 0; this->T = 0; this->total = 0;
  51       this->reference_base = reference_base;
  52     }
  53
  54     SNP(const SNP& h) {
  55       this->name = h.name;
  56       this->chr = h.chr;
  57       this->pos = h.pos;
  58       this->A = h.A; this->C = h.C; this->G = h.G; this->T = h.T; this->total = h.total;
  59       this->reference_base = h.reference_base;
  60     }
  61
  62     SNP& operator=(const SNP& h) {
  63       this->name = h.name;
  64       this->chr = h.chr;
  65       this->pos = h.pos;
  66       this->A = h.A; this->C = h.C; this->G = h.G; this->T = h.T; this->total = h.total;
  67       this->reference_base = h.reference_base;
  68       return *this;
  69     }
  70
  71     void add_read(char nuc) {
  72       switch(nuc) {
  73         case 'a':
  74         case 'A':
  75           A++; break;
  76         case 'c':
  77         case 'C':
  78           C++; break;
  79         case 'g':
  80         case 'G':
  81           G++; break;
  82         case 't':
  83         case 'T':
  84           T++; break;
  85         default:
  86           N++; break;
  87       }
  88       total++;
  89     }
  90 };
  91
  92 bool SNP_lt_Read(const SNP& feat, const Read& read) {
  93   //First compare based on chromosomes
  94   if(read.chr == feat.chr) { return feat.pos < read.pos; }
  95   else { return feat.chr < read.chr; }
  96 }
  97
  98 ostream &operator<<( ostream &out, const SNP &h ) {
  99   out << h.name.c_str() << "\t" << h.chr.c_str() << "\t" << h.pos << "\t" << h.reference_base << "\t" << h.A << "\t" << h.C << "\t" << h.G << "\t" << h.T << "\t" << h.total;
 100   return out;
 101 }
 102
 103 typedef vector<Read> Reads;
 104 typedef vector<SNP> SNPs;
 105
 106 bool compare_snps(const SNP &a, const SNP &b) { if(a.chr == b.chr) { return a.pos < b.pos; } else { return a.pos < b.pos; } }
 107 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; } }
 108
 109 void read_snps(const char* filename, SNPs& snps) {
 110   string delim("\t");
 111
 112   ifstream feat(filename);
 113   size_t N = 0;
 114   while(feat.peek() != EOF) {
 115     char line[1024];
 116     feat.getline(line,1024,'\n');
 117     N++;
 118     string line_str(line);
 119     vector<string> fields;
 120     split(line_str, delim, fields);
 121     if(fields.size() != 4) { cerr << "Error (" << filename << "): wrong number of fields in feature list (line " << N << " has " << fields.size() << " fields)\n"; }
 122
 123     string name = fields[0];
 124     string chr = fields[1];
 125     unsigned int pos = atoi(fields[2].c_str());
 126     char base = (fields[3])[0];
 127
 128     SNP snp(name,chr,pos,base);
 129     snps.push_back(snp);
 130   }
 131
 132   //sort the features so we can run through it once
 133   std::stable_sort(snps.begin(),snps.end(),compare_snps);
 134   feat.close();
 135
 136   cerr << "Found and sorted " << snps.size() << " snps." << endl;
 137 }
 138
 139 void read_align_file(char* filename, Reads& features) {
 140   string delim(" \n");
 141   string location_delim(":");
 142   char strand_str[2]; strand_str[1] = '\0';
 143   ifstream seqs(filename);
 144   string name("");
 145   while(seqs.peek() != EOF) {
 146     char line[2048];
 147     seqs.getline(line,2048,'\n');
 148
 149     string line_str(line);
 150     vector<string> fields;
 151     split(line_str, delim, fields);
 152     if(fields.size() != 7) { continue; }
 153
 154
 155     vector<string> location; split(fields[3], location_delim, location);
 156     string chr = location[0];
 157     if(chr == "NA") { continue; }
 158     int pos = atoi(location[1].c_str());
 159     bool strand = ((fields[4].c_str())[0] == 'F')?0:1;
 160
 161     string seq;
 162     if(strand == 0) {
 163       seq = fields[0];
 164     } else {
 165       seq = strrevcomp(fields[0]);
 166     }
 167     Read read(chr,pos,seq);
 168     features.push_back(read);
 169   }
 170   seqs.close();
 171
 172   //sort the data so we can run through it once
 173   std::sort(features.begin(),features.end(),compare_reads);
 174   cerr << "Found and sorted " << features.size() << " reads." << endl;
 175 }
 176
 177 void count_read_at_snps(SNPs& snps, Reads& data) {
 178   SNPs::iterator snp_it = snps.begin();
 179
 180   //assume, for now, that all reads have the same length
 181   unsigned int read_len = data[0].seq.length();
 182
 183   for(Reads::iterator i = data.begin(); i != data.end(); ++i) {
 184     //skip to first feature after read
 185     string start_chr = snp_it->chr;
 186     while(snp_it != snps.end() && SNP_lt_Read(*snp_it, *i)) {
 187       snp_it++;
 188     }
 189
 190     //stop if we have run out of features.
 191     if(snp_it == snps.end()) { break; }
 192
 193     if(i->pos + read_len > snp_it->pos && i->pos <= snp_it->pos) {
 194       cout << i->chr.c_str() << "\t" << i->pos << "\t" << i->seq.c_str() << endl;
 195       snp_it->add_read(i->seq[snp_it->pos - i->pos]);
 196     }
 197   }
 198 }
 199
 200 int main(int argc, char** argv) {
 201   if(argc != 3) { cerr << "Usage: " << argv[0] << " snp_file read_file\n"; exit(1); }
 202
 203   char snp_filename[1024]; strcpy(snp_filename,argv[1]);
 204   char read_filename[1024]; strcpy(read_filename,argv[2]);
 205
 206   SNPs snps; read_snps(snp_filename, snps);
 207   Reads reads; read_align_file(read_filename, reads);
 208
 209   count_read_at_snps(snps, reads);
 210
 211   for(SNPs::iterator i = snps.begin(); i != snps.end(); ++i) {
 212     cout << (*i) << endl;
 213   }
 214 }
 215
 216 void split (const string& text, const string& separators, vector<string>& words) {
 217
 218     size_t n     = text.length ();
 219     size_t start = text.find_first_not_of (separators);
 220
 221     while (start < n) {
 222         size_t stop = text.find_first_of (separators, start);
 223         if (stop > n) stop = n;
 224         words.push_back (text.substr (start, stop-start));
 225         start = text.find_first_not_of (separators, stop+1);
 226     }
 227 }
 228
 229 char *strrevcomp(const string& input)
 230 {
 231   char* str = new char[input.length()];
 232   strcpy(str,input.c_str());
 233
 234   char *p1, *p2;
 235
 236   if (! str || ! *str)
 237     return str;
 238
 239   for (p1 = str, p2 = str + strlen(str) - 1; p2 > p1; ++p1, --p2) {
 240     *p1 ^= *p2;
 241     *p2 ^= *p1;
 242     *p1 ^= *p2;
 243   }
 244
 245   for (p1 = str; p1 < str + strlen(str); ++p1) {
 246     if(*p1 == 'a' || *p1 == 'A') { *p1 = 'T'; }
 247     else if(*p1 == 'c' || *p1 == 'C') { *p1 = 'G'; }
 248     else if(*p1 == 'g' || *p1 == 'G') { *p1 = 'C'; }
 249     else if(*p1 == 't' || *p1 == 'T') { *p1 = 'A'; }
 250   }
 251
 252   return str;
 253 }
 254