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 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
  69   assert(fin.is_open());
  70
  71   names.clear(); names.push_back("");
  72   seqs.clear(); seqs.push_back("");
  73
  74   getline(fin, line);
  75   while ((fin) && (line[0] == '>')) {
  76     tag = line.substr(1);
  77     rawseq = "";
  78     while((getline(fin, line)) && (line[0] != '>')) {
  79       rawseq += line;
  80     }
  81     if (rawseq.size() <= 0) {
  82       printf("Warning: Fasta entry %s has an empty sequence! It is omitted!\n", tag.c_str());
  83       continue;
  84     }
  85     names.push_back(tag);
  86     seqs.push_back(convert(rawseq));
  87   }
  88
  89   fin.close();
  90
  91   M = names.size() - 1;
  92
  93   printf("The reference is loaded.\n");
  94 }
  95
  96 int main(int argc, char* argv[]) {
  97   if (argc != 4) {
  98     printf("Usage: rsem-for-ebseq-calculate-clustering-info k input_reference_fasta_file output_file\n");
  99     exit(-1);
 100   }
 101
 102   k = atoi(argv[1]);
 103   loadRef(argv[2]);
 104
 105   cands.clear();
 106   effL.assign(M + 1, 0);
 107   for (INTEGER i = 1; i <= M; i++) {
 108     effL[i] = seqs[i].length() - k + 1;
 109     if (effL[i] <= 0) effL[i] = 0; // effL should be non-negative
 110     for (INTEGER j = 0; j < effL[i]; j++)
 111       cands.push_back(ReadType(i, j));
 112   }
 113   printf("All possbile %d mers are generated.\n", k);
 114
 115   sort(cands.begin(), cands.end());
 116   printf("All %d mers are sorted.\n", k);
 117
 118   size_t p = 0;
 119   clusteringInfo.assign(M + 1, 0.0);
 120
 121   for (size_t i = 1; i <= cands.size(); i++)
 122     if (i == cands.size() || !cands[p].seq_equal(cands[i])) {
 123       size_t denominator = i - p;
 124       size_t q = p;
 125       for (size_t j = p + 1; j <= i; j++)
 126         if (j == i || cands[q].tid != cands[j].tid) {
 127           size_t numerator = j - q;
 128           //double prob = numerator * 1.0 / denominator;
 129           //clusteringInfo[cands[q].tid] += (double)numerator * prob * (1.0 - prob);
 130           if (numerator < denominator) clusteringInfo[cands[q].tid] += numerator;
 131           q = j;
 132         }
 133       p = i;
 134     }
 135
 136   for (INTEGER i = 1; i <= M; i++)
 137     if (effL[i] == 0) clusteringInfo[i] = -1.0;
 138     else clusteringInfo[i] /= effL[i];
 139
 140   printf("Clustering information is calculated.\n");
 141
 142
 143   ofstream fout(argv[3]);
 144   for (INTEGER i = 1; i <= M; i++) fout<<names[i]<<"\t"<<setprecision(6)<<clusteringInfo[i]<<endl;
 145   fout.close();
 146
 147   return 0;
 148 }