12 #include "my_assert.h"
16 #include "SingleModel.h"
17 #include "SingleQModel.h"
18 #include "PairedEndModel.h"
19 #include "PairedEndQModel.h"
22 #include "GroupInfo.h"
23 #include "WriteResults.h"
38 int start_gene_id, end_gene_id;
42 float lb, ub; // the interval is [lb, ub]
44 CIType() { lb = ub = 0.0; }
51 int nCV, nSpC, nSamples; // nCV: number of count vectors; nSpC: number of theta vectors sampled per count vector; nSamples: nCV * nSpC
56 char cvsF[STRLEN], tmpF[STRLEN], command[STRLEN];
59 CIType *iso_tpm = NULL, *iso_fpkm = NULL;
60 CIType *gene_tpm, *gene_fpkm;
65 char refName[STRLEN], imdName[STRLEN], statName[STRLEN];
66 char modelF[STRLEN], groupF[STRLEN], refF[STRLEN];
72 vector<double> eel; //expected effective lengths
83 CIParams *ciParamsArray;
85 void* sample_theta_from_c(void* arg) {
90 gamma_generator **rgs;
92 Params *params = (Params*)arg;
93 FILE *fi = params->fi;
94 const double *mw = params->mw;
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
104 while (fscanf(fi, "%d", &cvec[0]) == 1) {
105 for (int j = 1; j <= M; j++) assert(fscanf(fi, "%d", &cvec[j]) == 1);
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]);
114 for (int i = 0; i < nSpC; i++) {
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);
120 assert(sum >= EPSILON);
121 for (int j = 0; j <= M; j++) theta[j] /= sum;
125 for (int j = 1; j <= M; j++)
126 if (eel[j] >= EPSILON) {
127 tpm[j] = theta[j] / eel[j];
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
137 for (int j = 0; j <= M; j++) {
142 if (verbose && cnt % 100 == 0) { printf("Thread %d, %d count vectors are processed!\n", params->no, cnt); }
154 template<class ModelType>
155 void sample_theta_vectors_from_count_vectors() {
158 calcExpectedEffectiveLengths<ModelType>(M, refs, model, eel);
160 int num_threads = min(nThreads, nCV);
162 buffer = new Buffer(nMB, nSamples, M, l_bars, tmpF);
164 paramsArray = new Params[num_threads];
165 threads = new pthread_t[num_threads];
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();
176 /* set thread attribute to be joinable */
177 pthread_attr_init(&attr);
178 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
180 for (int i = 0; i < num_threads; i++) {
181 rc = pthread_create(&threads[i], &attr, &sample_theta_from_c, (void*)(¶msArray[i]));
182 pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0) in sample_theta_vectors_from_count_vectors!");
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!");
189 /* destroy attribute */
190 pthread_attr_destroy(&attr);
193 for (int i = 0; i < num_threads; i++) {
194 fclose(paramsArray[i].fi);
195 delete paramsArray[i].engine;
197 delete[] paramsArray;
199 delete buffer; // Must delete here, force the content left in the buffer be written into the disk
201 if (verbose) { printf("Sampling is finished!\n"); }
204 void calcCI(int nSamples, float *samples, float &lb, float &ub) {
205 int p, q; // p pointer for lb, q pointer for ub;
207 int threshold = nSamples - (int(confidence * nSamples - 1e-8) + 1);
210 sort(samples, samples + nSamples);
212 p = 0; q = nSamples - 1;
216 while (newq > 0 && samples[newq - 1] == samples[newq]) newq--;
218 } while (newq >= 0 && nSamples - (newq + 1) <= threshold);
220 nOutside = nSamples - (q + 1);
222 lb = -1e30; ub = 1e30;
224 if (samples[q] - samples[p] < ub - lb) {
230 while (newp < nSamples - 1 && samples[newp] == samples[newp + 1]) newp++;
232 if (newp <= threshold) {
233 nOutside += newp - p;
235 while (nOutside > threshold && q < nSamples - 1) {
237 while (newq < nSamples - 1 && samples[newq] == samples[newq + 1]) newq++;
238 nOutside -= newq - q;
241 assert(nOutside <= threshold);
244 } while (p <= threshold);
247 void* calcCI_batch(void* arg) {
248 float *tsamples, *fsamples;
249 float *itsamples = NULL, *ifsamples = NULL, *gtsamples, *gfsamples;
251 CIParams *ciParams = (CIParams*)arg;
252 int curtid, curaid, tid;
254 tsamples = new float[nSamples];
255 fsamples = new float[nSamples];
257 itsamples = new float[nSamples];
258 ifsamples = new float[nSamples];
260 gtsamples = new float[nSamples];
261 gfsamples = new float[nSamples];
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);
270 curtid = curaid = -1;
271 memset(itsamples, 0, FLOATSIZE * nSamples);
272 memset(ifsamples, 0, FLOATSIZE * nSamples);
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++) {
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);
288 iso_tpm[curtid] = tpm[curaid];
289 iso_fpkm[curtid] = fpkm[curaid];
297 for (int k = 0; k < nSamples; k++) {
298 fin.read((char*)(&tsamples[k]), FLOATSIZE);
299 fsamples[k] = 1e3 / l_bars[k] * tsamples[k];
301 itsamples[k] += tsamples[k];
302 ifsamples[k] += fsamples[k];
304 gtsamples[k] += tsamples[k];
305 gfsamples[k] += fsamples[k];
307 calcCI(nSamples, tsamples, tpm[j].lb, tpm[j].ub);
308 calcCI(nSamples, fsamples, fpkm[j].lb, fpkm[j].ub);
312 calcCI(nSamples, gtsamples, gene_tpm[i].lb, gene_tpm[i].ub);
313 calcCI(nSamples, gfsamples, gene_fpkm[i].lb, gene_fpkm[i].ub);
316 gene_tpm[i] = tpm[b];
317 gene_fpkm[i] = fpkm[b];
321 if (verbose && cnt % 1000 == 0) { printf("In thread %d, %d genes are processed for CI calculation!\n", ciParams->no, cnt); }
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);
331 iso_tpm[curtid] = tpm[curaid];
332 iso_fpkm[curtid] = fpkm[curaid];
348 void calculate_credibility_intervals(char* imdName) {
351 int num_threads = nThreads;
353 tpm = new CIType[M + 1];
354 fpkm = new CIType[M + 1];
356 iso_tpm = new CIType[m_trans];
357 iso_fpkm = new CIType[m_trans];
359 gene_tpm = new CIType[m];
360 gene_fpkm = new CIType[m];
363 int quotient = M / num_threads;
364 if (quotient < 1) { num_threads = M; quotient = 1; }
366 int num_isoforms = 0;
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;
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);
380 ciParamsArray[i].end_gene_id = cur_gene_id;
383 threads = new pthread_t[num_threads];
385 /* set thread attribute to be joinable */
386 pthread_attr_init(&attr);
387 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
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!");
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!");
399 // releasing resources
401 /* destroy attribute */
402 pthread_attr_destroy(&attr);
405 delete[] ciParamsArray;
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'));
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'));
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'));
456 if (verbose) { printf("All credibility intervals are calculated!\n"); }
459 int main(int argc, char* argv[]) {
461 printf("Usage: rsem-calculate-credibility-intervals reference_name imdName statName confidence nCV nSpC nMB [-p #Threads] [-q]\n");
465 strcpy(refName, argv[1]);
466 strcpy(imdName, argv[2]);
467 strcpy(statName, argv[3]);
469 confidence = atof(argv[4]);
471 nSpC = atoi(argv[6]);
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;
482 sprintf(refF, "%s.seq", refName);
483 refs.loadRefs(refF, 1);
486 sprintf(groupF, "%s.grp", refName);
491 alleleS = isAlleleSpecific(refName, NULL, &ta);
492 if (alleleS) m_trans = ta.getm();
494 nSamples = nCV * nSpC;
495 assert(nSamples > 0 && M > 0); // for Buffter.h: (bufsize_type)nSamples
496 l_bars = new float[nSamples];
498 sprintf(tmpF, "%s.tmp", imdName);
499 sprintf(cvsF, "%s.countvectors", imdName);
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);
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;
516 calculate_credibility_intervals(imdName);