calcCI.cpp

   1 #include<ctime>
   2 #include<cstdio>
   3 #include<cstring>
   4 #include<cstdlib>
   5 #include<cassert>
   6 #include<fstream>
   7 #include<algorithm>
   8
   9 #include "boost/random.hpp"
  10
  11 #include "utils.h"
  12
  13 #include "Model.h"
  14 #include "SingleModel.h"
  15 #include "SingleQModel.h"
  16 #include "PairedEndModel.h"
  17 #include "PairedEndQModel.h"
  18
  19 #include "Refs.h"
  20 #include "GroupInfo.h"
  21
  22 using namespace std;
  23
  24 typedef unsigned long bufsize_type;
  25 typedef boost::mt19937 engine_type;
  26 typedef boost::gamma_distribution<> distribution_type;
  27 typedef boost::variate_generator<engine_type&, distribution_type> generator_type;
  28
  29 const int FLOATSIZE = sizeof(float);
  30
  31 struct CIType {
  32         float lb, ub; // the interval is [lb, ub]
  33
  34         CIType() { lb = ub = 0.0; }
  35 };
  36
  37 bool quiet;
  38 int model_type;
  39
  40 double confidence;
  41
  42 int nC, cvlen, nSpC, nSamples; // nSpC : number of sample theta vectors per count vector
  43 int fr, to; // each flush, sample fr .. to - 1
  44
  45 int nMB;
  46 bufsize_type size;
  47 float *buffer;
  48 char cvsF[STRLEN], tmpF[STRLEN], command[STRLEN];
  49 ofstream ftmpOut;
  50
  51 int *cvec;
  52 double *theta;
  53 CIType *iso_tau, *gene_tau;
  54
  55 engine_type engine(time(NULL));
  56 distribution_type **gammas;
  57 generator_type **rgs;
  58
  59 int M, m;
  60 Refs refs;
  61 GroupInfo gi;
  62 char imdName[STRLEN], statName[STRLEN];
  63 char modelF[STRLEN], groupF[STRLEN], refF[STRLEN];
  64
  65 vector<double> eel; //expected effective lengths
  66 double *tau_denoms; // denominators for tau values
  67
  68 template<class ModelType>
  69 void calcExpectedEffectiveLengths(ModelType& model) {
  70   int lb, ub, span;
  71   double *pdf = NULL, *cdf = NULL, *clen = NULL; // clen[i] = sigma_{j=1}^{i}pdf[i]*(lb+i)
  72
  73   model.getGLD().copyTo(pdf, cdf, lb, ub, span);
  74   clen = new double[span + 1];
  75   clen[0] = 0.0;
  76   for (int i = 1; i <= span; i++) {
  77     clen[i] = clen[i - 1] + pdf[i] * (lb + i);
  78   }
  79
  80   eel.clear();
  81   eel.resize(M + 1, 0.0);
  82   for (int i = 1; i <= M; i++) {
  83     int totLen = refs.getRef(i).getTotLen();
  84     int fullLen = refs.getRef(i).getFullLen();
  85     int pos1 = max(min(totLen - fullLen + 1, ub) - lb, 0);
  86     int pos2 = max(min(totLen, ub) - lb, 0);
  87
  88     if (pos2 == 0) { eel[i] = 0.0; continue; }
  89
  90     eel[i] = fullLen * cdf[pos1] + ((cdf[pos2] - cdf[pos1]) * (totLen + 1) - (clen[pos2] - clen[pos1]));
  91     assert(eel[i] >= 0);
  92     if (eel[i] < MINEEL) { eel[i] = 0.0; }
  93   }
  94
  95   delete[] pdf;
  96   delete[] cdf;
  97   delete[] clen;
  98 }
  99
 100 void flushToTempFile() {
 101         int gap1 = fr * FLOATSIZE;
 102         int gap2 = (nSamples - to) * FLOATSIZE;
 103         float *p = NULL;
 104
 105         ftmpOut.seekp(0, ios_base::beg);
 106         for (int i = 0; i < cvlen; i++) {
 107                 p = buffer + i;
 108                 ftmpOut.seekp(gap1, ios_base::cur);
 109                 for (int j = fr; j < to; j++) {
 110                         ftmpOut.write((char*)p, FLOATSIZE);
 111                         p += cvlen;
 112                 }
 113                 ftmpOut.seekp(gap2, ios_base::cur);
 114         }
 115 }
 116
 117 template<class ModelType>
 118 void sampling() {
 119         float *p, *ub;
 120         ifstream fin(cvsF);
 121         ModelType model;
 122
 123         model.read(modelF);
 124         calcExpectedEffectiveLengths<ModelType>(model);
 125
 126         ftmpOut.open(tmpF, ios::binary);
 127
 128         fin>>nC>>cvlen;
 129         assert(cvlen = M + 1);
 130
 131         nSamples = nC * nSpC;
 132
 133         fr = to = 0;
 134
 135         size = bufsize_type(nMB) * 1024 * 1024 / FLOATSIZE / cvlen;
 136         if (size > (bufsize_type)nSamples) size = nSamples;
 137         size *= cvlen;
 138         buffer = new float[size];
 139
 140         ub = buffer + size;
 141         p = buffer;
 142
 143         cvec = new int[cvlen];
 144         theta = new double[cvlen];
 145         gammas = new distribution_type*[cvlen];
 146         rgs = new generator_type*[cvlen];
 147
 148         tau_denoms = new double[nSamples];
 149         memset(tau_denoms, 0, sizeof(double) * nSamples);
 150
 151         double *mw = model.getMW();
 152         for (int i = 0; i < nC; i++) {
 153                 for (int j = 0; j < cvlen; j++) {
 154                         fin>>cvec[j];
 155                 }
 156
 157                 for (int j = 0; j < cvlen; j++) {
 158                         gammas[j] = new distribution_type(cvec[j]); // need change back before publishing
 159                         rgs[j] = new generator_type(engine, *gammas[j]);
 160                 }
 161
 162                 for (int j = 0; j < nSpC; j++) {
 163                         double sum = 0.0;
 164                         for (int k = 0; k < cvlen; k++) {
 165                                 theta[k] = (k == 0 || eel[k] > EPSILON ? (*rgs[k])() : 0.0);
 166                                 sum += theta[k];
 167                         }
 168                         assert(sum > 0.0);
 169                         for (int k = 0; k < cvlen; k++) theta[k] /= sum;
 170
 171                         sum = 0.0;
 172                         for (int k = 0; k < cvlen; k++) {
 173                           theta[k] = (mw[k] < EPSILON ? 0.0 : theta[k] / mw[k]);
 174                           sum += theta[k];
 175                         }
 176                         assert(sum >= EPSILON);
 177                         for (int k = 0; k < cvlen; k++) theta[k] /= sum;
 178
 179                         *p = (float)theta[0]; ++p;
 180                         assert(1.0 - theta[0] > 0.0);
 181                         for (int k = 1; k < cvlen; k++) {
 182                                 if (eel[k] > EPSILON) {
 183                                         theta[k] /= (1.0 - theta[0]);
 184                                         tau_denoms[to] += theta[k] / eel[k];
 185                                 }
 186                                 else {
 187                                         if (theta[k] != 0.0) { fprintf(stderr, "K=%d Theta_K=%lf\n", k, theta[k]); exit(-1); }
 188                                 }
 189
 190                                 *p = (float)theta[k];
 191                                 ++p;
 192                         }
 193                         ++to;
 194                         if (p == ub) {
 195                                 flushToTempFile();
 196                                 p = buffer;
 197                                 fr = to;
 198                                 if (verbose) { printf("%d vectors are sampled!\n", to); }
 199                         }
 200                 }
 201
 202                 for (int j = 0; j < cvlen; j++) {
 203                         delete gammas[j];
 204                         delete rgs[j];
 205                 }
 206         }
 207
 208         if (fr != to) { flushToTempFile(); }
 209
 210         fin.close();
 211         ftmpOut.close();
 212
 213         delete[] buffer;
 214
 215         delete[] cvec;
 216         delete[] theta;
 217         delete[] gammas;
 218         delete[] rgs;
 219
 220         if (verbose) { printf("Sampling is finished!\n"); }
 221 }
 222
 223 void calcCI(int nSamples, float *samples, float &lb, float &ub) {
 224         int p, q; // p pointer for lb, q pointer for ub;
 225         int newp, newq;
 226         int threshold = nSamples - (int(confidence * nSamples - 1e-8) + 1);
 227         int nOutside = 0;
 228
 229         sort(samples, samples + nSamples);
 230
 231         p = 0; q = nSamples - 1;
 232         newq = nSamples - 1;
 233         do {
 234                 q = newq;
 235                 while (newq > 0 && samples[newq - 1] == samples[newq]) newq--;
 236                 newq--;
 237         } while (newq >= 0 && nSamples - (newq + 1) <= threshold);
 238
 239         nOutside = nSamples - (q + 1);
 240
 241         lb = -1e30; ub = 1e30;
 242         do {
 243                 if (samples[q] - samples[p] < ub - lb) {
 244                         lb = samples[p];
 245                         ub = samples[q];
 246                 }
 247
 248                 newp = p;
 249                 while (newp < nSamples - 1 && samples[newp] == samples[newp + 1]) newp++;
 250                 newp++;
 251                 if (newp <= threshold) {
 252                         nOutside += newp - p;
 253                         p = newp;
 254                         while (nOutside > threshold && q < nSamples - 1) {
 255                                 newq = q + 1;
 256                                 while (newq < nSamples - 1 && samples[newq] == samples[newq + 1]) newq++;
 257                                 nOutside -= newq - q;
 258                                 q = newq;
 259                         }
 260                         assert(nOutside <= threshold);
 261                 }
 262                 else p = newp;
 263         } while (p <= threshold);
 264 }
 265
 266 void generateResults(char* imdName) {
 267         float *itsamples, *gtsamples;
 268         ifstream fin;
 269         FILE *fo;
 270         char outF[STRLEN];
 271
 272         iso_tau = new CIType[M + 1];
 273         gene_tau = new CIType[m];
 274
 275         itsamples = new float[nSamples];
 276         gtsamples = new float[nSamples];
 277
 278         fin.open(tmpF, ios::binary);
 279
 280         //read sampled theta0 values
 281         for (int k = 0; k < nSamples; k++) fin.read((char*)(&itsamples[k]), FLOATSIZE);
 282
 283         for (int i = 0; i < m; i++) {
 284                 int b = gi.spAt(i), e = gi.spAt(i + 1);
 285                 memset(gtsamples, 0, FLOATSIZE * nSamples);
 286                 for (int j = b; j < e; j++) {
 287                         for (int k = 0; k < nSamples; k++) {
 288                                 fin.read((char*)(&itsamples[k]), FLOATSIZE);
 289                                 if (eel[j] > EPSILON && tau_denoms[k] > EPSILON) { itsamples[k] = itsamples[k] / eel[j] / tau_denoms[k]; }
 290                                 else {
 291                                         if (itsamples[k] != 0.0) { fprintf(stderr, "Not equal to 0! K=%d, Sampled Theta Value=%lf\n", k, itsamples[k]); exit(-1); }
 292                                         itsamples[k] = 0.0;
 293                                 }
 294                                 gtsamples[k] += itsamples[k];
 295                         }
 296                         calcCI(nSamples, itsamples, iso_tau[j].lb, iso_tau[j].ub);
 297                 }
 298                 calcCI(nSamples, gtsamples, gene_tau[i].lb, gene_tau[i].ub);
 299
 300                 if (verbose && (i + 1) % 1000 == 0) { printf("%d genes are done!\n", i + 1); }
 301         }
 302
 303         fin.close();
 304
 305         //isoform level results
 306         sprintf(outF, "%s.iso_res", imdName);
 307         fo = fopen(outF, "a");
 308         for (int i = 1; i <= M; i++)
 309                 fprintf(fo, "%.6g%c", iso_tau[i].lb, (i < M ? '\t' : '\n'));
 310         for (int i = 1; i <= M; i++)
 311                 fprintf(fo, "%.6g%c", iso_tau[i].ub, (i < M ? '\t' : '\n'));
 312         fclose(fo);
 313
 314         //gene level results
 315         sprintf(outF, "%s.gene_res", imdName);
 316         fo = fopen(outF, "a");
 317         for (int i = 0; i < m; i++)
 318                 fprintf(fo, "%.6g%c", gene_tau[i].lb, (i < m - 1 ? '\t' : '\n'));
 319         for (int i = 0; i < m; i++)
 320                 fprintf(fo, "%.6g%c", gene_tau[i].ub, (i < m - 1 ? '\t' : '\n'));
 321         fclose(fo);
 322
 323         delete[] itsamples;
 324         delete[] gtsamples;
 325
 326         delete[] iso_tau;
 327         delete[] gene_tau;
 328
 329         if (verbose) { printf("All credibility intervals are calculated!\n"); }
 330
 331         sprintf(outF, "%s.tau_denoms", imdName);
 332         fo = fopen(outF, "w");
 333         fprintf(fo, "%d\n", nSamples);
 334         for (int i = 0; i < nSamples; i++) fprintf(fo, "%.15g ", tau_denoms[i]);
 335         fprintf(fo, "\n");
 336         fclose(fo);
 337 }
 338
 339 int main(int argc, char* argv[]) {
 340         if (argc < 7) {
 341                 printf("Usage: rsem-calculate-credibility-intervals reference_name sample_name sampleToken confidence nSpC nMB[-q]\n");
 342                 exit(-1);
 343         }
 344
 345         confidence = atof(argv[4]);
 346         nSpC = atoi(argv[5]);
 347         nMB = atoi(argv[6]);
 348
 349         quiet = false;
 350         if (argc > 7 && !strcmp(argv[7], "-q")) {
 351                 quiet = true;
 352         }
 353         verbose = !quiet;
 354
 355         sprintf(imdName, "%s.temp/%s", argv[2], argv[3]);
 356         sprintf(statName, "%s.stat/%s", argv[2], argv[3]);
 357
 358         sprintf(modelF, "%s.model", statName);
 359         FILE *fi = fopen(modelF, "r");
 360         if (fi == NULL) { fprintf(stderr, "Cannot open %s!\n", modelF); exit(-1); }
 361         fscanf(fi, "%d", &model_type);
 362         fclose(fi);
 363
 364         sprintf(refF, "%s.seq", argv[1]);
 365         refs.loadRefs(refF, 1);
 366         M = refs.getM();
 367         sprintf(groupF, "%s.grp", argv[1]);
 368         gi.load(groupF);
 369         m = gi.getm();
 370
 371         sprintf(tmpF, "%s.tmp", imdName);
 372         sprintf(cvsF, "%s.countvectors", imdName);
 373
 374         switch(model_type) {
 375         case 0 : sampling<SingleModel>(); break;
 376         case 1 : sampling<SingleQModel>(); break;
 377         case 2 : sampling<PairedEndModel>(); break;
 378         case 3 : sampling<PairedEndQModel>(); break;
 379         }
 380
 381         generateResults(imdName);
 382
 383         delete[] tau_denoms;
 384
 385         sprintf(command, "rm -f %s", tmpF);
 386         int status = system(command);
 387         if (status != 0) {
 388                 fprintf(stderr, "Cannot delete %s!\n", tmpF);
 389                 exit(-1);
 390         }
 391
 392         return 0;
 393 }