X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=EM.cpp;h=cf53fc3881b557631dca0ff6fdc4d0bfdd209b5d;hb=229ba9d68e0a801907631887640ab475d51c560c;hp=322a3063ff6596c9b1cce15d037065c8c2c1af85;hpb=86e650e9577999a7ba00ab454d1f6bf674b0ea70;p=rsem.git diff --git a/EM.cpp b/EM.cpp index 322a306..cf53fc3 100644 --- a/EM.cpp +++ b/EM.cpp @@ -10,6 +10,7 @@ #include #include "utils.h" +#include "sampling.h" #include "Read.h" #include "SingleRead.h" @@ -58,6 +59,7 @@ int nThreads; bool genBamF; // If user wants to generate bam file, true; otherwise, false. +bool bamSampling; // true if sampling from read posterior distribution when bam file is generated bool updateModel, calcExpectedWeights; bool genGibbsOut; // generate file for Gibbs sampler @@ -100,7 +102,7 @@ void init(ReadReader **&readers, HitContainer **&hitvs, doubl indices[i] = new ReadIndex(readFs[i]); } for (int i = 0; i < nThreads; i++) { - readers[i] = new ReadReader(s, readFs); + readers[i] = new ReadReader(s, readFs, refs.hasPolyA(), mparams.seedLen); // allow calculation of calc_lq() function readers[i]->setIndices(indices); } @@ -338,7 +340,7 @@ void writeResults(ModelType& model, double* counts) { fprintf(fo, "%.15g%c", tau[i], (i < M ? '\t' : '\n')); for (int i = 1; i <= M; i++) { const Transcript& transcript = transcripts.getTranscriptAt(i); - fprintf(fo, "%s%c", transcript.getLeft().c_str(), (i < M ? '\t' : '\n')); + fprintf(fo, "%s%c", transcript.getGeneID().c_str(), (i < M ? '\t' : '\n')); } fclose(fo); @@ -394,8 +396,6 @@ void release(ReadReader **readers, HitContainer **hitvs, doub delete[] mhps; } -int tmp_n; - inline bool doesUpdateModel(int ROUND) { // return ROUND <= 20 || ROUND % 100 == 0; return ROUND <= 10; @@ -465,14 +465,18 @@ void EM() { //E step for (int i = 0; i < nThreads; i++) { rc = pthread_create(&threads[i], &attr, E_STEP, (void*)(&fparams[i])); - if (rc != 0) { fprintf(stderr, "Cannot create thread %d at ROUND %d! (numbered from 0)\n", i, ROUND); exit(-1); } - //assert(rc == 0); + if (rc != 0) { + fprintf(stderr, "Cannot create thread %d at ROUND %d! (numbered from 0)", i, ROUND); + pthread_exception(rc); + } } for (int i = 0; i < nThreads; i++) { rc = pthread_join(threads[i], &status); - if (rc != 0) { fprintf(stderr, "Cannot join thread %d at ROUND %d! (numbered from 0)\n", i, ROUND); exit(-1); } - //assert(rc == 0); + if (rc != 0) { + fprintf(stderr, "Cannot join thread %d at ROUND %d! (numbered from 0)\n", i, ROUND); + pthread_exception(rc); + } } model.setNeedCalcConPrb(false); @@ -518,11 +522,17 @@ void EM() { if (model.getNeedCalcConPrb()) { for (int i = 0; i < nThreads; i++) { rc = pthread_create(&threads[i], &attr, calcConProbs, (void*)(&fparams[i])); - if (rc != 0) { fprintf(stderr, "Cannot create thread %d when generate files for Gibbs sampler! (numbered from 0)\n", i); exit(-1); } + if (rc != 0) { + fprintf(stderr, "Cannot create thread %d when generate files for Gibbs sampler! (numbered from 0)\n", i); + pthread_exception(rc); + } } for (int i = 0; i < nThreads; i++) { rc = pthread_join(threads[i], &status); - if (rc != 0) { fprintf(stderr, "Cannot join thread %d when generate files for Gibbs sampler! (numbered from 0)\n", i); exit(-1); } + if (rc != 0) { + fprintf(stderr, "Cannot join thread %d when generate files for Gibbs sampler! (numbered from 0)\n", i); + pthread_exception(rc); + } } } model.setNeedCalcConPrb(false); @@ -576,13 +586,17 @@ void EM() { for (int i = 0; i <= M; i++) probv[i] = theta[i]; for (int i = 0; i < nThreads; i++) { rc = pthread_create(&threads[i], &attr, E_STEP, (void*)(&fparams[i])); - if (rc != 0) { fprintf(stderr, "Cannot create thread %d when calculate expected weights! (numbered from 0)\n", i); exit(-1); } - //assert(rc == 0); + if (rc != 0) { + fprintf(stderr, "Cannot create thread %d when calculate expected weights! (numbered from 0)\n", i); + pthread_exception(rc); + } } for (int i = 0; i < nThreads; i++) { rc = pthread_join(threads[i], &status); - if (rc != 0) { fprintf(stderr, "Cannot join thread %d! (numbered from 0) when calculate expected weights!\n", i); exit(-1); } - //assert(rc == 0); + if (rc != 0) { + fprintf(stderr, "Cannot join thread %d! (numbered from 0) when calculate expected weights!\n", i); + pthread_exception(rc); + } } model.setNeedCalcConPrb(false); for (int i = 1; i < nThreads; i++) { @@ -596,7 +610,7 @@ void EM() { pthread_attr_destroy(&attr); //convert theta' to theta - double *mw = model.getMW(); + double *mw = model.getMW(); sum = 0.0; for (int i = 0; i <= M; i++) { theta[i] = (mw[i] < EPSILON ? 0.0 : theta[i] / mw[i]); @@ -614,15 +628,34 @@ void EM() { writeResults(model, countvs[0]); if (genBamF) { - sprintf(outBamF, "%s.bam", outName); - if (transcripts.getType() == 0) { - sprintf(chr_list, "%s.chrlist", refName); - pt_chr_list = (char*)(&chr_list); + sprintf(outBamF, "%s.transcript.bam", outName); + + if (bamSampling) { + int local_N; + int fr, to, len, id; + vector arr; + arr.clear(); + + if (verbose) printf("Begin to sample reads from their posteriors.\n"); + for (int i = 0; i < nThreads; i++) { + local_N = hitvs[i]->getN(); + for (int j = 0; j < local_N; j++) { + fr = hitvs[i]->getSAt(j); + to = hitvs[i]->getSAt(j + 1); + len = to - fr + 1; + arr.resize(len); + arr[0] = ncpvs[i][j]; + for (int k = fr; k < to; k++) arr[k - fr + 1] = arr[k - fr] + hitvs[i]->getHitAt(k).getConPrb(); + id = (arr[len - 1] < EPSILON ? -1 : sample(arr, len)); // if all entries in arr are 0, let id be -1 + for (int k = fr; k < to; k++) hitvs[i]->getHitAt(k).setConPrb(k - fr + 1 == id ? 1.0 : 0.0); + } + } + if (verbose) printf("Sampling is finished.\n"); } - BamWriter writer(inpSamType, inpSamF, pt_fn_list, outBamF, pt_chr_list); + BamWriter writer(inpSamType, inpSamF, pt_fn_list, outBamF, transcripts); HitWrapper wrapper(nThreads, hitvs); - writer.work(wrapper, transcripts); + writer.work(wrapper); } release(readers, hitvs, ncpvs, mhps); @@ -633,7 +666,7 @@ int main(int argc, char* argv[]) { bool quiet = false; if (argc < 5) { - printf("Usage : rsem-run-em refName read_type sampleName sampleToken [-p #Threads] [-b samInpType samInpF has_fn_list_? [fn_list]] [-q] [--gibbs-out]\n\n"); + printf("Usage : rsem-run-em refName read_type sampleName sampleToken [-p #Threads] [-b samInpType samInpF has_fn_list_? [fn_list]] [-q] [--gibbs-out] [--sampling]\n\n"); printf(" refName: reference name\n"); printf(" read_type: 0 single read without quality score; 1 single read with quality score; 2 paired-end read without quality score; 3 paired-end read with quality score.\n"); printf(" sampleName: sample's name, including the path\n"); @@ -641,7 +674,8 @@ int main(int argc, char* argv[]) { printf(" -p: number of threads which user wants to use. (default: 1)\n"); printf(" -b: produce bam format output file. (default: off)\n"); printf(" -q: set it quiet\n"); - printf(" --gibbs-out: generate output file use by Gibbs sampler. (default: off)\n"); + printf(" --gibbs-out: generate output file used by Gibbs sampler. (default: off)\n"); + printf(" --sampling: sample each read from its posterior distribution when bam file is generated. (default: off)\n"); printf("// model parameters should be in imdName.mparams.\n"); exit(-1); } @@ -657,6 +691,7 @@ int main(int argc, char* argv[]) { nThreads = 1; genBamF = false; + bamSampling = false; genGibbsOut = false; pt_fn_list = pt_chr_list = NULL; @@ -673,6 +708,7 @@ int main(int argc, char* argv[]) { } if (!strcmp(argv[i], "-q")) { quiet = true; } if (!strcmp(argv[i], "--gibbs-out")) { genGibbsOut = true; } + if (!strcmp(argv[i], "--sampling")) { bamSampling = true; } } if (nThreads <= 0) { fprintf(stderr, "Number of threads should be bigger than 0!\n"); exit(-1); } //assert(nThreads > 0);