X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=Gibbs.cpp;h=e7a1182ce4b8ca78e5a41f59a0cdf5ddf657020a;hb=810375ef69e0de93188c6fa3819cfc584a32ff34;hp=63d35a829c059e12ae79bd8f3246fcec5b50438b;hpb=509ffe2d71cf6a6f5cca8c39909f8ee10b8db899;p=rsem.git diff --git a/Gibbs.cpp b/Gibbs.cpp index 63d35a8..e7a1182 100644 --- a/Gibbs.cpp +++ b/Gibbs.cpp @@ -18,7 +18,9 @@ #include "PairedEndQModel.h" #include "Refs.h" + #include "GroupInfo.h" +#include "WriteResults.h" using namespace std; @@ -27,7 +29,7 @@ struct Params { FILE *fo; engine_type *engine; double *pme_c, *pve_c; //posterior mean and variance vectors on counts - double *pme_theta; + double *pme_tpm, *pme_fpkm; }; @@ -44,25 +46,25 @@ struct Item { int nThreads; int model_type; -int m, M; +int M; READ_INT_TYPE N0, N1; HIT_INT_TYPE nHits; double totc; int BURNIN, NSAMPLES, GAP; -char imdName[STRLEN], statName[STRLEN]; -char thetaF[STRLEN], ofgF[STRLEN], groupF[STRLEN], refF[STRLEN], modelF[STRLEN]; +char refName[STRLEN], imdName[STRLEN], statName[STRLEN]; +char thetaF[STRLEN], ofgF[STRLEN], refF[STRLEN], modelF[STRLEN]; char cvsF[STRLEN]; Refs refs; -GroupInfo gi; vector s; vector hits; -vector theta; +vector eel; +double *mw; vector pme_c, pve_c; //global posterior mean and variance vectors on counts -vector pme_theta, eel; +vector pme_tpm, pme_fpkm; bool var_opt; bool quiet; @@ -70,34 +72,18 @@ bool quiet; Params *paramsArray; pthread_t *threads; pthread_attr_t attr; -void *status; int rc; -void load_data(char* reference_name, char* statName, char* imdName) { +void load_data(char* refName, char* statName, char* imdName) { ifstream fin; string line; int tmpVal; //load reference file - sprintf(refF, "%s.seq", reference_name); + sprintf(refF, "%s.seq", refName); refs.loadRefs(refF, 1); M = refs.getM(); - //load groupF - sprintf(groupF, "%s.grp", reference_name); - gi.load(groupF); - m = gi.getm(); - - //load thetaF - sprintf(thetaF, "%s.theta",statName); - fin.open(thetaF); - general_assert(fin.is_open(), "Cannot open " + cstrtos(thetaF) + "!"); - fin>>tmpVal; - general_assert(tmpVal == M + 1, "Number of transcripts is not consistent in " + cstrtos(refF) + " and " + cstrtos(thetaF) + "!"); - theta.assign(M + 1, 0); - for (int i = 0; i <= M; i++) fin>>theta[i]; - fin.close(); - //load ofgF; sprintf(ofgF, "%s.ofg", imdName); fin.open(ofgF); @@ -128,6 +114,15 @@ void load_data(char* reference_name, char* statName, char* imdName) { if (verbose) { printf("Loading Data is finished!\n"); } } +template +void init_model_related(char* modelF) { + ModelType model; + model.read(modelF); + + calcExpectedEffectiveLengths(M, refs, model, eel); + memcpy(mw, model.getMW(), sizeof(double) * (M + 1)); // otherwise, after exiting this procedure, mw becomes undefined +} + // assign threads void init() { int quotient, left; @@ -154,8 +149,10 @@ void init() { memset(paramsArray[i].pme_c, 0, sizeof(double) * (M + 1)); paramsArray[i].pve_c = new double[M + 1]; memset(paramsArray[i].pve_c, 0, sizeof(double) * (M + 1)); - paramsArray[i].pme_theta = new double[M + 1]; - memset(paramsArray[i].pme_theta, 0, sizeof(double) * (M + 1)); + paramsArray[i].pme_tpm = new double[M + 1]; + memset(paramsArray[i].pme_tpm, 0, sizeof(double) * (M + 1)); + paramsArray[i].pme_fpkm = new double[M + 1]; + memset(paramsArray[i].pme_fpkm, 0, sizeof(double) * (M + 1)); } /* set thread attribute to be joinable */ @@ -193,7 +190,7 @@ void* Gibbs(void* arg) { HIT_INT_TYPE len, fr, to; Params *params = (Params*)arg; - vector theta; + vector theta, tpm, fpkm; vector z, counts; vector arr; @@ -238,10 +235,14 @@ void* Gibbs(void* arg) { if (ROUND > BURNIN) { if ((ROUND - BURNIN - 1) % GAP == 0) { writeCountVector(params->fo, counts); + for (int i = 0; i <= M; i++) theta[i] = counts[i] / totc; + polishTheta(M, theta, eel, mw); + calcExpressionValues(M, theta, eel, tpm, fpkm); for (int i = 0; i <= M; i++) { params->pme_c[i] += counts[i] - 1; params->pve_c[i] += (counts[i] - 1) * (counts[i] - 1); - params->pme_theta[i] += counts[i] / totc; + params->pme_tpm[i] += tpm[i]; + params->pme_fpkm[i] += fpkm[i]; } } } @@ -262,18 +263,21 @@ void release() { pme_c.assign(M + 1, 0); pve_c.assign(M + 1, 0); - pme_theta.assign(M + 1, 0); + pme_tpm.assign(M + 1, 0); + pme_fpkm.assign(M + 1, 0); for (int i = 0; i < nThreads; i++) { fclose(paramsArray[i].fo); delete paramsArray[i].engine; for (int j = 0; j <= M; j++) { pme_c[j] += paramsArray[i].pme_c[j]; pve_c[j] += paramsArray[i].pve_c[j]; - pme_theta[j] += paramsArray[i].pme_theta[j]; + pme_tpm[j] += paramsArray[i].pme_tpm[j]; + pme_fpkm[j] += paramsArray[i].pme_fpkm[j]; } delete[] paramsArray[i].pme_c; delete[] paramsArray[i].pve_c; - delete[] paramsArray[i].pme_theta; + delete[] paramsArray[i].pme_tpm; + delete[] paramsArray[i].pme_fpkm; } delete[] paramsArray; @@ -281,133 +285,18 @@ void release() { for (int i = 0; i <= M; i++) { pme_c[i] /= NSAMPLES; pve_c[i] = (pve_c[i] - NSAMPLES * pme_c[i] * pme_c[i]) / (NSAMPLES - 1); - pme_theta[i] /= NSAMPLES; + pme_tpm[i] /= NSAMPLES; + pme_fpkm[i] /= NSAMPLES; } } -template -void calcExpectedEffectiveLengths(ModelType& model) { - int lb, ub, span; - double *pdf = NULL, *cdf = NULL, *clen = NULL; // clen[i] = \sigma_{j=1}^{i}pdf[i]*(lb+i) - - model.getGLD().copyTo(pdf, cdf, lb, ub, span); - clen = new double[span + 1]; - clen[0] = 0.0; - for (int i = 1; i <= span; i++) { - clen[i] = clen[i - 1] + pdf[i] * (lb + i); - } - - eel.assign(M + 1, 0.0); - for (int i = 1; i <= M; i++) { - int totLen = refs.getRef(i).getTotLen(); - int fullLen = refs.getRef(i).getFullLen(); - int pos1 = max(min(totLen - fullLen + 1, ub) - lb, 0); - int pos2 = max(min(totLen, ub) - lb, 0); - - if (pos2 == 0) { eel[i] = 0.0; continue; } - - eel[i] = fullLen * cdf[pos1] + ((cdf[pos2] - cdf[pos1]) * (totLen + 1) - (clen[pos2] - clen[pos1])); - assert(eel[i] >= 0); - if (eel[i] < MINEEL) { eel[i] = 0.0; } - } - - delete[] pdf; - delete[] cdf; - delete[] clen; -} - -template -void writeEstimatedParameters(char* modelF, char* imdName) { - ModelType model; - double denom; - char outF[STRLEN]; - FILE *fo; - - model.read(modelF); - - calcExpectedEffectiveLengths(model); - - denom = pme_theta[0]; - for (int i = 1; i <= M; i++) - if (eel[i] < EPSILON) pme_theta[i] = 0.0; - else denom += pme_theta[i]; - - general_assert(denom >= EPSILON, "No Expected Effective Length is no less than " + ftos(MINEEL, 6) + "?!"); - - for (int i = 0; i <= M; i++) pme_theta[i] /= denom; - - denom = 0.0; - double *mw = model.getMW(); - for (int i = 0; i <= M; i++) { - pme_theta[i] = (mw[i] < EPSILON ? 0.0 : pme_theta[i] / mw[i]); - denom += pme_theta[i]; - } - assert(denom >= EPSILON); - for (int i = 0; i <= M; i++) pme_theta[i] /= denom; - - //calculate tau values - double *tau = new double[M + 1]; - memset(tau, 0, sizeof(double) * (M + 1)); - - denom = 0.0; - for (int i = 1; i <= M; i++) - if (eel[i] > EPSILON) { - tau[i] = pme_theta[i] / eel[i]; - denom += tau[i]; - } - - general_assert(denom >= EPSILON, "No alignable reads?!"); - - for (int i = 1; i <= M; i++) { - tau[i] /= denom; - } - - //isoform level results - sprintf(outF, "%s.iso_res", imdName); - fo = fopen(outF, "a"); - general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!"); - - for (int i = 1; i <= M; i++) - fprintf(fo, "%.2f%c", pme_c[i], (i < M ? '\t' : '\n')); - for (int i = 1; i <= M; i++) - fprintf(fo, "%.15g%c", tau[i], (i < M ? '\t' : '\n')); - - fclose(fo); - - //gene level results - sprintf(outF, "%s.gene_res", imdName); - fo = fopen(outF, "a"); - general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!"); - - for (int i = 0; i < m; i++) { - double sumC = 0.0; // sum of pme counts - int b = gi.spAt(i), e = gi.spAt(i + 1); - for (int j = b; j < e; j++) { - sumC += pme_c[j]; - } - fprintf(fo, "%.15g%c", sumC, (i < m - 1 ? '\t' : '\n')); - } - for (int i = 0; i < m; i++) { - double sumT = 0.0; // sum of tau values - int b = gi.spAt(i), e = gi.spAt(i + 1); - for (int j = b; j < e; j++) { - sumT += tau[j]; - } - fprintf(fo, "%.15g%c", sumT, (i < m - 1 ? '\t' : '\n')); - } - fclose(fo); - - delete[] tau; - - if (verbose) { printf("Gibbs based expression values are written!\n"); } -} - int main(int argc, char* argv[]) { if (argc < 7) { printf("Usage: rsem-run-gibbs reference_name imdName statName BURNIN NSAMPLES GAP [-p #Threads] [--var] [-q]\n"); exit(-1); } + strcpy(refName, argv[1]); strcpy(imdName, argv[2]); strcpy(statName, argv[3]); @@ -415,8 +304,6 @@ int main(int argc, char* argv[]) { NSAMPLES = atoi(argv[5]); GAP = atoi(argv[6]); - load_data(argv[1], statName, imdName); - nThreads = 1; var_opt = false; quiet = false; @@ -435,6 +322,23 @@ int main(int argc, char* argv[]) { printf("Warning: Number of samples is less than number of threads! Change the number of threads to %d!\n", nThreads); } + load_data(refName, statName, imdName); + + sprintf(modelF, "%s.model", statName); + FILE *fi = fopen(modelF, "r"); + general_assert(fi != NULL, "Cannot open " + cstrtos(modelF) + "!"); + assert(fscanf(fi, "%d", &model_type) == 1); + fclose(fi); + + mw = new double[M + 1]; // make an extra copy + + switch(model_type) { + case 0 : init_model_related(modelF); break; + case 1 : init_model_related(modelF); break; + case 2 : init_model_related(modelF); break; + case 3 : init_model_related(modelF); break; + } + if (verbose) printf("Gibbs started!\n"); init(); @@ -443,29 +347,26 @@ int main(int argc, char* argv[]) { pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0)!"); } for (int i = 0; i < nThreads; i++) { - rc = pthread_join(threads[i], &status); + rc = pthread_join(threads[i], NULL); pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0)!"); } release(); if (verbose) printf("Gibbs finished!\n"); - - sprintf(modelF, "%s.model", statName); - FILE *fi = fopen(modelF, "r"); - general_assert(fi != NULL, "Cannot open " + cstrtos(modelF) + "!"); - assert(fscanf(fi, "%d", &model_type) == 1); - fclose(fi); - - switch(model_type) { - case 0 : writeEstimatedParameters(modelF, imdName); break; - case 1 : writeEstimatedParameters(modelF, imdName); break; - case 2 : writeEstimatedParameters(modelF, imdName); break; - case 3 : writeEstimatedParameters(modelF, imdName); break; - } + + writeResultsGibbs(M, refName, imdName, pme_c, pme_fpkm, pme_tpm); if (var_opt) { char varF[STRLEN]; + // Load group info + int m; + GroupInfo gi; + char groupF[STRLEN]; + sprintf(groupF, "%s.grp", refName); + gi.load(groupF); + m = gi.getm(); + sprintf(varF, "%s.var", statName); FILE *fo = fopen(varF, "w"); general_assert(fo != NULL, "Cannot open " + cstrtos(varF) + "!"); @@ -477,6 +378,8 @@ int main(int argc, char* argv[]) { } fclose(fo); } + + delete mw; // delete the copy return 0; }