11 #include "my_assert.h"
15 #include "SingleModel.h"
16 #include "SingleQModel.h"
17 #include "PairedEndModel.h"
18 #include "PairedEndQModel.h"
22 #include "GroupInfo.h"
23 #include "WriteResults.h"
31 double *pme_c, *pve_c; //posterior mean and variance vectors on counts
32 double *pme_tpm, *pme_fpkm;
34 double *pve_c_genes, *pve_c_trans;
41 Item(int sid, double conprb) {
43 this->conprb = conprb;
54 int BURNIN, NSAMPLES, GAP;
55 char refName[STRLEN], imdName[STRLEN], statName[STRLEN];
56 char thetaF[STRLEN], ofgF[STRLEN], refF[STRLEN], modelF[STRLEN];
61 vector<HIT_INT_TYPE> s;
67 vector<int> pseudo_counts;
68 vector<double> pme_c, pve_c; //global posterior mean and variance vectors on counts
69 vector<double> pme_tpm, pme_fpkm;
88 vector<double> pve_c_genes, pve_c_trans;
90 void load_data(char* refName, char* statName, char* imdName) {
96 sprintf(refF, "%s.seq", refName);
97 refs.loadRefs(refF, 1);
101 sprintf(ofgF, "%s.ofg", imdName);
103 general_assert(fin.is_open(), "Cannot open " + cstrtos(ofgF) + "!");
105 general_assert(tmpVal == M, "M in " + cstrtos(ofgF) + " is not consistent with " + cstrtos(refF) + "!");
108 s.clear(); hits.clear();
110 while (getline(fin, line)) {
111 istringstream strin(line);
115 while (strin>>sid>>conprb) {
116 hits.push_back(Item(sid, conprb));
118 s.push_back(hits.size());
125 totc = N0 + N1 + (M + 1);
127 if (verbose) { printf("Loading data is finished!\n"); }
130 void load_group_info(char* refName) {
132 sprintf(groupF, "%s.grp", refName);
136 alleleS = isAlleleSpecific(refName, >, &ta); // if allele-specific
137 m_trans = (alleleS ? ta.getm() : 0);
139 if (verbose) { printf("Loading group information is finished!\n"); }
142 // Load imdName.omit and initialize the pseudo count vector.
143 void load_omit_info(const char* imdName) {
145 sprintf(omitF, "%s.omit", imdName);
146 FILE *fi = fopen(omitF, "r");
147 pseudo_counts.assign(M + 1, 1);
149 while (fscanf(fi, "%d", &tid) == 1) pseudo_counts[tid] = 0;
153 template<class ModelType>
154 void init_model_related(char* modelF) {
158 calcExpectedEffectiveLengths<ModelType>(M, refs, model, eel);
159 memcpy(mw, model.getMW(), sizeof(double) * (M + 1)); // otherwise, after exiting this procedure, mw becomes undefined
167 quotient = NSAMPLES / nThreads;
168 left = NSAMPLES % nThreads;
170 sprintf(cvsF, "%s.countvectors", imdName);
171 paramsArray = new Params[nThreads];
172 threads = new pthread_t[nThreads];
174 hasSeed ? engineFactory::init(seed) : engineFactory::init();
175 for (int i = 0; i < nThreads; i++) {
176 paramsArray[i].no = i;
178 paramsArray[i].nsamples = quotient;
179 if (i < left) paramsArray[i].nsamples++;
181 sprintf(outF, "%s%d", cvsF, i);
182 paramsArray[i].fo = fopen(outF, "w");
184 paramsArray[i].engine = engineFactory::new_engine();
185 paramsArray[i].pme_c = new double[M + 1];
186 memset(paramsArray[i].pme_c, 0, sizeof(double) * (M + 1));
187 paramsArray[i].pve_c = new double[M + 1];
188 memset(paramsArray[i].pve_c, 0, sizeof(double) * (M + 1));
189 paramsArray[i].pme_tpm = new double[M + 1];
190 memset(paramsArray[i].pme_tpm, 0, sizeof(double) * (M + 1));
191 paramsArray[i].pme_fpkm = new double[M + 1];
192 memset(paramsArray[i].pme_fpkm, 0, sizeof(double) * (M + 1));
194 paramsArray[i].pve_c_genes = new double[m];
195 memset(paramsArray[i].pve_c_genes, 0, sizeof(double) * m);
197 paramsArray[i].pve_c_trans = NULL;
199 paramsArray[i].pve_c_trans = new double[m_trans];
200 memset(paramsArray[i].pve_c_trans, 0, sizeof(double) * m_trans);
203 engineFactory::finish();
205 /* set thread attribute to be joinable */
206 pthread_attr_init(&attr);
207 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
209 if (verbose) { printf("Initialization finished!\n"); }
212 //sample theta from Dir(1)
213 void sampleTheta(engine_type& engine, vector<double>& theta) {
215 gamma_generator gmg(engine, gm);
218 theta.assign(M + 1, 0);
220 for (int i = 0; i <= M; i++) {
221 theta[i] = (pseudo_counts[i] > 0 ? gmg() : 0.0);
224 assert(denom > EPSILON);
225 for (int i = 0; i <= M; i++) theta[i] /= denom;
228 void writeCountVector(FILE* fo, vector<int>& counts) {
229 for (int i = 0; i < M; i++) {
230 fprintf(fo, "%d ", counts[i]);
232 fprintf(fo, "%d\n", counts[M]);
235 void* Gibbs(void* arg) {
237 HIT_INT_TYPE len, fr, to;
238 Params *params = (Params*)arg;
240 vector<double> theta, tpm, fpkm;
241 vector<int> z, counts(pseudo_counts);
244 uniform_01_generator rg(*params->engine, uniform_01_dist());
246 // generate initial state
247 sampleTheta(*params->engine, theta);
252 for (READ_INT_TYPE i = 0; i < N1; i++) {
253 fr = s[i]; to = s[i + 1];
256 for (HIT_INT_TYPE j = fr; j < to; j++) {
257 arr[j - fr] = theta[hits[j].sid] * hits[j].conprb;
258 if (j > fr) arr[j - fr] += arr[j - fr - 1]; // cumulative
260 z[i] = hits[fr + sample(rg, arr, len)].sid;
265 CHAINLEN = 1 + (params->nsamples - 1) * GAP;
266 for (int ROUND = 1; ROUND <= BURNIN + CHAINLEN; ROUND++) {
268 for (READ_INT_TYPE i = 0; i < N1; i++) {
270 fr = s[i]; to = s[i + 1]; len = to - fr;
272 for (HIT_INT_TYPE j = fr; j < to; j++) {
273 arr[j - fr] = counts[hits[j].sid] * hits[j].conprb;
274 if (j > fr) arr[j - fr] += arr[j - fr - 1]; //cumulative
276 z[i] = hits[fr + sample(rg, arr, len)].sid;
280 if (ROUND > BURNIN) {
281 if ((ROUND - BURNIN - 1) % GAP == 0) {
282 writeCountVector(params->fo, counts);
283 for (int i = 0; i <= M; i++) theta[i] = counts[i] / totc;
284 polishTheta(M, theta, eel, mw);
285 calcExpressionValues(M, theta, eel, tpm, fpkm);
286 for (int i = 0; i <= M; i++) {
287 params->pme_c[i] += counts[i] - pseudo_counts[i];
288 params->pve_c[i] += double(counts[i] - pseudo_counts[i]) * (counts[i] - pseudo_counts[i]);
289 params->pme_tpm[i] += tpm[i];
290 params->pme_fpkm[i] += fpkm[i];
293 for (int i = 0; i < m; i++) {
294 int b = gi.spAt(i), e = gi.spAt(i + 1);
296 for (int j = b; j < e; j++) count += counts[j] - pseudo_counts[j];
297 params->pve_c_genes[i] += count * count;
301 for (int i = 0; i < m_trans; i++) {
302 int b = ta.spAt(i), e = ta.spAt(i + 1);
304 for (int j = b; j < e; j++) count += counts[j] - pseudo_counts[j];
305 params->pve_c_trans[i] += count * count;
310 if (verbose && ROUND % 100 == 0) { printf("Thread %d, ROUND %d is finished!\n", params->no, ROUND); }
317 // char inpF[STRLEN], command[STRLEN];
320 /* destroy attribute */
321 pthread_attr_destroy(&attr);
324 pme_c.assign(M + 1, 0);
325 pve_c.assign(M + 1, 0);
326 pme_tpm.assign(M + 1, 0);
327 pme_fpkm.assign(M + 1, 0);
329 pve_c_genes.assign(m, 0);
331 if (alleleS) pve_c_trans.assign(m_trans, 0);
333 for (int i = 0; i < nThreads; i++) {
334 fclose(paramsArray[i].fo);
335 delete paramsArray[i].engine;
336 for (int j = 0; j <= M; j++) {
337 pme_c[j] += paramsArray[i].pme_c[j];
338 pve_c[j] += paramsArray[i].pve_c[j];
339 pme_tpm[j] += paramsArray[i].pme_tpm[j];
340 pme_fpkm[j] += paramsArray[i].pme_fpkm[j];
343 for (int j = 0; j < m; j++)
344 pve_c_genes[j] += paramsArray[i].pve_c_genes[j];
347 for (int j = 0; j < m_trans; j++)
348 pve_c_trans[j] += paramsArray[i].pve_c_trans[j];
350 delete[] paramsArray[i].pme_c;
351 delete[] paramsArray[i].pve_c;
352 delete[] paramsArray[i].pme_tpm;
353 delete[] paramsArray[i].pme_fpkm;
355 delete[] paramsArray[i].pve_c_genes;
356 if (alleleS) delete[] paramsArray[i].pve_c_trans;
358 delete[] paramsArray;
360 for (int i = 0; i <= M; i++) {
361 pme_c[i] /= NSAMPLES;
362 pve_c[i] = (pve_c[i] - double(NSAMPLES) * pme_c[i] * pme_c[i]) / double(NSAMPLES - 1);
363 if (pve_c[i] < 0.0) pve_c[i] = 0.0;
364 pme_tpm[i] /= NSAMPLES;
365 pme_fpkm[i] /= NSAMPLES;
368 for (int i = 0; i < m; i++) {
369 int b = gi.spAt(i), e = gi.spAt(i + 1);
370 double pme_c_gene = 0.0;
371 for (int j = b; j < e; j++) pme_c_gene += pme_c[j];
372 pve_c_genes[i] = (pve_c_genes[i] - double(NSAMPLES) * pme_c_gene * pme_c_gene) / double(NSAMPLES - 1);
373 if (pve_c_genes[i] < 0.0) pve_c_genes[i] = 0.0;
377 for (int i = 0; i < m_trans; i++) {
378 int b = ta.spAt(i), e = ta.spAt(i + 1);
379 double pme_c_tran = 0.0;
380 for (int j = b; j < e; j++) pme_c_tran += pme_c[j];
381 pve_c_trans[i] = (pve_c_trans[i] - double(NSAMPLES) * pme_c_tran * pme_c_tran) / double(NSAMPLES - 1);
382 if (pve_c_trans[i] < 0.0) pve_c_trans[i] = 0.0;
386 int main(int argc, char* argv[]) {
388 printf("Usage: rsem-run-gibbs reference_name imdName statName BURNIN NSAMPLES GAP [-p #Threads] [--seed seed] [-q]\n");
392 strcpy(refName, argv[1]);
393 strcpy(imdName, argv[2]);
394 strcpy(statName, argv[3]);
396 BURNIN = atoi(argv[4]);
397 NSAMPLES = atoi(argv[5]);
404 for (int i = 7; i < argc; i++) {
405 if (!strcmp(argv[i], "-p")) nThreads = atoi(argv[i + 1]);
406 if (!strcmp(argv[i], "--seed")) {
408 int len = strlen(argv[i + 1]);
410 for (int k = 0; k < len; k++) seed = seed * 10 + (argv[i + 1][k] - '0');
412 if (!strcmp(argv[i], "-q")) quiet = true;
416 assert(NSAMPLES > 1); // Otherwise, we cannot calculate posterior variance
418 if (nThreads > NSAMPLES) {
420 printf("Warning: Number of samples is less than number of threads! Change the number of threads to %d!\n", nThreads);
423 load_data(refName, statName, imdName);
424 load_group_info(refName);
425 load_omit_info(imdName);
427 sprintf(modelF, "%s.model", statName);
428 FILE *fi = fopen(modelF, "r");
429 general_assert(fi != NULL, "Cannot open " + cstrtos(modelF) + "!");
430 assert(fscanf(fi, "%d", &model_type) == 1);
433 mw = new double[M + 1]; // make an extra copy
436 case 0 : init_model_related<SingleModel>(modelF); break;
437 case 1 : init_model_related<SingleQModel>(modelF); break;
438 case 2 : init_model_related<PairedEndModel>(modelF); break;
439 case 3 : init_model_related<PairedEndQModel>(modelF); break;
442 if (verbose) printf("Gibbs started!\n");
445 for (int i = 0; i < nThreads; i++) {
446 rc = pthread_create(&threads[i], &attr, Gibbs, (void*)(¶msArray[i]));
447 pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0)!");
449 for (int i = 0; i < nThreads; i++) {
450 rc = pthread_join(threads[i], NULL);
451 pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0)!");
455 if (verbose) printf("Gibbs finished!\n");
457 writeResultsGibbs(M, m, m_trans, gi, gt, ta, alleleS, imdName, pme_c, pme_fpkm, pme_tpm, pve_c, pve_c_genes, pve_c_trans);
459 delete mw; // delete the copy