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