9 #include "boost/random.hpp"
14 #include "SingleModel.h"
15 #include "SingleQModel.h"
16 #include "PairedEndModel.h"
17 #include "PairedEndQModel.h"
20 #include "GroupInfo.h"
24 typedef unsigned long bufsize_type;
25 typedef boost::mt19937 engine_type;
26 typedef boost::gamma_distribution<> distribution_type;
27 typedef boost::variate_generator<engine_type&, distribution_type> generator_type;
29 const int FLOATSIZE = sizeof(float);
32 float lb, ub; // the interval is [lb, ub]
34 CIType() { lb = ub = 0.0; }
42 int nC, cvlen, nSpC, nSamples; // nSpC : number of sample theta vectors per count vector
43 int fr, to; // each flush, sample fr .. to - 1
48 char cvsF[STRLEN], tmpF[STRLEN], command[STRLEN];
53 CIType *iso_tau, *gene_tau;
55 engine_type engine(time(NULL));
56 distribution_type **gammas;
62 char imdName[STRLEN], statName[STRLEN];
63 char modelF[STRLEN], groupF[STRLEN], refF[STRLEN];
65 vector<double> eel; //expected effective lengths
66 double *tau_denoms; // denominators for tau values
68 template<class ModelType>
69 void calcExpectedEffectiveLengths(ModelType& model) {
71 double *pdf = NULL, *cdf = NULL, *clen = NULL; // clen[i] = sigma_{j=1}^{i}pdf[i]*(lb+i)
73 model.getGLD().copyTo(pdf, cdf, lb, ub, span);
74 clen = new double[span + 1];
76 for (int i = 1; i <= span; i++) {
77 clen[i] = clen[i - 1] + pdf[i] * (lb + i);
81 eel.resize(M + 1, 0.0);
82 for (int i = 1; i <= M; i++) {
83 int totLen = refs.getRef(i).getTotLen();
84 int fullLen = refs.getRef(i).getFullLen();
85 int pos1 = max(min(totLen - fullLen + 1, ub) - lb, 0);
86 int pos2 = max(min(totLen, ub) - lb, 0);
88 if (pos2 == 0) { eel[i] = 0.0; continue; }
90 eel[i] = fullLen * cdf[pos1] + ((cdf[pos2] - cdf[pos1]) * (totLen + 1) - (clen[pos2] - clen[pos1]));
92 if (eel[i] < MINEEL) { eel[i] = 0.0; }
100 void flushToTempFile() {
101 int gap1 = fr * FLOATSIZE;
102 int gap2 = (nSamples - to) * FLOATSIZE;
105 ftmpOut.seekp(0, ios_base::beg);
106 for (int i = 0; i < cvlen; i++) {
108 ftmpOut.seekp(gap1, ios_base::cur);
109 for (int j = fr; j < to; j++) {
110 ftmpOut.write((char*)p, FLOATSIZE);
113 ftmpOut.seekp(gap2, ios_base::cur);
117 template<class ModelType>
124 calcExpectedEffectiveLengths<ModelType>(model);
126 ftmpOut.open(tmpF, ios::binary);
129 assert(cvlen = M + 1);
131 nSamples = nC * nSpC;
135 size = bufsize_type(nMB) * 1024 * 1024 / FLOATSIZE / cvlen;
136 if (size > (bufsize_type)nSamples) size = nSamples;
138 buffer = new float[size];
143 cvec = new int[cvlen];
144 theta = new double[cvlen];
145 gammas = new distribution_type*[cvlen];
146 rgs = new generator_type*[cvlen];
148 tau_denoms = new double[nSamples];
149 memset(tau_denoms, 0, sizeof(double) * nSamples);
151 double *mw = model.getMW();
152 for (int i = 0; i < nC; i++) {
153 for (int j = 0; j < cvlen; j++) {
157 for (int j = 0; j < cvlen; j++) {
158 gammas[j] = new distribution_type(cvec[j]); // need change back before publishing
159 rgs[j] = new generator_type(engine, *gammas[j]);
162 for (int j = 0; j < nSpC; j++) {
164 for (int k = 0; k < cvlen; k++) {
165 theta[k] = (k == 0 || eel[k] > EPSILON ? (*rgs[k])() : 0.0);
169 for (int k = 0; k < cvlen; k++) theta[k] /= sum;
172 for (int k = 0; k < cvlen; k++) {
173 theta[k] = (mw[k] < EPSILON ? 0.0 : theta[k] / mw[k]);
176 assert(sum >= EPSILON);
177 for (int k = 0; k < cvlen; k++) theta[k] /= sum;
179 *p = (float)theta[0]; ++p;
180 assert(1.0 - theta[0] > 0.0);
181 for (int k = 1; k < cvlen; k++) {
182 if (eel[k] > EPSILON) {
183 theta[k] /= (1.0 - theta[0]);
184 tau_denoms[to] += theta[k] / eel[k];
187 if (theta[k] != 0.0) { fprintf(stderr, "K=%d Theta_K=%lf\n", k, theta[k]); exit(-1); }
190 *p = (float)theta[k];
198 if (verbose) { printf("%d vectors are sampled!\n", to); }
202 for (int j = 0; j < cvlen; j++) {
208 if (fr != to) { flushToTempFile(); }
220 if (verbose) { printf("Sampling is finished!\n"); }
223 void calcCI(int nSamples, float *samples, float &lb, float &ub) {
224 int p, q; // p pointer for lb, q pointer for ub;
226 int threshold = nSamples - (int(confidence * nSamples - 1e-8) + 1);
229 sort(samples, samples + nSamples);
231 p = 0; q = nSamples - 1;
235 while (newq > 0 && samples[newq - 1] == samples[newq]) newq--;
237 } while (newq >= 0 && nSamples - (newq + 1) <= threshold);
239 nOutside = nSamples - (q + 1);
241 lb = -1e30; ub = 1e30;
243 if (samples[q] - samples[p] < ub - lb) {
249 while (newp < nSamples - 1 && samples[newp] == samples[newp + 1]) newp++;
251 if (newp <= threshold) {
252 nOutside += newp - p;
254 while (nOutside > threshold && q < nSamples - 1) {
256 while (newq < nSamples - 1 && samples[newq] == samples[newq + 1]) newq++;
257 nOutside -= newq - q;
260 assert(nOutside <= threshold);
263 } while (p <= threshold);
266 void generateResults(char* imdName) {
267 float *itsamples, *gtsamples;
272 iso_tau = new CIType[M + 1];
273 gene_tau = new CIType[m];
275 itsamples = new float[nSamples];
276 gtsamples = new float[nSamples];
278 fin.open(tmpF, ios::binary);
280 //read sampled theta0 values
281 for (int k = 0; k < nSamples; k++) fin.read((char*)(&itsamples[k]), FLOATSIZE);
283 for (int i = 0; i < m; i++) {
284 int b = gi.spAt(i), e = gi.spAt(i + 1);
285 memset(gtsamples, 0, FLOATSIZE * nSamples);
286 for (int j = b; j < e; j++) {
287 for (int k = 0; k < nSamples; k++) {
288 fin.read((char*)(&itsamples[k]), FLOATSIZE);
289 if (eel[j] > EPSILON && tau_denoms[k] > EPSILON) { itsamples[k] = itsamples[k] / eel[j] / tau_denoms[k]; }
291 if (itsamples[k] != 0.0) { fprintf(stderr, "Not equal to 0! K=%d, Sampled Theta Value=%lf\n", k, itsamples[k]); exit(-1); }
294 gtsamples[k] += itsamples[k];
296 calcCI(nSamples, itsamples, iso_tau[j].lb, iso_tau[j].ub);
298 calcCI(nSamples, gtsamples, gene_tau[i].lb, gene_tau[i].ub);
300 if (verbose && (i + 1) % 1000 == 0) { printf("%d genes are done!\n", i + 1); }
305 //isoform level results
306 sprintf(outF, "%s.iso_res", imdName);
307 fo = fopen(outF, "a");
308 for (int i = 1; i <= M; i++)
309 fprintf(fo, "%.6g%c", iso_tau[i].lb, (i < M ? '\t' : '\n'));
310 for (int i = 1; i <= M; i++)
311 fprintf(fo, "%.6g%c", iso_tau[i].ub, (i < M ? '\t' : '\n'));
315 sprintf(outF, "%s.gene_res", imdName);
316 fo = fopen(outF, "a");
317 for (int i = 0; i < m; i++)
318 fprintf(fo, "%.6g%c", gene_tau[i].lb, (i < m - 1 ? '\t' : '\n'));
319 for (int i = 0; i < m; i++)
320 fprintf(fo, "%.6g%c", gene_tau[i].ub, (i < m - 1 ? '\t' : '\n'));
329 if (verbose) { printf("All credibility intervals are calculated!\n"); }
331 sprintf(outF, "%s.tau_denoms", imdName);
332 fo = fopen(outF, "w");
333 fprintf(fo, "%d\n", nSamples);
334 for (int i = 0; i < nSamples; i++) fprintf(fo, "%.15g ", tau_denoms[i]);
339 int main(int argc, char* argv[]) {
341 printf("Usage: rsem-calculate-credibility-intervals reference_name sample_name sampleToken confidence nSpC nMB[-q]\n");
345 confidence = atof(argv[4]);
346 nSpC = atoi(argv[5]);
350 if (argc > 7 && !strcmp(argv[7], "-q")) {
355 sprintf(imdName, "%s.temp/%s", argv[2], argv[3]);
356 sprintf(statName, "%s.stat/%s", argv[2], argv[3]);
358 sprintf(modelF, "%s.model", statName);
359 FILE *fi = fopen(modelF, "r");
360 if (fi == NULL) { fprintf(stderr, "Cannot open %s!\n", modelF); exit(-1); }
361 fscanf(fi, "%d", &model_type);
364 sprintf(refF, "%s.seq", argv[1]);
365 refs.loadRefs(refF, 1);
367 sprintf(groupF, "%s.grp", argv[1]);
371 sprintf(tmpF, "%s.tmp", imdName);
372 sprintf(cvsF, "%s.countvectors", imdName);
375 case 0 : sampling<SingleModel>(); break;
376 case 1 : sampling<SingleQModel>(); break;
377 case 2 : sampling<PairedEndModel>(); break;
378 case 3 : sampling<PairedEndQModel>(); break;
381 generateResults(imdName);
386 sprintf(command, "rm -f %s", tmpF);
387 int status = system(command);
389 fprintf(stderr, "Cannot delete %s!\n", tmpF);