X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=calcCI.cpp;h=c9ccce5909a2be59baabe8b3345232256e8a2da4;hb=53e8bbb15e0bfed6a0caae7b5ba6777a9a942266;hp=581e43470dda10032600caa319bc269c700f5d1f;hpb=a95154919f950f86de9104b2b9dcf1f0c7e83387;p=rsem.git diff --git a/calcCI.cpp b/calcCI.cpp index 581e434..c9ccce5 100644 --- a/calcCI.cpp +++ b/calcCI.cpp @@ -50,7 +50,7 @@ ofstream ftmpOut; int *cvec; double *theta; -CIType *iso_nrf, *gene_nrf, *iso_tau, *gene_tau; +CIType *iso_tau, *gene_tau; engine_type engine(time(NULL)); distribution_type **gammas; @@ -59,6 +59,7 @@ generator_type **rgs; int M, m; Refs refs; GroupInfo gi; +char imdName[STRLEN], statName[STRLEN]; char modelF[STRLEN], groupF[STRLEN], refF[STRLEN]; vector eel; //expected effective lengths @@ -263,45 +264,37 @@ void calcCI(int nSamples, float *samples, float &lb, float &ub) { } void generateResults(char* imdName) { - float *izsamples, *gzsamples, *itsamples, *gtsamples; + float *itsamples, *gtsamples; ifstream fin; FILE *fo; char outF[STRLEN]; - iso_nrf = new CIType[M + 1]; iso_tau = new CIType[M + 1]; - gene_nrf = new CIType[m]; gene_tau = new CIType[m]; - izsamples = new float[nSamples]; itsamples = new float[nSamples]; - gzsamples = new float[nSamples]; gtsamples = new float[nSamples]; fin.open(tmpF, ios::binary); - for (int k = 0; k < nSamples; k++) fin.read((char*)(&izsamples[k]), FLOATSIZE); - calcCI(nSamples, izsamples, iso_nrf[0].lb, iso_nrf[0].ub); + //read sampled theta0 values + for (int k = 0; k < nSamples; k++) fin.read((char*)(&itsamples[k]), FLOATSIZE); for (int i = 0; i < m; i++) { int b = gi.spAt(i), e = gi.spAt(i + 1); - memset(gzsamples, 0, FLOATSIZE * nSamples); memset(gtsamples, 0, FLOATSIZE * nSamples); for (int j = b; j < e; j++) { for (int k = 0; k < nSamples; k++) { - fin.read((char*)(&izsamples[k]), FLOATSIZE); - if (eel[j] > EPSILON && tau_denoms[k] > EPSILON) { itsamples[k] = izsamples[k] / eel[j] / tau_denoms[k]; } + fin.read((char*)(&itsamples[k]), FLOATSIZE); + if (eel[j] > EPSILON && tau_denoms[k] > EPSILON) { itsamples[k] = itsamples[k] / eel[j] / tau_denoms[k]; } else { - if (izsamples[k] != 0.0) { fprintf(stderr, "K=%d, IZSAMPLES_K=%lf\n", k, izsamples[k]); exit(-1); } + if (itsamples[k] != 0.0) { fprintf(stderr, "Not equal to 0! K=%d, Sampled Theta Value=%lf\n", k, itsamples[k]); exit(-1); } itsamples[k] = 0.0; } - gzsamples[k] += izsamples[k]; gtsamples[k] += itsamples[k]; } - calcCI(nSamples, izsamples, iso_nrf[j].lb, iso_nrf[j].ub); calcCI(nSamples, itsamples, iso_tau[j].lb, iso_tau[j].ub); } - calcCI(nSamples, gzsamples, gene_nrf[i].lb, gene_nrf[i].ub); calcCI(nSamples, gtsamples, gene_tau[i].lb, gene_tau[i].ub); if (verbose && (i + 1) % 1000 == 0) { printf("%d genes are done!\n", i + 1); } @@ -312,10 +305,6 @@ void generateResults(char* imdName) { //isoform level results sprintf(outF, "%s.iso_res", imdName); fo = fopen(outF, "a"); - for (int i = 1; i <= M; i++) - fprintf(fo, "%.6g%c", iso_nrf[i].lb, (i < M ? '\t' : '\n')); - for (int i = 1; i <= M; i++) - fprintf(fo, "%.6g%c", iso_nrf[i].ub, (i < M ? '\t' : '\n')); for (int i = 1; i <= M; i++) fprintf(fo, "%.6g%c", iso_tau[i].lb, (i < M ? '\t' : '\n')); for (int i = 1; i <= M; i++) @@ -325,42 +314,31 @@ void generateResults(char* imdName) { //gene level results sprintf(outF, "%s.gene_res", imdName); fo = fopen(outF, "a"); - for (int i = 0; i < m; i++) - fprintf(fo, "%.6g%c", gene_nrf[i].lb, (i < m - 1 ? '\t' : '\n')); - for (int i = 0; i < m; i++) - fprintf(fo, "%.6g%c", gene_nrf[i].ub, (i < m - 1 ? '\t' : '\n')); for (int i = 0; i < m; i++) fprintf(fo, "%.6g%c", gene_tau[i].lb, (i < m - 1 ? '\t' : '\n')); for (int i = 0; i < m; i++) fprintf(fo, "%.6g%c", gene_tau[i].ub, (i < m - 1 ? '\t' : '\n')); fclose(fo); - printf("CI of noise isoform is [%.6g, %.6g]\n", iso_nrf[0].lb, iso_nrf[0].ub); - - delete[] izsamples; delete[] itsamples; - delete[] gzsamples; delete[] gtsamples; - delete[] iso_nrf; delete[] iso_tau; - delete[] gene_nrf; delete[] gene_tau; if (verbose) { printf("All credibility intervals are calculated!\n"); } - sprintf(outF, "%s.tau_denoms", imdName); - fo = fopen(outF, "w"); - fprintf(fo, "%d\n", nSamples); - for (int i = 0; i < nSamples; i++) fprintf(fo, "%.15g ", tau_denoms[i]); - fprintf(fo, "\n"); - fclose(fo); - + sprintf(outF, "%s.tau_denoms", imdName); + fo = fopen(outF, "w"); + fprintf(fo, "%d\n", nSamples); + for (int i = 0; i < nSamples; i++) fprintf(fo, "%.15g ", tau_denoms[i]); + fprintf(fo, "\n"); + fclose(fo); } int main(int argc, char* argv[]) { if (argc < 7) { - printf("Usage: rsem-calculate-credibility-intervals reference_name sample_name imdName confidence nSpC nMB[-q]\n"); + printf("Usage: rsem-calculate-credibility-intervals reference_name sample_name sampleToken confidence nSpC nMB[-q]\n"); exit(-1); } @@ -374,7 +352,10 @@ int main(int argc, char* argv[]) { } verbose = !quiet; - sprintf(modelF, "%s.model", argv[2]); + sprintf(imdName, "%s.temp/%s", argv[2], argv[3]); + sprintf(statName, "%s.stat/%s", argv[2], argv[3]); + + sprintf(modelF, "%s.model", statName); FILE *fi = fopen(modelF, "r"); if (fi == NULL) { fprintf(stderr, "Cannot open %s!\n", modelF); exit(-1); } fscanf(fi, "%d", &model_type); @@ -387,8 +368,8 @@ int main(int argc, char* argv[]) { gi.load(groupF); m = gi.getm(); - sprintf(tmpF, "%s.tmp", argv[3]); - sprintf(cvsF, "%s.countvectors", argv[3]); + sprintf(tmpF, "%s.tmp", imdName); + sprintf(cvsF, "%s.countvectors", imdName); switch(model_type) { case 0 : sampling(); break; @@ -397,15 +378,15 @@ int main(int argc, char* argv[]) { case 3 : sampling(); break; } - generateResults(argv[3]); + generateResults(imdName); delete[] tau_denoms; sprintf(command, "rm -f %s", tmpF); int status = system(command); if (status != 0) { - fprintf(stderr, "Cannot delete %s!\n", tmpF); - exit(-1); + fprintf(stderr, "Cannot delete %s!\n", tmpF); + exit(-1); } return 0;