12 #include "my_assert.h"
14 #include "SingleRead.h"
15 #include "SingleReadQ.h"
16 #include "PairedEndRead.h"
17 #include "PairedEndReadQ.h"
20 #include "SingleModel.h"
21 #include "SingleQModel.h"
22 #include "PairedEndModel.h"
23 #include "PairedEndQModel.h"
26 #include "GroupInfo.h"
27 #include "Transcript.h"
28 #include "Transcripts.h"
34 const int OFFSITE = 5;
41 Transcripts transcripts;
44 vector<double> theta, counts;
48 char outReadF[2][STRLEN];
50 char refF[STRLEN], groupF[STRLEN], tiF[STRLEN];
51 char geneResF[STRLEN], isoResF[STRLEN];
55 void genOutReadStreams(int type, char *outFN) {
59 sprintf(outReadF[0], "%s.fa", outFN);
63 sprintf(outReadF[0], "%s.fq", outFN);
67 for (int i = 0; i < n_os; i++)
68 sprintf(outReadF[i], "%s_%d.fa", outFN, i + 1);
72 for (int i = 0; i < n_os; i++)
73 sprintf(outReadF[i], "%s_%d.fq", outFN, i + 1);
77 for (int i = 0; i < n_os; i++)
78 os[i] = new ofstream(outReadF[i]);
81 template<class ModelType>
82 void calcExpectedEffectiveLengths(ModelType& model) {
84 double *pdf = NULL, *cdf = NULL, *clen = NULL; // clen[i] = sigma_{j=1}^{i}pdf[i]*(lb+i)
86 model.getGLD().copyTo(pdf, cdf, lb, ub, span);
87 clen = new double[span + 1];
89 for (int i = 1; i <= span; i++) {
90 clen[i] = clen[i - 1] + pdf[i] * (lb + i);
93 eel.assign(M + 1, 0.0);
94 for (int i = 1; i <= M; i++) {
95 int totLen = refs.getRef(i).getTotLen();
96 int fullLen = refs.getRef(i).getFullLen();
97 int pos1 = max(min(totLen - fullLen + 1, ub) - lb, 0);
98 int pos2 = max(min(totLen, ub) - lb, 0);
100 if (pos2 == 0) { eel[i] = 0.0; continue; }
102 eel[i] = fullLen * cdf[pos1] + ((cdf[pos2] - cdf[pos1]) * (totLen + 1) - (clen[pos2] - clen[pos1]));
104 if (eel[i] < MINEEL) { eel[i] = 0.0; }
112 template<class ReadType, class ModelType>
113 void simulate(char* modelF, char* resultsF) {
114 ModelType model(&refs);
121 calcExpectedEffectiveLengths<ModelType>(model);
123 //generate theta vector
124 ifstream fin(resultsF);
128 getline(fin, line); // read the first line, which is just column names
129 for (int i = 1; i <= M; i++) {
132 for (int j = 0; j < OFFSITE; j++) pos = line.find_first_of('\t', pos) + 1;
133 size_t pos2 = line.find_first_of('\t', pos);
134 if (pos2 == string::npos) pos2 = line.length();
135 tpm = atof(line.substr(pos, pos2 - pos).c_str());
136 theta[i] = tpm * eel[i];
139 assert(denom > EPSILON);
141 for (int i = 1; i <= M; i++) theta[i] = theta[i] / denom * (1.0 - theta[0]);
143 READ_INT_TYPE resimulation_count = 0;
146 model.startSimulation(&sampler, theta);
147 for (READ_INT_TYPE i = 0; i < N; i++) {
148 while (!model.simulate(i, read, sid)) { ++resimulation_count; }
149 read.write(n_os, os);
151 if ((i + 1) % 1000000 == 0 && verbose) cout<<"GEN "<< i + 1<< endl;
153 model.finishSimulation();
155 cout<< "Total number of resimulation is "<< resimulation_count<< endl;
158 void calcExpressionValues(const vector<double>& theta, const vector<double>& eel, vector<double>& tpm, vector<double>& fpkm) {
162 //calculate fraction of count over all mappabile reads
164 frac.assign(M + 1, 0.0);
165 for (int i = 1; i <= M; i++)
166 if (eel[i] >= EPSILON) {
170 general_assert(denom > 0, "No alignable reads?!");
171 for (int i = 1; i <= M; i++) frac[i] /= denom;
174 fpkm.assign(M + 1, 0.0);
175 for (int i = 1; i <= M; i++)
176 if (eel[i] >= EPSILON) fpkm[i] = frac[i] * 1e9 / eel[i];
179 tpm.assign(M + 1, 0.0);
181 for (int i = 1; i <= M; i++) denom += fpkm[i];
182 for (int i = 1; i <= M; i++) tpm[i] = fpkm[i] / denom * 1e6;
185 void writeResFiles(char* outFN) {
188 vector<double> fpkm, tpm, isopct;
189 vector<double> glens, gene_eels, gene_counts, gene_tpm, gene_fpkm;
191 for (int i = 1; i <= M; i++)
192 general_assert(eel[i] > EPSILON || counts[i] <= EPSILON, "An isoform whose effecitve length < " + ftos(MINEEL, 6) + " got sampled!");
194 calcExpressionValues(counts, eel, tpm, fpkm);
196 //calculate IsoPct, etc.
197 isopct.assign(M + 1, 0.0);
198 tlens.assign(M + 1, 0);
200 glens.assign(m, 0.0); gene_eels.assign(m, 0.0);
201 gene_counts.assign(m, 0.0); gene_tpm.assign(m, 0.0); gene_fpkm.assign(m, 0.0);
203 for (int i = 0; i < m; i++) {
204 int b = gi.spAt(i), e = gi.spAt(i + 1);
205 for (int j = b; j < e; j++) {
206 const Transcript& transcript = transcripts.getTranscriptAt(j);
207 tlens[j] = transcript.getLength();
209 glens[i] += tlens[j] * tpm[j];
210 gene_eels[i] += eel[j] * tpm[j];
211 gene_counts[i] += counts[j];
212 gene_tpm[i] += tpm[j];
213 gene_fpkm[i] += fpkm[j];
216 if (gene_tpm[i] < EPSILON) continue;
218 for (int j = b; j < e; j++)
219 isopct[j] = tpm[j] / gene_tpm[i];
220 glens[i] /= gene_tpm[i];
221 gene_eels[i] /= gene_tpm[i];
225 sprintf(isoResF, "%s.sim.isoforms.results", outFN);
226 fo = fopen(isoResF, "w");
227 fprintf(fo, "transcript_id\tgene_id\tlength\teffective_length\tcount\tTPM\tFPKM\tIsoPct\n");
228 for (int i = 1; i <= M; i++) {
229 const Transcript& transcript = transcripts.getTranscriptAt(i);
230 fprintf(fo, "%s\t%s\t%d\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n", transcript.getTranscriptID().c_str(), transcript.getGeneID().c_str(), tlens[i],
231 eel[i], counts[i], tpm[i], fpkm[i], isopct[i] * 1e2);
236 sprintf(geneResF, "%s.sim.genes.results", outFN);
237 fo = fopen(geneResF, "w");
238 fprintf(fo, "gene_id\ttranscript_id(s)\tlength\teffective_length\tcount\tTPM\tFPKM\n");
239 for (int i = 0; i < m; i++) {
240 int b = gi.spAt(i), e = gi.spAt(i + 1);
241 const string& gene_id = transcripts.getTranscriptAt(b).getGeneID();
242 fprintf(fo, "%s\t", gene_id.c_str());
243 for (int j = b; j < e; j++) {
244 fprintf(fo, "%s%c", transcripts.getTranscriptAt(j).getTranscriptID().c_str(), (j < e - 1 ? ',' : '\t'));
246 fprintf(fo, "%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n", glens[i], gene_eels[i], gene_counts[i], gene_tpm[i], gene_fpkm[i]);
251 void releaseOutReadStreams() {
252 for (int i = 0; i < n_os; i++) {
253 ((ofstream*)os[i])->close();
258 int main(int argc, char* argv[]) {
262 if (argc != 7 && argc != 8) {
263 printf("Usage: rsem-simulate-reads reference_name estimated_model_file estimated_isoform_results theta0 N output_name [-q]\n");
267 if (argc == 8 && !strcmp(argv[7], "-q")) quiet = true;
271 sprintf(refF, "%s.seq", argv[1]);
274 sprintf(groupF, "%s.grp", argv[1]);
277 sprintf(tiF, "%s.ti", argv[1]);
278 transcripts.readFrom(tiF);
280 //read model type from modelF
281 fi = fopen(argv[2], "r");
282 if (fi == NULL) { fprintf(stderr, "Cannot open %s! It may not exist.\n", argv[2]); exit(-1); }
283 assert(fscanf(fi, "%d", &model_type) == 1);
286 theta.assign(M + 1, 0.0);
287 theta[0] = atof(argv[4]);
290 genOutReadStreams(model_type, argv[6]);
292 counts.assign(M + 1, 0.0);
295 case 0: simulate<SingleRead, SingleModel>(argv[2], argv[3]); break;
296 case 1: simulate<SingleReadQ, SingleQModel>(argv[2], argv[3]); break;
297 case 2: simulate<PairedEndRead, PairedEndModel>(argv[2], argv[3]); break;
298 case 3: simulate<PairedEndReadQ, PairedEndQModel>(argv[2], argv[3]); break;
301 writeResFiles(argv[6]);
302 releaseOutReadStreams();