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