X-Git-Url: https://git.donarmstrong.com/?p=rsem.git;a=blobdiff_plain;f=Gibbs.cpp;h=7563dd49f508129992d6a52aa70ed3a8c3130b40;hp=63d35a829c059e12ae79bd8f3246fcec5b50438b;hb=92b24279a3ecc72946e7e7c23149ad0d181f373a;hpb=509ffe2d71cf6a6f5cca8c39909f8ee10b8db899 diff --git a/Gibbs.cpp b/Gibbs.cpp index 63d35a8..7563dd4 100644 --- a/Gibbs.cpp +++ b/Gibbs.cpp @@ -18,19 +18,22 @@ #include "PairedEndQModel.h" #include "Refs.h" + #include "GroupInfo.h" +#include "WriteResults.h" using namespace std; struct Params { - int no, nsamples; - FILE *fo; - engine_type *engine; - double *pme_c, *pve_c; //posterior mean and variance vectors on counts - double *pme_theta; + int no, nsamples; + FILE *fo; + engine_type *engine; + double *pme_c, *pve_c; //posterior mean and variance vectors on counts + double *pme_tpm, *pme_fpkm; + + double *pve_c_genes, *pve_c_trans; }; - struct Item { int sid; double conprb; @@ -44,60 +47,56 @@ 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 pseudo_counts; 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; Params *paramsArray; pthread_t *threads; pthread_attr_t attr; -void *status; int rc; -void load_data(char* reference_name, char* statName, char* imdName) { +bool hasSeed; +seedType seed; + +int m; +char groupF[STRLEN]; +GroupInfo gi; + +bool alleleS; +int m_trans; +GroupInfo gt, ta; +vector pve_c_genes, pve_c_trans; + +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); @@ -125,7 +124,39 @@ void load_data(char* reference_name, char* statName, char* imdName) { totc = N0 + N1 + (M + 1); - if (verbose) { printf("Loading Data is finished!\n"); } + if (verbose) { printf("Loading data is finished!\n"); } +} + +void load_group_info(char* refName) { + // Load group info + sprintf(groupF, "%s.grp", refName); + gi.load(groupF); + m = gi.getm(); + + alleleS = isAlleleSpecific(refName, >, &ta); // if allele-specific + m_trans = (alleleS ? ta.getm() : 0); + + if (verbose) { printf("Loading group information is finished!\n"); } +} + +// Load imdName.omit and initialize the pseudo count vector. +void load_omit_info(const char* imdName) { + char omitF[STRLEN]; + sprintf(omitF, "%s.omit", imdName); + FILE *fi = fopen(omitF, "r"); + pseudo_counts.assign(M + 1, 1); + int tid; + while (fscanf(fi, "%d", &tid) == 1) pseudo_counts[tid] = 0; + fclose(fi); +} + +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 @@ -140,6 +171,7 @@ void init() { paramsArray = new Params[nThreads]; threads = new pthread_t[nThreads]; + hasSeed ? engineFactory::init(seed) : engineFactory::init(); for (int i = 0; i < nThreads; i++) { paramsArray[i].no = i; @@ -154,9 +186,21 @@ 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)); + + paramsArray[i].pve_c_genes = new double[m]; + memset(paramsArray[i].pve_c_genes, 0, sizeof(double) * m); + + paramsArray[i].pve_c_trans = NULL; + if (alleleS) { + paramsArray[i].pve_c_trans = new double[m_trans]; + memset(paramsArray[i].pve_c_trans, 0, sizeof(double) * m_trans); + } } + engineFactory::finish(); /* set thread attribute to be joinable */ pthread_attr_init(&attr); @@ -174,7 +218,7 @@ void sampleTheta(engine_type& engine, vector& theta) { theta.assign(M + 1, 0); denom = 0.0; for (int i = 0; i <= M; i++) { - theta[i] = gmg(); + theta[i] = (pseudo_counts[i] > 0 ? gmg() : 0.0); denom += theta[i]; } assert(denom > EPSILON); @@ -193,18 +237,16 @@ void* Gibbs(void* arg) { HIT_INT_TYPE len, fr, to; Params *params = (Params*)arg; - vector theta; - vector z, counts; + vector theta, tpm, fpkm; + vector z, counts(pseudo_counts); vector arr; - uniform01 rg(*params->engine); + uniform_01_generator rg(*params->engine, uniform_01_dist()); // generate initial state sampleTheta(*params->engine, theta); z.assign(N1, 0); - - counts.assign(M + 1, 1); // 1 pseudo count counts[0] += N0; for (READ_INT_TYPE i = 0; i < N1; i++) { @@ -238,11 +280,30 @@ 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_c[i] += counts[i] - pseudo_counts[i]; + params->pve_c[i] += double(counts[i] - pseudo_counts[i]) * (counts[i] - pseudo_counts[i]); + params->pme_tpm[i] += tpm[i]; + params->pme_fpkm[i] += fpkm[i]; } + + for (int i = 0; i < m; i++) { + int b = gi.spAt(i), e = gi.spAt(i + 1); + double count = 0.0; + for (int j = b; j < e; j++) count += counts[j] - pseudo_counts[j]; + params->pve_c_genes[i] += count * count; + } + + if (alleleS) + for (int i = 0; i < m_trans; i++) { + int b = ta.spAt(i), e = ta.spAt(i + 1); + double count = 0.0; + for (int j = b; j < e; j++) count += counts[j] - pseudo_counts[j]; + params->pve_c_trans[i] += count * count; + } } } @@ -262,152 +323,73 @@ 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); + + pve_c_genes.assign(m, 0); + pve_c_trans.clear(); + if (alleleS) pve_c_trans.assign(m_trans, 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]; } + + for (int j = 0; j < m; j++) + pve_c_genes[j] += paramsArray[i].pve_c_genes[j]; + + if (alleleS) + for (int j = 0; j < m_trans; j++) + pve_c_trans[j] += paramsArray[i].pve_c_trans[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[i].pve_c_genes; + if (alleleS) delete[] paramsArray[i].pve_c_trans; } delete[] paramsArray; - 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; - } -} - -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); + pve_c[i] = (pve_c[i] - double(NSAMPLES) * pme_c[i] * pme_c[i]) / double(NSAMPLES - 1); + if (pve_c[i] < 0.0) pve_c[i] = 0.0; + pme_tpm[i] /= NSAMPLES; + pme_fpkm[i] /= NSAMPLES; } - 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')); + int b = gi.spAt(i), e = gi.spAt(i + 1); + double pme_c_gene = 0.0; + for (int j = b; j < e; j++) pme_c_gene += pme_c[j]; + pve_c_genes[i] = (pve_c_genes[i] - double(NSAMPLES) * pme_c_gene * pme_c_gene) / double(NSAMPLES - 1); + if (pve_c_genes[i] < 0.0) pve_c_genes[i] = 0.0; } - fclose(fo); - delete[] tau; - - if (verbose) { printf("Gibbs based expression values are written!\n"); } + if (alleleS) + for (int i = 0; i < m_trans; i++) { + int b = ta.spAt(i), e = ta.spAt(i + 1); + double pme_c_tran = 0.0; + for (int j = b; j < e; j++) pme_c_tran += pme_c[j]; + pve_c_trans[i] = (pve_c_trans[i] - double(NSAMPLES) * pme_c_tran * pme_c_tran) / double(NSAMPLES - 1); + if (pve_c_trans[i] < 0.0) pve_c_trans[i] = 0.0; + } } 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"); + printf("Usage: rsem-run-gibbs reference_name imdName statName BURNIN NSAMPLES GAP [-p #Threads] [--seed seed] [-q]\n"); exit(-1); } + strcpy(refName, argv[1]); strcpy(imdName, argv[2]); strcpy(statName, argv[3]); @@ -415,15 +397,18 @@ 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; + hasSeed = false; quiet = false; for (int i = 7; i < argc; i++) { if (!strcmp(argv[i], "-p")) nThreads = atoi(argv[i + 1]); - if (!strcmp(argv[i], "--var")) var_opt = true; + if (!strcmp(argv[i], "--seed")) { + hasSeed = true; + int len = strlen(argv[i + 1]); + seed = 0; + for (int k = 0; k < len; k++) seed = seed * 10 + (argv[i + 1][k] - '0'); + } if (!strcmp(argv[i], "-q")) quiet = true; } verbose = !quiet; @@ -435,6 +420,25 @@ 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); + load_group_info(refName); + load_omit_info(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,40 +447,16 @@ 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"); + + writeResultsGibbs(M, m, m_trans, gi, gt, ta, alleleS, imdName, pme_c, pme_fpkm, pme_tpm, pve_c, pve_c_genes, pve_c_trans); - 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; - } - - if (var_opt) { - char varF[STRLEN]; - - sprintf(varF, "%s.var", statName); - FILE *fo = fopen(varF, "w"); - general_assert(fo != NULL, "Cannot open " + cstrtos(varF) + "!"); - for (int i = 0; i < m; i++) { - int b = gi.spAt(i), e = gi.spAt(i + 1), number_of_isoforms = e - b; - for (int j = b; j < e; j++) { - fprintf(fo, "%s\t%d\t%.15g\t%.15g\n", refs.getRef(j).getName().c_str(), number_of_isoforms, pme_c[j], pve_c[j]); - } - } - fclose(fo); - } + delete mw; // delete the copy return 0; }