X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=EM.cpp;h=f9449a85cd3226a3f29356f10fa018e4ca908cba;hb=4a435fbee229af1dd5f76b4feb0ca71398ac8796;hp=ff257680ead4f4f2087515131623a6e61fb498c1;hpb=58f823f5be6dfbe00896fc44ac3aac5e881e9c5c;p=rsem.git diff --git a/EM.cpp b/EM.cpp index ff25768..f9449a8 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); @@ -395,8 +397,8 @@ void release(ReadReader **readers, HitContainer **hitvs, doub } inline bool doesUpdateModel(int ROUND) { - //return false; // never update, for debugging only - return ROUND <= 20 || ROUND % 100 == 0; + // return ROUND <= 20 || ROUND % 100 == 0; + return ROUND <= 10; } //Including initialize, algorithm and results saving @@ -506,7 +508,7 @@ void EM() { } if (verbose) printf("ROUND = %d, SUM = %.15g, bChange = %f, totNum = %d\n", ROUND, sum, bChange, totNum); - } while (ROUND < MIN_ROUND || totNum > 0 && ROUND < MAX_ROUND); + } while (ROUND < MIN_ROUND || (totNum > 0 && ROUND < MAX_ROUND)); //while (ROUND < MAX_ROUND); if (totNum > 0) fprintf(stderr, "Warning: RSEM reaches %d iterations before meeting the convergence criteria.\n", MAX_ROUND); @@ -571,6 +573,7 @@ void EM() { //calculate expected weights and counts using learned parameters updateModel = false; calcExpectedWeights = true; + 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); } @@ -611,15 +614,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); @@ -630,15 +652,16 @@ 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(" samplePath: sample path.\n"); - printf(" sampleName: sample name.\n"); + printf(" sampleName: sample's name, including the path\n"); + printf(" sampleToken: sampleName excludes the path\n"); 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); } @@ -654,6 +677,7 @@ int main(int argc, char* argv[]) { nThreads = 1; genBamF = false; + bamSampling = false; genGibbsOut = false; pt_fn_list = pt_chr_list = NULL; @@ -670,6 +694,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);