X-Git-Url: https://git.donarmstrong.com/?p=rsem.git;a=blobdiff_plain;f=Gibbs.cpp;h=7563dd49f508129992d6a52aa70ed3a8c3130b40;hp=bf8563e16ea673b467e22b7d51343b7bd63b27e3;hb=HEAD;hpb=9c2e46183a19d661f0a618a8eabe8ce1f6a8e2d6 diff --git a/Gibbs.cpp b/Gibbs.cpp index bf8563e..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 + 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,17 +47,16 @@ 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; @@ -62,10 +64,10 @@ vector hits; vector eel; double *mw; +vector pseudo_counts; vector pme_c, pve_c; //global posterior mean and variance vectors on counts vector pme_tpm, pme_fpkm; -bool var_opt; bool quiet; Params *paramsArray; @@ -73,21 +75,28 @@ pthread_t *threads; pthread_attr_t attr; 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 ofgF; sprintf(ofgF, "%s.ofg", imdName); fin.open(ofgF); @@ -115,38 +124,30 @@ 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"); } } -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) +void load_group_info(char* refName) { + // Load group info + sprintf(groupF, "%s.grp", refName); + gi.load(groupF); + m = gi.getm(); - 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); - } + alleleS = isAlleleSpecific(refName, >, &ta); // if allele-specific + m_trans = (alleleS ? ta.getm() : 0); - 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; + 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 @@ -154,7 +155,7 @@ void init_model_related(char* modelF) { ModelType model; model.read(modelF); - calcExpectedEffectiveLengths(model); + calcExpectedEffectiveLengths(M, refs, model, eel); memcpy(mw, model.getMW(), sizeof(double) * (M + 1)); // otherwise, after exiting this procedure, mw becomes undefined } @@ -170,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; @@ -188,7 +190,17 @@ void init() { 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); @@ -206,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); @@ -220,72 +232,21 @@ void writeCountVector(FILE* fo, vector& counts) { fprintf(fo, "%d\n", counts[M]); } -void polishTheta(vector& theta, const vector& eel, const double* mw) { - double sum = 0.0; - - /* The reason that for noise gene, mw value is 1 is : - * currently, all masked positions are for poly(A) sites, which in theory should be filtered out. - * So the theta0 does not containing reads from any masked position - */ - - for (int i = 0; i <= M; i++) { - // i == 0, mw[i] == 1 - if (i > 0 && (mw[i] < EPSILON || eel[i] < EPSILON)) { - theta[i] = 0.0; - continue; - } - theta[i] = theta[i] / mw[i]; - sum += theta[i]; - } - // currently is OK, since no transcript should be masked totally, only the poly(A) tail related part will be masked - general_assert(sum >= EPSILON, "No effective length is no less than" + ftos(MINEEL, 6) + " !"); - for (int i = 0; i <= M; i++) theta[i] /= sum; -} - -void calcExpressionValues(const vector& theta, const vector& eel, vector& tpm, vector& fpkm) { - double denom; - vector frac; - - //calculate fraction of count over all mappabile reads - denom = 0.0; - frac.assign(M + 1, 0.0); - for (int i = 1; i <= M; i++) - if (eel[i] >= EPSILON) { - frac[i] = theta[i]; - denom += frac[i]; - } - general_assert(denom > 0, "No alignable reads?!"); - for (int i = 1; i <= M; i++) frac[i] /= denom; - - //calculate FPKM - fpkm.assign(M + 1, 0.0); - for (int i = 1; i <= M; i++) - if (eel[i] >= EPSILON) fpkm[i] = frac[i] * 1e9 / eel[i]; - - //calculate TPM - tpm.assign(M + 1, 0.0); - denom = 0.0; - for (int i = 1; i <= M; i++) denom += fpkm[i]; - for (int i = 1; i <= M; i++) tpm[i] = fpkm[i] / denom * 1e6; -} - void* Gibbs(void* arg) { int CHAINLEN; HIT_INT_TYPE len, fr, to; Params *params = (Params*)arg; vector theta, tpm, fpkm; - vector z, counts; + 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++) { @@ -320,14 +281,29 @@ void* Gibbs(void* arg) { if ((ROUND - BURNIN - 1) % GAP == 0) { writeCountVector(params->fo, counts); for (int i = 0; i <= M; i++) theta[i] = counts[i] / totc; - polishTheta(theta, eel, mw); - calcExpressionValues(theta, eel, tpm, fpkm); + 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_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; + } } } @@ -349,6 +325,11 @@ void release() { pve_c.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; @@ -358,82 +339,57 @@ void release() { 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_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); + 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; } -} - -void writeResults(char* imdName) { - char outF[STRLEN]; - FILE *fo; - - vector isopct; - vector gene_counts, gene_tpm, gene_fpkm; - - //calculate IsoPct, etc. - isopct.assign(M + 1, 0.0); - gene_counts.assign(m, 0.0); gene_tpm.assign(m, 0.0); gene_fpkm.assign(m, 0.0); for (int i = 0; i < m; i++) { - int b = gi.spAt(i), e = gi.spAt(i + 1); - for (int j = b; j < e; j++) { - gene_counts[i] += pme_c[j]; - gene_tpm[i] += pme_tpm[j]; - gene_fpkm[i] += pme_fpkm[j]; - } - if (gene_tpm[i] < EPSILON) continue; - for (int j = b; j < e; j++) - isopct[j] = pme_tpm[j] / gene_tpm[i]; + 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; } - //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, "%.2f%c", pme_tpm[i], (i < M ? '\t' : '\n')); - for (int i = 1; i <= M; i++) - fprintf(fo, "%.2f%c", pme_fpkm[i], (i < M ? '\t' : '\n')); - for (int i = 1; i <= M; i++) - fprintf(fo, "%.2f%c", isopct[i] * 1e2, (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++) - fprintf(fo, "%.2f%c", gene_counts[i], (i < m - 1 ? '\t' : '\n')); - for (int i = 0; i < m; i++) - fprintf(fo, "%.2f%c", gene_tpm[i], (i < m - 1 ? '\t' : '\n')); - for (int i = 0; i < m; i++) - fprintf(fo, "%.2f%c", gene_fpkm[i], (i < m - 1 ? '\t' : '\n')); - fclose(fo); - - 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]); @@ -442,12 +398,17 @@ int main(int argc, char* argv[]) { GAP = atoi(argv[6]); 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; @@ -459,7 +420,9 @@ 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(argv[1], statName, imdName); + load_data(refName, statName, imdName); + load_group_info(refName); + load_omit_info(imdName); sprintf(modelF, "%s.model", statName); FILE *fi = fopen(modelF, "r"); @@ -491,22 +454,7 @@ int main(int argc, char* argv[]) { if (verbose) printf("Gibbs finished!\n"); - writeResults(imdName); - - 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); - } + writeResultsGibbs(M, m, m_trans, gi, gt, ta, alleleS, imdName, pme_c, pme_fpkm, pme_tpm, pve_c, pve_c_genes, pve_c_trans); delete mw; // delete the copy