simulation.cpp

   1 #include<cmath>
   2 #include<cstdio>
   3 #include<cstring>
   4 #include<cstdlib>
   5 #include<cassert>
   6 #include<string>
   7 #include<iostream>
   8 #include<fstream>
   9 #include<vector>
  10
  11 #include "utils.h"
  12
  13 #include "Read.h"
  14 #include "SingleRead.h"
  15 #include "SingleReadQ.h"
  16 #include "PairedEndRead.h"
  17 #include "PairedEndReadQ.h"
  18
  19 #include "Model.h"
  20 #include "SingleModel.h"
  21 #include "SingleQModel.h"
  22 #include "PairedEndModel.h"
  23 #include "PairedEndQModel.h"
  24
  25 #include "Refs.h"
  26 #include "GroupInfo.h"
  27 #include "Transcript.h"
  28 #include "Transcripts.h"
  29
  30 #include "simul.h"
  31
  32 using namespace std;
  33
  34 READ_INT_TYPE N;
  35 int model_type, M, m;
  36
  37 Refs refs;
  38 GroupInfo gi;
  39 Transcripts transcripts;
  40
  41 double *theta, *counts;
  42 vector<double> eel;
  43
  44 int n_os;
  45 ostream *os[2];
  46 char outReadF[2][STRLEN];
  47
  48 char refF[STRLEN], groupF[STRLEN], tiF[STRLEN];
  49 char geneResF[STRLEN], isoResF[STRLEN];
  50
  51 simul sampler;
  52
  53 void genOutReadStreams(int type, char *outFN) {
  54         switch(type) {
  55         case 0 :
  56                 n_os = 1;
  57                 sprintf(outReadF[0], "%s.fa", outFN);
  58                 break;
  59         case 1 :
  60                 n_os = 1;
  61                 sprintf(outReadF[0], "%s.fq", outFN);
  62                 break;
  63         case 2 :
  64                 n_os = 2;
  65                 for (int i = 0; i < n_os; i++)
  66                         sprintf(outReadF[i], "%s_%d.fa", outFN, i + 1);
  67                 break;
  68         case 3 :
  69                 n_os = 2;
  70                 for (int i = 0; i < n_os; i++)
  71                         sprintf(outReadF[i], "%s_%d.fq", outFN, i + 1);
  72                 break;
  73         }
  74
  75         for (int i = 0; i < n_os; i++)
  76                 os[i] = new ofstream(outReadF[i]);
  77 }
  78
  79 template<class ModelType>
  80 void calcExpectedEffectiveLengths(ModelType& model) {
  81   int lb, ub, span;
  82   double *pdf = NULL, *cdf = NULL, *clen = NULL; // clen[i] = sigma_{j=1}^{i}pdf[i]*(lb+i)
  83
  84   model.getGLD().copyTo(pdf, cdf, lb, ub, span);
  85   clen = new double[span + 1];
  86   clen[0] = 0.0;
  87   for (int i = 1; i <= span; i++) {
  88     clen[i] = clen[i - 1] + pdf[i] * (lb + i);
  89   }
  90
  91   eel.clear();
  92   eel.resize(M + 1, 0.0);
  93   for (int i = 1; i <= M; i++) {
  94     int totLen = refs.getRef(i).getTotLen();
  95     int fullLen = refs.getRef(i).getFullLen();
  96     int pos1 = max(min(totLen - fullLen + 1, ub) - lb, 0);
  97     int pos2 = max(min(totLen, ub) - lb, 0);
  98
  99     if (pos2 == 0) { eel[i] = 0.0; continue; }
 100
 101     eel[i] = fullLen * cdf[pos1] + ((cdf[pos2] - cdf[pos1]) * (totLen + 1) - (clen[pos2] - clen[pos1]));
 102     assert(eel[i] >= 0);
 103     if (eel[i] < MINEEL) { eel[i] = 0.0; }
 104   }
 105
 106   delete[] pdf;
 107   delete[] cdf;
 108   delete[] clen;
 109 }
 110
 111 template<class ReadType, class ModelType>
 112 void simulate(char* modelF, char* resultsF) {
 113         ModelType model(&refs);
 114         ReadType read;
 115         int sid;
 116
 117         model.read(modelF);
 118
 119         //calculate eel
 120         calcExpectedEffectiveLengths<ModelType>(model);
 121
 122         //generate theta vector
 123         ifstream fin(resultsF);
 124         string line;
 125         double tau;
 126         double denom = 0.0;
 127         for (int i = 1; i <= M; i++) {
 128           getline(fin, line);
 129           size_t pos = 0;
 130           for (int j = 0; j < 2; j++) pos = line.find_first_of('\t', pos) + 1;
 131           size_t pos2 = line.find_first_of('\t', pos);
 132           if (pos2 == string::npos) pos2 = line.length();
 133           tau = atof(line.substr(pos, pos2 - pos).c_str());
 134           theta[i] = tau * eel[i];
 135           denom += theta[i];
 136         }
 137         assert(denom > EPSILON);
 138         fin.close();
 139         for (int i = 1; i <= M; i++) theta[i] = theta[i] / denom * (1.0 - theta[0]);
 140
 141         READ_INT_TYPE resimulation_count = 0;
 142
 143         //simulating...
 144         model.startSimulation(&sampler, theta);
 145         for (READ_INT_TYPE i = 0; i < N; i++) {
 146                 while (!model.simulate(i, read, sid)) { ++resimulation_count; }
 147                 read.write(n_os, os);
 148                 ++counts[sid];
 149                 if ((i + 1) % 1000000 == 0 && verbose) cout<<"GEN "<< i + 1<< endl;
 150         }
 151         model.finishSimulation();
 152
 153         cout<< "Total number of resimulation is "<< resimulation_count<< endl;
 154 }
 155
 156 void writeResFiles(char* outFN) {
 157         FILE *fo;
 158         double denom;
 159
 160         //calculate tau values
 161         double *tau = new double[M + 1];
 162         memset(tau, 0, sizeof(double) * (M + 1));
 163         denom = 0.0;
 164         for (int i = 1; i <= M; i++)
 165                 if (eel[i] > EPSILON) {
 166                         tau[i] = counts[i] / eel[i];
 167                         denom += tau[i];
 168                 }
 169                 else {
 170                     if (counts[i] > EPSILON) { printf("Warning: An isoform which EEL is less than %.6g gets sampled!\n", MINEEL); }
 171                 }
 172         assert(denom > 0.0);
 173         for (int i = 1; i <= M; i++) tau[i] /= denom;
 174
 175         //isoform level
 176         sprintf(isoResF, "%s.sim.isoforms.results", outFN);
 177         fo = fopen(isoResF, "w");
 178         for (int i = 1; i <= M; i++) {
 179                 const Transcript& transcript = transcripts.getTranscriptAt(i);
 180                 fprintf(fo, "%s\t%.2f\t%.15g", transcript.getTranscriptID().c_str(), counts[i], tau[i]);
 181
 182                 if (transcript.getLeft() != "") { fprintf(fo, "\t%s", transcript.getLeft().c_str()); }
 183                 fprintf(fo, "\n");
 184         }
 185         fclose(fo);
 186
 187         //gene level
 188         sprintf(geneResF, "%s.sim.genes.results", outFN);
 189         fo = fopen(geneResF, "w");
 190         for (int i = 0; i < m; i++) {
 191           double sum_c = 0.0, sum_tau = 0.0;
 192                 int b = gi.spAt(i), e = gi.spAt(i + 1);
 193                 for (int j = b; j < e; j++) {
 194                         sum_c += counts[j];
 195                         sum_tau += tau[j];
 196                 }
 197                 const string& gene_id = transcripts.getTranscriptAt(b).getGeneID();
 198                 fprintf(fo, "%s\t%.2f\t%.15g\t", gene_id.c_str(), sum_c, sum_tau);
 199                 for (int j = b; j < e; j++) {
 200                         fprintf(fo, "%s%c", transcripts.getTranscriptAt(j).getTranscriptID().c_str(), (j < e - 1 ? ',' : '\n'));
 201                 }
 202         }
 203         fclose(fo);
 204
 205         delete[] tau;
 206 }
 207
 208 void releaseOutReadStreams() {
 209         for (int i = 0; i < n_os; i++) {
 210                 ((ofstream*)os[i])->close();
 211                 delete os[i];
 212         }
 213 }
 214
 215 int main(int argc, char* argv[]) {
 216         bool quiet = false;
 217         FILE *fi = NULL;
 218
 219         if (argc != 7 && argc != 8) {
 220                 printf("Usage: rsem-simulate-reads reference_name estimated_model_file estimated_isoform_results theta0 N output_name [-q]\n");
 221                 exit(-1);
 222         }
 223
 224         if (argc == 8 && !strcmp(argv[7], "-q")) quiet = true;
 225         verbose = !quiet;
 226
 227         //load basic files
 228         sprintf(refF, "%s.seq", argv[1]);
 229         refs.loadRefs(refF);
 230         M = refs.getM();
 231         sprintf(groupF, "%s.grp", argv[1]);
 232         gi.load(groupF);
 233         m = gi.getm();
 234         sprintf(tiF, "%s.ti", argv[1]);
 235         transcripts.readFrom(tiF);
 236
 237         //read model type from modelF
 238         fi = fopen(argv[2], "r");
 239         if (fi == NULL) { fprintf(stderr, "Cannot open %s! It may not exist.\n", argv[2]); exit(-1); }
 240         assert(fscanf(fi, "%d", &model_type) == 1);
 241         fclose(fi);
 242
 243         theta = new double[M + 1];
 244         theta[0] = atof(argv[4]);
 245         N = atoi(argv[5]);
 246
 247         genOutReadStreams(model_type, argv[6]);
 248
 249         counts = new double[M + 1];
 250         memset(counts, 0, sizeof(double) * (M + 1));
 251
 252         switch(model_type) {
 253         case 0: simulate<SingleRead, SingleModel>(argv[2], argv[3]); break;
 254         case 1: simulate<SingleReadQ, SingleQModel>(argv[2], argv[3]); break;
 255         case 2: simulate<PairedEndRead, PairedEndModel>(argv[2], argv[3]); break;
 256         case 3: simulate<PairedEndReadQ, PairedEndQModel>(argv[2], argv[3]); break;
 257         }
 258
 259         writeResFiles(argv[6]);
 260         releaseOutReadStreams();
 261
 262         delete[] theta;
 263         delete[] counts;
 264
 265         return 0;
 266 }