calcCI.cpp

   1 #include<ctime>
   2 #include<cstdio>
   3 #include<cstring>
   4 #include<cstdlib>
   5 #include<cassert>
   6 #include<fstream>
   7 #include<algorithm>
   8 #include<vector>
   9 #include<pthread.h>
  10
  11 #include "utils.h"
  12 #include "my_assert.h"
  13 #include "sampling.h"
  14
  15 #include "Model.h"
  16 #include "SingleModel.h"
  17 #include "SingleQModel.h"
  18 #include "PairedEndModel.h"
  19 #include "PairedEndQModel.h"
  20
  21 #include "Refs.h"
  22 #include "GroupInfo.h"
  23 #include "WriteResults.h"
  24
  25 #include "Buffer.h"
  26
  27 using namespace std;
  28
  29 struct Params {
  30         int no;
  31         FILE *fi;
  32         engine_type *engine;
  33         const double *mw;
  34 };
  35
  36 struct CIParams {
  37         int no;
  38         int start_gene_id, end_gene_id;
  39 };
  40
  41 struct CIType {
  42         float lb, ub; // the interval is [lb, ub]
  43
  44         CIType() { lb = ub = 0.0; }
  45 };
  46
  47 int model_type;
  48
  49 int nMB;
  50 double confidence;
  51 int nCV, nSpC, nSamples; // nCV: number of count vectors; nSpC: number of theta vectors sampled per count vector; nSamples: nCV * nSpC
  52 int nThreads;
  53
  54 float *l_bars;
  55
  56 char cvsF[STRLEN], tmpF[STRLEN], command[STRLEN];
  57
  58 CIType *tpm, *fpkm;
  59 CIType *iso_tpm = NULL, *iso_fpkm = NULL;
  60 CIType *gene_tpm, *gene_fpkm;
  61
  62 int M, m;
  63 Refs refs;
  64 GroupInfo gi;
  65 char refName[STRLEN], imdName[STRLEN], statName[STRLEN];
  66 char modelF[STRLEN], groupF[STRLEN], refF[STRLEN];
  67
  68 bool alleleS;
  69 int m_trans;
  70 GroupInfo ta;
  71
  72 vector<double> eel; //expected effective lengths
  73
  74 Buffer *buffer;
  75
  76 bool quiet;
  77
  78 Params *paramsArray;
  79 pthread_t *threads;
  80 pthread_attr_t attr;
  81 int rc;
  82
  83 bool hasSeed;
  84 seedType seed;
  85
  86 CIParams *ciParamsArray;
  87
  88 void* sample_theta_from_c(void* arg) {
  89         int *cvec;
  90         double *theta;
  91         float *tpm;
  92         gamma_dist **gammas;
  93         gamma_generator **rgs;
  94
  95         Params *params = (Params*)arg;
  96         FILE *fi = params->fi;
  97         const double *mw = params->mw;
  98
  99         cvec = new int[M + 1];
 100         theta = new double[M + 1];
 101         gammas = new gamma_dist*[M + 1];
 102         rgs = new gamma_generator*[M + 1];
 103         tpm = new float[M + 1];
 104         float l_bar; // the mean transcript length over the sample
 105
 106         int cnt = 0;
 107         while (fscanf(fi, "%d", &cvec[0]) == 1) {
 108                 for (int j = 1; j <= M; j++) assert(fscanf(fi, "%d", &cvec[j]) == 1);
 109
 110                 ++cnt;
 111
 112                 for (int j = 0; j <= M; j++) {
 113                         gammas[j] = new gamma_dist(cvec[j]);
 114                         rgs[j] = new gamma_generator(*(params->engine), *gammas[j]);
 115                 }
 116
 117                 for (int i = 0; i < nSpC; i++) {
 118                         double sum = 0.0;
 119                         for (int j = 0; j <= M; j++) {
 120                                 theta[j] = ((j == 0 || (eel[j] >= EPSILON && mw[j] >= EPSILON)) ? (*rgs[j])() / mw[j] : 0.0);
 121                                 sum += theta[j];
 122                         }
 123                         assert(sum >= EPSILON);
 124                         for (int j = 0; j <= M; j++) theta[j] /= sum;
 125
 126                         sum = 0.0;
 127                         tpm[0] = 0.0;
 128                         for (int j = 1; j <= M; j++)
 129                                 if (eel[j] >= EPSILON) {
 130                                         tpm[j] = theta[j] / eel[j];
 131                                         sum += tpm[j];
 132                                 }
 133                                 else assert(theta[j] < EPSILON);
 134                         assert(sum >= EPSILON);
 135                         l_bar = 0.0; // store mean effective length of the sample
 136                         for (int j = 1; j <= M; j++) { tpm[j] /= sum; l_bar += tpm[j] * eel[j]; tpm[j] *= 1e6; }
 137                         buffer->write(l_bar, tpm + 1); // ommit the first element in tpm
 138                 }
 139
 140                 for (int j = 0; j <= M; j++) {
 141                         delete gammas[j];
 142                         delete rgs[j];
 143                 }
 144
 145                 if (verbose && cnt % 100 == 0) { printf("Thread %d, %d count vectors are processed!\n", params->no, cnt); }
 146         }
 147
 148         delete[] cvec;
 149         delete[] theta;
 150         delete[] gammas;
 151         delete[] rgs;
 152         delete[] tpm;
 153
 154         return NULL;
 155 }
 156
 157 template<class ModelType>
 158 void sample_theta_vectors_from_count_vectors() {
 159         ModelType model;
 160         model.read(modelF);
 161         calcExpectedEffectiveLengths<ModelType>(M, refs, model, eel);
 162
 163         int num_threads = min(nThreads, nCV);
 164
 165         buffer = new Buffer(nMB, nSamples, M, l_bars, tmpF);
 166
 167         paramsArray = new Params[num_threads];
 168         threads = new pthread_t[num_threads];
 169
 170         char inpF[STRLEN];
 171         hasSeed ? engineFactory::init(seed) : engineFactory::init();
 172         for (int i = 0; i < num_threads; i++) {
 173                 paramsArray[i].no = i;
 174                 sprintf(inpF, "%s%d", cvsF, i);
 175                 paramsArray[i].fi = fopen(inpF, "r");
 176                 paramsArray[i].engine = engineFactory::new_engine();
 177                 paramsArray[i].mw = model.getMW();
 178         }
 179         engineFactory::finish();
 180
 181         /* set thread attribute to be joinable */
 182         pthread_attr_init(&attr);
 183         pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
 184
 185         for (int i = 0; i < num_threads; i++) {
 186                 rc = pthread_create(&threads[i], &attr, &sample_theta_from_c, (void*)(&paramsArray[i]));
 187                 pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0) in sample_theta_vectors_from_count_vectors!");
 188         }
 189         for (int i = 0; i < num_threads; i++) {
 190                 rc = pthread_join(threads[i], NULL);
 191                 pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0) in sample_theta_vectors_from_count_vectors!");
 192         }
 193
 194         /* destroy attribute */
 195         pthread_attr_destroy(&attr);
 196         delete[] threads;
 197
 198         for (int i = 0; i < num_threads; i++) {
 199                 fclose(paramsArray[i].fi);
 200                 delete paramsArray[i].engine;
 201         }
 202         delete[] paramsArray;
 203
 204         delete buffer; // Must delete here, force the content left in the buffer be written into the disk
 205
 206         if (verbose) { printf("Sampling is finished!\n"); }
 207 }
 208
 209 void calcCI(int nSamples, float *samples, float &lb, float &ub) {
 210         int p, q; // p pointer for lb, q pointer for ub;
 211         int newp, newq;
 212         int threshold = nSamples - (int(confidence * nSamples - 1e-8) + 1);
 213         int nOutside = 0;
 214
 215         sort(samples, samples + nSamples);
 216
 217         p = 0; q = nSamples - 1;
 218         newq = nSamples - 1;
 219         do {
 220                 q = newq;
 221                 while (newq > 0 && samples[newq - 1] == samples[newq]) newq--;
 222                 newq--;
 223         } while (newq >= 0 && nSamples - (newq + 1) <= threshold);
 224
 225         nOutside = nSamples - (q + 1);
 226
 227         lb = -1e30; ub = 1e30;
 228         do {
 229                 if (samples[q] - samples[p] < ub - lb) {
 230                         lb = samples[p];
 231                         ub = samples[q];
 232                 }
 233
 234                 newp = p;
 235                 while (newp < nSamples - 1 && samples[newp] == samples[newp + 1]) newp++;
 236                 newp++;
 237                 if (newp <= threshold) {
 238                         nOutside += newp - p;
 239                         p = newp;
 240                         while (nOutside > threshold && q < nSamples - 1) {
 241                                 newq = q + 1;
 242                                 while (newq < nSamples - 1 && samples[newq] == samples[newq + 1]) newq++;
 243                                 nOutside -= newq - q;
 244                                 q = newq;
 245                         }
 246                         assert(nOutside <= threshold);
 247                 }
 248                 else p = newp;
 249         } while (p <= threshold);
 250 }
 251
 252 void* calcCI_batch(void* arg) {
 253         float *tsamples, *fsamples;
 254         float *itsamples = NULL, *ifsamples = NULL, *gtsamples, *gfsamples;
 255         ifstream fin;
 256         CIParams *ciParams = (CIParams*)arg;
 257         int curtid, curaid, tid;
 258
 259         tsamples = new float[nSamples];
 260         fsamples = new float[nSamples];
 261         if (alleleS) {
 262           itsamples = new float[nSamples];
 263           ifsamples = new float[nSamples];
 264         }
 265         gtsamples = new float[nSamples];
 266         gfsamples = new float[nSamples];
 267
 268         fin.open(tmpF, ios::binary);
 269         // minus 1 here for that theta0 is not written!
 270         streampos pos = streampos(gi.spAt(ciParams->start_gene_id) - 1) * nSamples * FLOATSIZE;
 271         fin.seekg(pos, ios::beg);
 272
 273         int cnt = 0;
 274         if (alleleS) {
 275           curtid = curaid = -1;
 276           memset(itsamples, 0, FLOATSIZE * nSamples);
 277           memset(ifsamples, 0, FLOATSIZE * nSamples);
 278         }
 279         for (int i = ciParams->start_gene_id; i < ciParams->end_gene_id; i++) {
 280                 int b = gi.spAt(i), e = gi.spAt(i + 1);
 281                 memset(gtsamples, 0, FLOATSIZE * nSamples);
 282                 memset(gfsamples, 0, FLOATSIZE * nSamples);
 283                 for (int j = b; j < e; j++) {
 284                         if (alleleS) {
 285                           tid = ta.gidAt(j);
 286                           if (curtid != tid) {
 287                             if (curtid >= 0) {
 288                               if (j - curaid > 1) {
 289                                 calcCI(nSamples, itsamples, iso_tpm[curtid].lb, iso_tpm[curtid].ub);
 290                                 calcCI(nSamples, ifsamples, iso_fpkm[curtid].lb, iso_fpkm[curtid].ub);
 291                               }
 292                               else {
 293                                 iso_tpm[curtid] = tpm[curaid];
 294                                 iso_fpkm[curtid] = fpkm[curaid];
 295                               }
 296                             }
 297                             curtid = tid;
 298                             curaid = j;
 299                           }
 300                         }
 301
 302                         for (int k = 0; k < nSamples; k++) {
 303                                 fin.read((char*)(&tsamples[k]), FLOATSIZE);
 304                                 fsamples[k] = 1e3 / l_bars[k] * tsamples[k];
 305                                 if (alleleS) {
 306                                   itsamples[k] += tsamples[k];
 307                                   ifsamples[k] += fsamples[k];
 308                                 }
 309                                 gtsamples[k] += tsamples[k];
 310                                 gfsamples[k] += fsamples[k];
 311                         }
 312                         calcCI(nSamples, tsamples, tpm[j].lb, tpm[j].ub);
 313                         calcCI(nSamples, fsamples, fpkm[j].lb, fpkm[j].ub);
 314                 }
 315
 316                 if (e - b > 1) {
 317                         calcCI(nSamples, gtsamples, gene_tpm[i].lb, gene_tpm[i].ub);
 318                         calcCI(nSamples, gfsamples, gene_fpkm[i].lb, gene_fpkm[i].ub);
 319                 }
 320                 else {
 321                         gene_tpm[i] = tpm[b];
 322                         gene_fpkm[i] = fpkm[b];
 323                 }
 324
 325                 ++cnt;
 326                 if (verbose && cnt % 1000 == 0) { printf("In thread %d, %d genes are processed for CI calculation!\n", ciParams->no, cnt); }
 327         }
 328         fin.close();
 329
 330         if (alleleS && (curtid >= 0)) {
 331           if (gi.spAt(ciParams->end_gene_id) - curaid > 1) {
 332             calcCI(nSamples, itsamples, iso_tpm[curtid].lb, iso_tpm[curtid].ub);
 333             calcCI(nSamples, ifsamples, iso_fpkm[curtid].lb, iso_fpkm[curtid].ub);
 334           }
 335           else {
 336             iso_tpm[curtid] = tpm[curaid];
 337             iso_fpkm[curtid] = fpkm[curaid];
 338           }
 339         }
 340
 341         delete[] tsamples;
 342         delete[] fsamples;
 343         if (alleleS) {
 344           delete[] itsamples;
 345           delete[] ifsamples;
 346         }
 347         delete[] gtsamples;
 348         delete[] gfsamples;
 349
 350         return NULL;
 351 }
 352
 353 void calculate_credibility_intervals(char* imdName) {
 354         FILE *fo;
 355         char outF[STRLEN];
 356         int num_threads = nThreads;
 357
 358         tpm = new CIType[M + 1];
 359         fpkm = new CIType[M + 1];
 360         if (alleleS) {
 361           iso_tpm = new CIType[m_trans];
 362           iso_fpkm = new CIType[m_trans];
 363         }
 364         gene_tpm = new CIType[m];
 365         gene_fpkm = new CIType[m];
 366
 367         assert(M > 0);
 368         int quotient = M / num_threads;
 369         if (quotient < 1) { num_threads = M; quotient = 1; }
 370         int cur_gene_id = 0;
 371         int num_isoforms = 0;
 372
 373         // A just so so strategy for paralleling
 374         ciParamsArray = new CIParams[num_threads];
 375         for (int i = 0; i < num_threads; i++) {
 376                 ciParamsArray[i].no = i;
 377                 ciParamsArray[i].start_gene_id = cur_gene_id;
 378                 num_isoforms = 0;
 379
 380                 while ((m - cur_gene_id > num_threads - i - 1) && (i == num_threads - 1 || num_isoforms < quotient)) {
 381                         num_isoforms += gi.spAt(cur_gene_id + 1) - gi.spAt(cur_gene_id);
 382                         ++cur_gene_id;
 383                 }
 384
 385                 ciParamsArray[i].end_gene_id = cur_gene_id;
 386         }
 387
 388         threads = new pthread_t[num_threads];
 389
 390         /* set thread attribute to be joinable */
 391         pthread_attr_init(&attr);
 392         pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
 393
 394         // paralleling
 395         for (int i = 0; i < num_threads; i++) {
 396                 rc = pthread_create(&threads[i], &attr, &calcCI_batch, (void*)(&ciParamsArray[i]));
 397                 pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0) in calculate_credibility_intervals!");
 398         }
 399         for (int i = 0; i < num_threads; i++) {
 400                 rc = pthread_join(threads[i], NULL);
 401                 pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0) in calculate_credibility_intervals!");
 402         }
 403
 404         // releasing resources
 405
 406         /* destroy attribute */
 407         pthread_attr_destroy(&attr);
 408         delete[] threads;
 409
 410         delete[] ciParamsArray;
 411
 412         alleleS ? sprintf(outF, "%s.allele_res", imdName) : sprintf(outF, "%s.iso_res", imdName);
 413         fo = fopen(outF, "a");
 414         for (int i = 1; i <= M; i++)
 415           fprintf(fo, "%.6g%c", tpm[i].lb, (i < M ? '\t' : '\n'));
 416         for (int i = 1; i <= M; i++)
 417           fprintf(fo, "%.6g%c", tpm[i].ub, (i < M ? '\t' : '\n'));
 418         for (int i = 1; i <= M; i++)
 419           fprintf(fo, "%.6g%c", fpkm[i].lb, (i < M ? '\t' : '\n'));
 420         for (int i = 1; i <= M; i++)
 421           fprintf(fo, "%.6g%c", fpkm[i].ub, (i < M ? '\t' : '\n'));
 422         fclose(fo);
 423
 424         if (alleleS) {
 425           //isoform level results
 426           sprintf(outF, "%s.iso_res", imdName);
 427           fo = fopen(outF, "a");
 428           for (int i = 0; i < m_trans; i++)
 429             fprintf(fo, "%.6g%c", iso_tpm[i].lb, (i < m_trans - 1 ? '\t' : '\n'));
 430           for (int i = 0; i < m_trans; i++)
 431             fprintf(fo, "%.6g%c", iso_tpm[i].ub, (i < m_trans - 1 ? '\t' : '\n'));
 432           for (int i = 0; i < m_trans; i++)
 433             fprintf(fo, "%.6g%c", iso_fpkm[i].lb, (i < m_trans - 1 ? '\t' : '\n'));
 434           for (int i = 0; i < m_trans; i++)
 435             fprintf(fo, "%.6g%c", iso_fpkm[i].ub, (i < m_trans - 1 ? '\t' : '\n'));
 436           fclose(fo);
 437         }
 438
 439         //gene level results
 440         sprintf(outF, "%s.gene_res", imdName);
 441         fo = fopen(outF, "a");
 442         for (int i = 0; i < m; i++)
 443                 fprintf(fo, "%.6g%c", gene_tpm[i].lb, (i < m - 1 ? '\t' : '\n'));
 444         for (int i = 0; i < m; i++)
 445                 fprintf(fo, "%.6g%c", gene_tpm[i].ub, (i < m - 1 ? '\t' : '\n'));
 446         for (int i = 0; i < m; i++)
 447                 fprintf(fo, "%.6g%c", gene_fpkm[i].lb, (i < m - 1 ? '\t' : '\n'));
 448         for (int i = 0; i < m; i++)
 449                 fprintf(fo, "%.6g%c", gene_fpkm[i].ub, (i < m - 1 ? '\t' : '\n'));
 450         fclose(fo);
 451
 452         delete[] tpm;
 453         delete[] fpkm;
 454         if (alleleS) {
 455           delete[] iso_tpm;
 456           delete[] iso_fpkm;
 457         }
 458         delete[] gene_tpm;
 459         delete[] gene_fpkm;
 460
 461         if (verbose) { printf("All credibility intervals are calculated!\n"); }
 462 }
 463
 464 int main(int argc, char* argv[]) {
 465         if (argc < 8) {
 466                 printf("Usage: rsem-calculate-credibility-intervals reference_name imdName statName confidence nCV nSpC nMB [-p #Threads] [--seed seed] [-q]\n");
 467                 exit(-1);
 468         }
 469
 470         strcpy(refName, argv[1]);
 471         strcpy(imdName, argv[2]);
 472         strcpy(statName, argv[3]);
 473
 474         confidence = atof(argv[4]);
 475         nCV = atoi(argv[5]);
 476         nSpC = atoi(argv[6]);
 477         nMB = atoi(argv[7]);
 478
 479         nThreads = 1;
 480         quiet = false;
 481         hasSeed = false;
 482         for (int i = 8; i < argc; i++) {
 483                 if (!strcmp(argv[i], "-p")) nThreads = atoi(argv[i + 1]);
 484                 if (!strcmp(argv[i], "--seed")) {
 485                   hasSeed = true;
 486                   int len = strlen(argv[i + 1]);
 487                   seed = 0;
 488                   for (int k = 0; k < len; k++) seed = seed * 10 + (argv[i + 1][k] - '0');
 489                 }
 490                 if (!strcmp(argv[i], "-q")) quiet = true;
 491         }
 492         verbose = !quiet;
 493
 494         sprintf(refF, "%s.seq", refName);
 495         refs.loadRefs(refF, 1);
 496         M = refs.getM();
 497
 498         sprintf(groupF, "%s.grp", refName);
 499         gi.load(groupF);
 500         m = gi.getm();
 501
 502         // allele-specific
 503         alleleS = isAlleleSpecific(refName, NULL, &ta);
 504         if (alleleS) m_trans = ta.getm();
 505
 506         nSamples = nCV * nSpC;
 507         assert(nSamples > 0 && M > 0); // for Buffter.h: (bufsize_type)nSamples
 508         l_bars = new float[nSamples];
 509
 510         sprintf(tmpF, "%s.tmp", imdName);
 511         sprintf(cvsF, "%s.countvectors", imdName);
 512
 513         sprintf(modelF, "%s.model", statName);
 514         FILE *fi = fopen(modelF, "r");
 515         general_assert(fi != NULL, "Cannot open " + cstrtos(modelF) + "!");
 516         assert(fscanf(fi, "%d", &model_type) == 1);
 517         fclose(fi);
 518
 519         // Phase I
 520         switch(model_type) {
 521         case 0 : sample_theta_vectors_from_count_vectors<SingleModel>(); break;
 522         case 1 : sample_theta_vectors_from_count_vectors<SingleQModel>(); break;
 523         case 2 : sample_theta_vectors_from_count_vectors<PairedEndModel>(); break;
 524         case 3 : sample_theta_vectors_from_count_vectors<PairedEndQModel>(); break;
 525         }
 526
 527         // Phase II
 528         calculate_credibility_intervals(imdName);
 529
 530         delete l_bars;
 531
 532         return 0;
 533 }