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