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
  24 #include "Buffer.h"
  25 using namespace std;
  26
  27 struct Params {
  28         int no;
  29         FILE *fi;
  30         engine_type *engine;
  31         double *mw;
  32 };
  33
  34 struct CIType {
  35         float lb, ub; // the interval is [lb, ub]
  36
  37         CIType() { lb = ub = 0.0; }
  38 };
  39
  40 struct CIParams {
  41         int no;
  42         int start_gene_id, end_gene_id;
  43 };
  44
  45 int model_type;
  46
  47 int nMB;
  48 double confidence;
  49 int nCV, nSpC, nSamples; // nCV: number of count vectors; nSpC: number of theta vectors sampled per count vector; nSamples: nCV * nSpC
  50 int nThreads;
  51 int cvlen;
  52
  53 char cvsF[STRLEN], tmpF[STRLEN], command[STRLEN];
  54
  55 CIType *iso_tau, *gene_tau;
  56
  57 int M, m;
  58 Refs refs;
  59 GroupInfo gi;
  60 char imdName[STRLEN], statName[STRLEN];
  61 char modelF[STRLEN], groupF[STRLEN], refF[STRLEN];
  62
  63 vector<double> eel; //expected effective lengths
  64
  65 Buffer *buffer;
  66
  67 bool quiet;
  68
  69 Params *paramsArray;
  70 pthread_t *threads;
  71 pthread_attr_t attr;
  72 void *status;
  73 int rc;
  74
  75 CIParams *ciParamsArray;
  76
  77 template<class ModelType>
  78 void calcExpectedEffectiveLengths(ModelType& model) {
  79         int lb, ub, span;
  80         double *pdf = NULL, *cdf = NULL, *clen = NULL; // clen[i] = \sigma_{j=1}^{i}pdf[i]*(lb+i)
  81
  82         model.getGLD().copyTo(pdf, cdf, lb, ub, span);
  83         clen = new double[span + 1];
  84         clen[0] = 0.0;
  85         for (int i = 1; i <= span; i++) {
  86                 clen[i] = clen[i - 1] + pdf[i] * (lb + i);
  87         }
  88
  89         eel.assign(M + 1, 0.0);
  90         for (int i = 1; i <= M; i++) {
  91                 int totLen = refs.getRef(i).getTotLen();
  92                 int fullLen = refs.getRef(i).getFullLen();
  93                 int pos1 = max(min(totLen - fullLen + 1, ub) - lb, 0);
  94                 int pos2 = max(min(totLen, ub) - lb, 0);
  95
  96                 if (pos2 == 0) { eel[i] = 0.0; continue; }
  97
  98                 eel[i] = fullLen * cdf[pos1] + ((cdf[pos2] - cdf[pos1]) * (totLen + 1) - (clen[pos2] - clen[pos1]));
  99                 assert(eel[i] >= 0);
 100                 if (eel[i] < MINEEL) { eel[i] = 0.0; }
 101         }
 102
 103         delete[] pdf;
 104         delete[] cdf;
 105         delete[] clen;
 106 }
 107
 108 void* sample_theta_from_c(void* arg) {
 109
 110         int *cvec;
 111         double *theta;
 112         gamma_dist **gammas;
 113         gamma_generator **rgs;
 114
 115         Params *params = (Params*)arg;
 116         FILE *fi = params->fi;
 117         double *mw = params->mw;
 118
 119         cvec = new int[cvlen];
 120         theta = new double[cvlen];
 121         gammas = new gamma_dist*[cvlen];
 122         rgs = new gamma_generator*[cvlen];
 123
 124         float **vecs = new float*[nSpC];
 125         for (int i = 0; i < nSpC; i++) vecs[i] = new float[cvlen];
 126
 127         int cnt = 0;
 128         while (fscanf(fi, "%d", &cvec[0]) == 1) {
 129                 for (int j = 1; j < cvlen; j++) assert(fscanf(fi, "%d", &cvec[j]) == 1);
 130
 131                 ++cnt;
 132
 133                 for (int j = 0; j < cvlen; j++) {
 134                         gammas[j] = new gamma_dist(cvec[j]);
 135                         rgs[j] = new gamma_generator(*(params->engine), *gammas[j]);
 136                 }
 137
 138                 for (int i = 0; i < nSpC; i++) {
 139                         double sum = 0.0;
 140                         for (int j = 0; j < cvlen; j++) {
 141                                 theta[j] = ((j == 0 || eel[j] >= EPSILON) ? (*rgs[j])() : 0.0);
 142                                 sum += theta[j];
 143                         }
 144                         assert(sum >= EPSILON);
 145                         for (int j = 0; j < cvlen; j++) theta[j] /= sum;
 146
 147                         sum = 0.0;
 148                         for (int j = 0; j < cvlen; j++) {
 149                                 theta[j] = (mw[j] < EPSILON ? 0.0 : theta[j] / mw[j]);
 150                                 sum += theta[j];
 151                         }
 152                         assert(sum >= EPSILON);
 153                         for (int j = 0; j < cvlen; j++) theta[j] /= sum;
 154
 155
 156                         sum = 0.0;
 157                         vecs[i][0] = theta[0];
 158                         for (int j = 1; j < cvlen; j++)
 159                                 if (eel[j] >= EPSILON) {
 160                                         vecs[i][j] = theta[j] / eel[j];
 161                                         sum += vecs[i][j];
 162                                 }
 163                                 else assert(theta[j] < EPSILON);
 164
 165                         assert(sum >= EPSILON);
 166                         for (int j = 1; j < cvlen; j++) vecs[i][j] /= sum;
 167                 }
 168
 169                 buffer->write(nSpC, vecs);
 170
 171                 for (int j = 0; j < cvlen; j++) {
 172                         delete gammas[j];
 173                         delete rgs[j];
 174                 }
 175
 176                 if (verbose && cnt % 100 == 0) { printf("Thread %d, %d count vectors are processed!\n", params->no, cnt); }
 177         }
 178
 179         delete[] cvec;
 180         delete[] theta;
 181         delete[] gammas;
 182         delete[] rgs;
 183
 184         for (int i = 0; i < nSpC; i++) delete[] vecs[i];
 185         delete[] vecs;
 186
 187         return NULL;
 188 }
 189
 190 template<class ModelType>
 191 void sample_theta_vectors_from_count_vectors() {
 192         ModelType model;
 193         model.read(modelF);
 194         calcExpectedEffectiveLengths<ModelType>(model);
 195
 196         buffer = new Buffer(nMB, nSamples, cvlen, tmpF);
 197
 198         paramsArray = new Params[nThreads];
 199         threads = new pthread_t[nThreads];
 200
 201         char inpF[STRLEN];
 202         for (int i = 0; i < nThreads; i++) {
 203                 paramsArray[i].no = i;
 204                 sprintf(inpF, "%s%d", cvsF, i);
 205                 paramsArray[i].fi = fopen(inpF, "r");
 206                 paramsArray[i].engine = engineFactory::new_engine();
 207                 paramsArray[i].mw = model.getMW();
 208         }
 209
 210         /* set thread attribute to be joinable */
 211         pthread_attr_init(&attr);
 212         pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
 213
 214         for (int i = 0; i < nThreads; i++) {
 215                 rc = pthread_create(&threads[i], &attr, &sample_theta_from_c, (void*)(&paramsArray[i]));
 216                 pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0) in sample_theta_vectors_from_count_vectors!");
 217         }
 218         for (int i = 0; i < nThreads; i++) {
 219                 rc = pthread_join(threads[i], &status);
 220                 pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0) in sample_theta_vectors_from_count_vectors!");
 221         }
 222
 223         /* destroy attribute */
 224         pthread_attr_destroy(&attr);
 225         delete[] threads;
 226
 227         for (int i = 0; i < nThreads; i++) {
 228                 fclose(paramsArray[i].fi);
 229                 delete paramsArray[i].engine;
 230         }
 231         delete[] paramsArray;
 232
 233         delete buffer; // Must delete here, force the content left in the buffer be written into the disk
 234
 235         if (verbose) { printf("Sampling is finished!\n"); }
 236 }
 237
 238 void calcCI(int nSamples, float *samples, float &lb, float &ub) {
 239         int p, q; // p pointer for lb, q pointer for ub;
 240         int newp, newq;
 241         int threshold = nSamples - (int(confidence * nSamples - 1e-8) + 1);
 242         int nOutside = 0;
 243
 244         sort(samples, samples + nSamples);
 245
 246         p = 0; q = nSamples - 1;
 247         newq = nSamples - 1;
 248         do {
 249                 q = newq;
 250                 while (newq > 0 && samples[newq - 1] == samples[newq]) newq--;
 251                 newq--;
 252         } while (newq >= 0 && nSamples - (newq + 1) <= threshold);
 253
 254         nOutside = nSamples - (q + 1);
 255
 256         lb = -1e30; ub = 1e30;
 257         do {
 258                 if (samples[q] - samples[p] < ub - lb) {
 259                         lb = samples[p];
 260                         ub = samples[q];
 261                 }
 262
 263                 newp = p;
 264                 while (newp < nSamples - 1 && samples[newp] == samples[newp + 1]) newp++;
 265                 newp++;
 266                 if (newp <= threshold) {
 267                         nOutside += newp - p;
 268                         p = newp;
 269                         while (nOutside > threshold && q < nSamples - 1) {
 270                                 newq = q + 1;
 271                                 while (newq < nSamples - 1 && samples[newq] == samples[newq + 1]) newq++;
 272                                 nOutside -= newq - q;
 273                                 q = newq;
 274                         }
 275                         assert(nOutside <= threshold);
 276                 }
 277                 else p = newp;
 278         } while (p <= threshold);
 279 }
 280
 281 void* calcCI_batch(void* arg) {
 282         float *itsamples, *gtsamples;
 283         ifstream fin;
 284         CIParams *ciParams = (CIParams*)arg;
 285
 286         itsamples = new float[nSamples];
 287         gtsamples = new float[nSamples];
 288
 289         fin.open(tmpF, ios::binary);
 290         streampos pos = streampos(gi.spAt(ciParams->start_gene_id)) * nSamples * FLOATSIZE;
 291         fin.seekg(pos, ios::beg);
 292
 293         int cnt = 0;
 294         for (int i = ciParams->start_gene_id; i < ciParams->end_gene_id; i++) {
 295                 int b = gi.spAt(i), e = gi.spAt(i + 1);
 296                 memset(gtsamples, 0, FLOATSIZE * nSamples);
 297                 for (int j = b; j < e; j++) {
 298                         for (int k = 0; k < nSamples; k++) {
 299                                 fin.read((char*)(&itsamples[k]), FLOATSIZE);
 300                                 gtsamples[k] += itsamples[k];
 301                         }
 302                         calcCI(nSamples, itsamples, iso_tau[j].lb, iso_tau[j].ub);
 303                 }
 304                 calcCI(nSamples, gtsamples, gene_tau[i].lb, gene_tau[i].ub);
 305
 306                 ++cnt;
 307                 if (verbose && cnt % 1000 == 0) { printf("In thread %d, %d genes are processed for CI calculation!\n", ciParams->no, cnt); }
 308         }
 309
 310         fin.close();
 311
 312         delete[] itsamples;
 313         delete[] gtsamples;
 314
 315         return NULL;
 316 }
 317
 318 void calculate_credibility_intervals(char* imdName) {
 319         FILE *fo;
 320         char outF[STRLEN];
 321
 322         iso_tau = new CIType[M + 1];
 323         gene_tau = new CIType[m];
 324
 325         assert(M > 0);
 326         int quotient = M / nThreads;
 327         if (quotient < 1) { nThreads = M; quotient = 1; }
 328         int cur_gene_id = 0;
 329         int num_isoforms = 0;
 330
 331         // A just so so strategy for paralleling
 332         ciParamsArray = new CIParams[nThreads];
 333         for (int i = 0; i < nThreads; i++) {
 334                 ciParamsArray[i].no = i;
 335                 ciParamsArray[i].start_gene_id = cur_gene_id;
 336                 num_isoforms = 0;
 337
 338                 while ((m - cur_gene_id > nThreads - i - 1) && (i == nThreads - 1 || num_isoforms < quotient)) {
 339                         num_isoforms += gi.spAt(cur_gene_id + 1) - gi.spAt(cur_gene_id);
 340                         ++cur_gene_id;
 341                 }
 342
 343                 ciParamsArray[i].end_gene_id = cur_gene_id;
 344         }
 345
 346         threads = new pthread_t[nThreads];
 347
 348         /* set thread attribute to be joinable */
 349         pthread_attr_init(&attr);
 350         pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
 351
 352         // paralleling
 353         for (int i = 0; i < nThreads; i++) {
 354                 rc = pthread_create(&threads[i], &attr, &calcCI_batch, (void*)(&ciParamsArray[i]));
 355                 pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0) in calculate_credibility_intervals!");
 356         }
 357         for (int i = 0; i < nThreads; i++) {
 358                 rc = pthread_join(threads[i], &status);
 359                 pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0) in calculate_credibility_intervals!");
 360         }
 361
 362         // releasing resources
 363
 364         /* destroy attribute */
 365         pthread_attr_destroy(&attr);
 366         delete[] threads;
 367
 368         delete[] ciParamsArray;
 369
 370         //isoform level results
 371         sprintf(outF, "%s.iso_res", imdName);
 372         fo = fopen(outF, "a");
 373         for (int i = 1; i <= M; i++)
 374                 fprintf(fo, "%.6g%c", iso_tau[i].lb, (i < M ? '\t' : '\n'));
 375         for (int i = 1; i <= M; i++)
 376                 fprintf(fo, "%.6g%c", iso_tau[i].ub, (i < M ? '\t' : '\n'));
 377         fclose(fo);
 378
 379         //gene level results
 380         sprintf(outF, "%s.gene_res", imdName);
 381         fo = fopen(outF, "a");
 382         for (int i = 0; i < m; i++)
 383                 fprintf(fo, "%.6g%c", gene_tau[i].lb, (i < m - 1 ? '\t' : '\n'));
 384         for (int i = 0; i < m; i++)
 385                 fprintf(fo, "%.6g%c", gene_tau[i].ub, (i < m - 1 ? '\t' : '\n'));
 386         fclose(fo);
 387
 388         delete[] iso_tau;
 389         delete[] gene_tau;
 390
 391         if (verbose) { printf("All credibility intervals are calculated!\n"); }
 392 }
 393
 394 int main(int argc, char* argv[]) {
 395         if (argc < 8) {
 396                 printf("Usage: rsem-calculate-credibility-intervals reference_name sample_name sampleToken confidence nCV nSpC nMB [-p #Threads] [-q]\n");
 397                 exit(-1);
 398         }
 399
 400         confidence = atof(argv[4]);
 401         nCV = atoi(argv[5]);
 402         nSpC = atoi(argv[6]);
 403         nMB = atoi(argv[7]);
 404
 405         nThreads = 1;
 406         quiet = false;
 407         for (int i = 8; i < argc; i++) {
 408                 if (!strcmp(argv[i], "-p")) nThreads = atoi(argv[i + 1]);
 409                 if (!strcmp(argv[i], "-q")) quiet = true;
 410         }
 411         verbose = !quiet;
 412
 413         if (nThreads > nCV) {
 414                 nThreads = nCV;
 415                 printf("Warning: Number of count vectors is less than number of threads! Change the number of threads to %d!\n", nThreads);
 416         }
 417
 418         sprintf(refF, "%s.seq", argv[1]);
 419         refs.loadRefs(refF, 1);
 420         M = refs.getM();
 421         sprintf(groupF, "%s.grp", argv[1]);
 422         gi.load(groupF);
 423         m = gi.getm();
 424
 425         nSamples = nCV * nSpC;
 426         cvlen = M + 1;
 427         assert(nSamples > 0 && cvlen > 1); // for Buffter.h: (bufsize_type)nSamples
 428
 429         sprintf(imdName, "%s.temp/%s", argv[2], argv[3]);
 430         sprintf(statName, "%s.stat/%s", argv[2], argv[3]);
 431         sprintf(tmpF, "%s.tmp", imdName);
 432         sprintf(cvsF, "%s.countvectors", imdName);
 433
 434         sprintf(modelF, "%s.model", statName);
 435         FILE *fi = fopen(modelF, "r");
 436         general_assert(fi != NULL, "Cannot open " + cstrtos(modelF) + "!");
 437         assert(fscanf(fi, "%d", &model_type) == 1);
 438         fclose(fi);
 439
 440         // Phase I
 441         switch(model_type) {
 442         case 0 : sample_theta_vectors_from_count_vectors<SingleModel>(); break;
 443         case 1 : sample_theta_vectors_from_count_vectors<SingleQModel>(); break;
 444         case 2 : sample_theta_vectors_from_count_vectors<PairedEndModel>(); break;
 445         case 3 : sample_theta_vectors_from_count_vectors<PairedEndQModel>(); break;
 446         }
 447
 448         // Phase II
 449         calculate_credibility_intervals(imdName);
 450
 451         /*
 452         sprintf(command, "rm -f %s", tmpF);
 453         int status = system(command);
 454         if (status != 0) {
 455                 fprintf(stderr, "Cannot delete %s!\n", tmpF);
 456                 exit(-1);
 457         }
 458         */
 459
 460         return 0;
 461 }