EBSeq/calcClusteringInfo.cpp

   1 #include<cstdio>
   2 #include<cctype>
   3 #include<cstring>
   4 #include<cstdlib>
   5 #include<cassert>
   6 #include<fstream>
   7 #include<iomanip>
   8 #include<string>
   9 #include<vector>
  10 #include<algorithm>
  11 using namespace std;
  12
  13 const int STRLEN = 1005;
  14
  15 int M;
  16 int k; // k-mer size
  17 vector<string> names;
  18 vector<string> seqs;
  19 vector<int> effL;
  20
  21 // tid starts from 1
  22 struct ReadType {
  23   int tid, pos;
  24
  25   ReadType(int tid, int pos) {
  26     this->tid = tid;
  27     this->pos = pos;
  28   }
  29
  30   bool operator< (const ReadType& o) const {
  31     string& a = seqs[tid];
  32     string& b = seqs[o.tid];
  33     for (int i = 0; i < k; i++) {
  34       if (a[pos + i] != b[o.pos + i]) {
  35         return a[pos + i] < b[o.pos + i];
  36       }
  37     }
  38     return tid < o.tid;
  39   }
  40
  41   bool seq_equal(const ReadType& o) const {
  42     string& a = seqs[tid];
  43     string& b = seqs[o.tid];
  44     for (int i = 0; i < k; i++)
  45       if (a[pos + i] != b[o.pos + i]) return false;
  46     return true;
  47   }
  48 };
  49
  50 vector<ReadType> cands;
  51 vector<double> clusteringInfo;
  52
  53 string convert(const string& rawseq) {
  54   int size = (int)rawseq.size();
  55   string seq = rawseq;
  56   for (int i = 0; i < size; i++) {
  57     seq[i] = toupper(rawseq[i]);
  58     if (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G' && seq[i] != 'T') seq[i] = 'N';
  59   }
  60   return seq;
  61 }
  62
  63 void loadRef(char* inpF) {
  64   ifstream fin(inpF);
  65   string tag, line, rawseq;
  66
  67   assert(fin.is_open());
  68
  69   names.clear(); names.push_back("");
  70   seqs.clear(); seqs.push_back("");
  71
  72   getline(fin, line);
  73   while ((fin) && (line[0] == '>')) {
  74     tag = line.substr(1);
  75     rawseq = "";
  76     while((getline(fin, line)) && (line[0] != '>')) {
  77       rawseq += line;
  78     }
  79     if (rawseq.size() <= 0) {
  80       printf("Warning: Fasta entry %s has an empty sequence! It is omitted!\n", tag.c_str());
  81       continue;
  82     }
  83     names.push_back(tag);
  84     seqs.push_back(convert(rawseq));
  85   }
  86
  87   fin.close();
  88
  89   M = names.size() - 1;
  90
  91   printf("The reference is loaded.\n");
  92 }
  93
  94 int main(int argc, char* argv[]) {
  95   if (argc != 4) {
  96     printf("Usage: rsem-for-ebseq-calculate-clustering-info k input_reference_fasta_file output_file\n");
  97     exit(-1);
  98   }
  99
 100   k = atoi(argv[1]);
 101   loadRef(argv[2]);
 102
 103   cands.clear();
 104   effL.assign(M + 1, 0);
 105   for (int i = 1; i <= M; i++) {
 106     effL[i] = int(seqs[i].length()) - k + 1;
 107     if (effL[i] <= 0) effL[i] = 0; // effL should be non-negative
 108     for (int j = 0; j < effL[i]; j++)
 109       cands.push_back(ReadType(i, j));
 110   }
 111   printf("All possbile %d mers are generated.\n", k);
 112
 113   sort(cands.begin(), cands.end());
 114   printf("All %d mers are sorted.\n", k);
 115
 116   size_t p = 0;
 117   clusteringInfo.assign(M + 1, 0.0);
 118
 119   for (size_t i = 1; i <= cands.size(); i++)
 120     if (i == cands.size() || !cands[p].seq_equal(cands[i])) {
 121       size_t denominator = i - p;
 122       size_t q = p;
 123       for (size_t j = p + 1; j <= i; j++)
 124         if (j == i || cands[q].tid != cands[j].tid) {
 125           size_t numerator = j - q;
 126           //double prob = numerator * 1.0 / denominator;
 127           //clusteringInfo[cands[q].tid] += (double)numerator * prob * (1.0 - prob);
 128           if (numerator < denominator) clusteringInfo[cands[q].tid] += numerator;
 129           q = j;
 130         }
 131       p = i;
 132     }
 133
 134   for (int i = 1; i <= M; i++)
 135     if (effL[i] == 0) clusteringInfo[i] = -1.0;
 136     else clusteringInfo[i] /= effL[i];
 137
 138   printf("Clustering information is calculated.\n");
 139
 140
 141   ofstream fout(argv[3]);
 142   for (int i = 1; i <= M; i++) fout<<names[i]<<"\t"<<setprecision(6)<<clusteringInfo[i]<<endl;
 143   fout.close();
 144
 145   return 0;
 146 }