#include "PairedEndQModel.h"
#include "Refs.h"
+
#include "GroupInfo.h"
+#include "WriteResults.h"
using namespace std;
struct Params {
- int no, nsamples;
- FILE *fo;
- engine_type *engine;
- double *pme_c, *pve_c; //posterior mean and variance vectors on counts
+ int no, nsamples;
+ FILE *fo;
+ engine_type *engine;
+ double *pme_c, *pve_c; //posterior mean and variance vectors on counts
double *pme_tpm, *pme_fpkm;
+
+ double *pve_c_genes, *pve_c_trans;
};
-
struct Item {
int sid;
double conprb;
int nThreads;
int model_type;
-int m, M;
+int M;
READ_INT_TYPE N0, N1;
HIT_INT_TYPE nHits;
double totc;
int BURNIN, NSAMPLES, GAP;
-char imdName[STRLEN], statName[STRLEN];
-char thetaF[STRLEN], ofgF[STRLEN], groupF[STRLEN], refF[STRLEN], modelF[STRLEN];
+char refName[STRLEN], imdName[STRLEN], statName[STRLEN];
+char thetaF[STRLEN], ofgF[STRLEN], refF[STRLEN], modelF[STRLEN];
char cvsF[STRLEN];
Refs refs;
-GroupInfo gi;
vector<HIT_INT_TYPE> s;
vector<Item> hits;
vector<double> eel;
double *mw;
+vector<int> pseudo_counts;
vector<double> pme_c, pve_c; //global posterior mean and variance vectors on counts
vector<double> pme_tpm, pme_fpkm;
-bool var_opt;
bool quiet;
Params *paramsArray;
pthread_attr_t attr;
int rc;
-void load_data(char* reference_name, char* statName, char* imdName) {
+bool hasSeed;
+seedType seed;
+
+int m;
+char groupF[STRLEN];
+GroupInfo gi;
+
+bool alleleS;
+int m_trans;
+GroupInfo gt, ta;
+vector<double> pve_c_genes, pve_c_trans;
+
+void load_data(char* refName, char* statName, char* imdName) {
ifstream fin;
string line;
int tmpVal;
//load reference file
- sprintf(refF, "%s.seq", reference_name);
+ sprintf(refF, "%s.seq", refName);
refs.loadRefs(refF, 1);
M = refs.getM();
- //load groupF
- sprintf(groupF, "%s.grp", reference_name);
- gi.load(groupF);
- m = gi.getm();
-
//load ofgF;
sprintf(ofgF, "%s.ofg", imdName);
fin.open(ofgF);
totc = N0 + N1 + (M + 1);
- if (verbose) { printf("Loading Data is finished!\n"); }
+ if (verbose) { printf("Loading data is finished!\n"); }
}
-template<class ModelType>
-void calcExpectedEffectiveLengths(ModelType& model) {
- int lb, ub, span;
- double *pdf = NULL, *cdf = NULL, *clen = NULL; // clen[i] = sigma_{j=1}^{i}pdf[i]*(lb+i)
+void load_group_info(char* refName) {
+ // Load group info
+ sprintf(groupF, "%s.grp", refName);
+ gi.load(groupF);
+ m = gi.getm();
- model.getGLD().copyTo(pdf, cdf, lb, ub, span);
- clen = new double[span + 1];
- clen[0] = 0.0;
- for (int i = 1; i <= span; i++) {
- clen[i] = clen[i - 1] + pdf[i] * (lb + i);
- }
+ alleleS = isAlleleSpecific(refName, >, &ta); // if allele-specific
+ m_trans = (alleleS ? ta.getm() : 0);
- eel.assign(M + 1, 0.0);
- for (int i = 1; i <= M; i++) {
- int totLen = refs.getRef(i).getTotLen();
- int fullLen = refs.getRef(i).getFullLen();
- int pos1 = max(min(totLen - fullLen + 1, ub) - lb, 0);
- int pos2 = max(min(totLen, ub) - lb, 0);
-
- if (pos2 == 0) { eel[i] = 0.0; continue; }
-
- eel[i] = fullLen * cdf[pos1] + ((cdf[pos2] - cdf[pos1]) * (totLen + 1) - (clen[pos2] - clen[pos1]));
- assert(eel[i] >= 0);
- if (eel[i] < MINEEL) { eel[i] = 0.0; }
- }
-
- delete[] pdf;
- delete[] cdf;
- delete[] clen;
+ if (verbose) { printf("Loading group information is finished!\n"); }
+}
+
+// Load imdName.omit and initialize the pseudo count vector.
+void load_omit_info(const char* imdName) {
+ char omitF[STRLEN];
+ sprintf(omitF, "%s.omit", imdName);
+ FILE *fi = fopen(omitF, "r");
+ pseudo_counts.assign(M + 1, 1);
+ int tid;
+ while (fscanf(fi, "%d", &tid) == 1) pseudo_counts[tid] = 0;
+ fclose(fi);
}
template<class ModelType>
ModelType model;
model.read(modelF);
- calcExpectedEffectiveLengths<ModelType>(model);
+ calcExpectedEffectiveLengths<ModelType>(M, refs, model, eel);
memcpy(mw, model.getMW(), sizeof(double) * (M + 1)); // otherwise, after exiting this procedure, mw becomes undefined
}
paramsArray = new Params[nThreads];
threads = new pthread_t[nThreads];
+ hasSeed ? engineFactory::init(seed) : engineFactory::init();
for (int i = 0; i < nThreads; i++) {
paramsArray[i].no = i;
memset(paramsArray[i].pme_tpm, 0, sizeof(double) * (M + 1));
paramsArray[i].pme_fpkm = new double[M + 1];
memset(paramsArray[i].pme_fpkm, 0, sizeof(double) * (M + 1));
+
+ paramsArray[i].pve_c_genes = new double[m];
+ memset(paramsArray[i].pve_c_genes, 0, sizeof(double) * m);
+
+ paramsArray[i].pve_c_trans = NULL;
+ if (alleleS) {
+ paramsArray[i].pve_c_trans = new double[m_trans];
+ memset(paramsArray[i].pve_c_trans, 0, sizeof(double) * m_trans);
+ }
}
+ engineFactory::finish();
/* set thread attribute to be joinable */
pthread_attr_init(&attr);
theta.assign(M + 1, 0);
denom = 0.0;
for (int i = 0; i <= M; i++) {
- theta[i] = gmg();
+ theta[i] = (pseudo_counts[i] > 0 ? gmg() : 0.0);
denom += theta[i];
}
assert(denom > EPSILON);
fprintf(fo, "%d\n", counts[M]);
}
-void polishTheta(vector<double>& theta, const vector<double>& eel, const double* mw) {
- double sum = 0.0;
-
- /* The reason that for noise gene, mw value is 1 is :
- * currently, all masked positions are for poly(A) sites, which in theory should be filtered out.
- * So the theta0 does not containing reads from any masked position
- */
-
- for (int i = 0; i <= M; i++) {
- // i == 0, mw[i] == 1
- if (i > 0 && (mw[i] < EPSILON || eel[i] < EPSILON)) {
- theta[i] = 0.0;
- continue;
- }
- theta[i] = theta[i] / mw[i];
- sum += theta[i];
- }
- // currently is OK, since no transcript should be masked totally, only the poly(A) tail related part will be masked
- general_assert(sum >= EPSILON, "No effective length is no less than" + ftos(MINEEL, 6) + " !");
- for (int i = 0; i <= M; i++) theta[i] /= sum;
-}
-
-void calcExpressionValues(const vector<double>& theta, const vector<double>& eel, vector<double>& tpm, vector<double>& fpkm) {
- double denom;
- vector<double> frac;
-
- //calculate fraction of count over all mappabile reads
- denom = 0.0;
- frac.assign(M + 1, 0.0);
- for (int i = 1; i <= M; i++)
- if (eel[i] >= EPSILON) {
- frac[i] = theta[i];
- denom += frac[i];
- }
- general_assert(denom >= EPSILON, "No alignable reads?!");
- for (int i = 1; i <= M; i++) frac[i] /= denom;
-
- //calculate FPKM
- fpkm.assign(M + 1, 0.0);
- for (int i = 1; i <= M; i++)
- if (eel[i] >= EPSILON) fpkm[i] = frac[i] * 1e9 / eel[i];
-
- //calculate TPM
- tpm.assign(M + 1, 0.0);
- denom = 0.0;
- for (int i = 1; i <= M; i++) denom += fpkm[i];
- for (int i = 1; i <= M; i++) tpm[i] = fpkm[i] / denom * 1e6;
-}
-
void* Gibbs(void* arg) {
int CHAINLEN;
HIT_INT_TYPE len, fr, to;
Params *params = (Params*)arg;
vector<double> theta, tpm, fpkm;
- vector<int> z, counts;
+ vector<int> z, counts(pseudo_counts);
vector<double> arr;
- uniform01 rg(*params->engine);
+ uniform_01_generator rg(*params->engine, uniform_01_dist());
// generate initial state
sampleTheta(*params->engine, theta);
z.assign(N1, 0);
-
- counts.assign(M + 1, 1); // 1 pseudo count
counts[0] += N0;
for (READ_INT_TYPE i = 0; i < N1; i++) {
if ((ROUND - BURNIN - 1) % GAP == 0) {
writeCountVector(params->fo, counts);
for (int i = 0; i <= M; i++) theta[i] = counts[i] / totc;
- polishTheta(theta, eel, mw);
- calcExpressionValues(theta, eel, tpm, fpkm);
+ polishTheta(M, theta, eel, mw);
+ calcExpressionValues(M, theta, eel, tpm, fpkm);
for (int i = 0; i <= M; i++) {
- params->pme_c[i] += counts[i] - 1;
- params->pve_c[i] += (counts[i] - 1) * (counts[i] - 1);
+ params->pme_c[i] += counts[i] - pseudo_counts[i];
+ params->pve_c[i] += double(counts[i] - pseudo_counts[i]) * (counts[i] - pseudo_counts[i]);
params->pme_tpm[i] += tpm[i];
params->pme_fpkm[i] += fpkm[i];
}
+
+ for (int i = 0; i < m; i++) {
+ int b = gi.spAt(i), e = gi.spAt(i + 1);
+ double count = 0.0;
+ for (int j = b; j < e; j++) count += counts[j] - pseudo_counts[j];
+ params->pve_c_genes[i] += count * count;
+ }
+
+ if (alleleS)
+ for (int i = 0; i < m_trans; i++) {
+ int b = ta.spAt(i), e = ta.spAt(i + 1);
+ double count = 0.0;
+ for (int j = b; j < e; j++) count += counts[j] - pseudo_counts[j];
+ params->pve_c_trans[i] += count * count;
+ }
}
}
pve_c.assign(M + 1, 0);
pme_tpm.assign(M + 1, 0);
pme_fpkm.assign(M + 1, 0);
+
+ pve_c_genes.assign(m, 0);
+ pve_c_trans.clear();
+ if (alleleS) pve_c_trans.assign(m_trans, 0);
+
for (int i = 0; i < nThreads; i++) {
fclose(paramsArray[i].fo);
delete paramsArray[i].engine;
pme_tpm[j] += paramsArray[i].pme_tpm[j];
pme_fpkm[j] += paramsArray[i].pme_fpkm[j];
}
+
+ for (int j = 0; j < m; j++)
+ pve_c_genes[j] += paramsArray[i].pve_c_genes[j];
+
+ if (alleleS)
+ for (int j = 0; j < m_trans; j++)
+ pve_c_trans[j] += paramsArray[i].pve_c_trans[j];
+
delete[] paramsArray[i].pme_c;
delete[] paramsArray[i].pve_c;
delete[] paramsArray[i].pme_tpm;
delete[] paramsArray[i].pme_fpkm;
+
+ delete[] paramsArray[i].pve_c_genes;
+ if (alleleS) delete[] paramsArray[i].pve_c_trans;
}
delete[] paramsArray;
-
for (int i = 0; i <= M; i++) {
pme_c[i] /= NSAMPLES;
- pve_c[i] = (pve_c[i] - NSAMPLES * pme_c[i] * pme_c[i]) / (NSAMPLES - 1);
+ pve_c[i] = (pve_c[i] - double(NSAMPLES) * pme_c[i] * pme_c[i]) / double(NSAMPLES - 1);
+ if (pve_c[i] < 0.0) pve_c[i] = 0.0;
pme_tpm[i] /= NSAMPLES;
pme_fpkm[i] /= NSAMPLES;
}
-}
-
-void writeResults(char* imdName) {
- char outF[STRLEN];
- FILE *fo;
-
- vector<double> isopct;
- vector<double> gene_counts, gene_tpm, gene_fpkm;
-
- //calculate IsoPct, etc.
- isopct.assign(M + 1, 0.0);
- gene_counts.assign(m, 0.0); gene_tpm.assign(m, 0.0); gene_fpkm.assign(m, 0.0);
for (int i = 0; i < m; i++) {
- int b = gi.spAt(i), e = gi.spAt(i + 1);
- for (int j = b; j < e; j++) {
- gene_counts[i] += pme_c[j];
- gene_tpm[i] += pme_tpm[j];
- gene_fpkm[i] += pme_fpkm[j];
- }
- if (gene_tpm[i] < EPSILON) continue;
- for (int j = b; j < e; j++)
- isopct[j] = pme_tpm[j] / gene_tpm[i];
+ int b = gi.spAt(i), e = gi.spAt(i + 1);
+ double pme_c_gene = 0.0;
+ for (int j = b; j < e; j++) pme_c_gene += pme_c[j];
+ pve_c_genes[i] = (pve_c_genes[i] - double(NSAMPLES) * pme_c_gene * pme_c_gene) / double(NSAMPLES - 1);
+ if (pve_c_genes[i] < 0.0) pve_c_genes[i] = 0.0;
}
- //isoform level results
- sprintf(outF, "%s.iso_res", imdName);
- fo = fopen(outF, "a");
- general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
-
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.2f%c", pme_c[i], (i < M ? '\t' : '\n'));
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.2f%c", pme_tpm[i], (i < M ? '\t' : '\n'));
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.2f%c", pme_fpkm[i], (i < M ? '\t' : '\n'));
- for (int i = 1; i <= M; i++)
- fprintf(fo, "%.2f%c", isopct[i] * 1e2, (i < M ? '\t' : '\n'));
- fclose(fo);
-
- //gene level results
- sprintf(outF, "%s.gene_res", imdName);
- fo = fopen(outF, "a");
- general_assert(fo != NULL, "Cannot open " + cstrtos(outF) + "!");
-
- for (int i = 0; i < m; i++)
- fprintf(fo, "%.2f%c", gene_counts[i], (i < m - 1 ? '\t' : '\n'));
- for (int i = 0; i < m; i++)
- fprintf(fo, "%.2f%c", gene_tpm[i], (i < m - 1 ? '\t' : '\n'));
- for (int i = 0; i < m; i++)
- fprintf(fo, "%.2f%c", gene_fpkm[i], (i < m - 1 ? '\t' : '\n'));
- fclose(fo);
-
- if (verbose) { printf("Gibbs based expression values are written!\n"); }
+ if (alleleS)
+ for (int i = 0; i < m_trans; i++) {
+ int b = ta.spAt(i), e = ta.spAt(i + 1);
+ double pme_c_tran = 0.0;
+ for (int j = b; j < e; j++) pme_c_tran += pme_c[j];
+ pve_c_trans[i] = (pve_c_trans[i] - double(NSAMPLES) * pme_c_tran * pme_c_tran) / double(NSAMPLES - 1);
+ if (pve_c_trans[i] < 0.0) pve_c_trans[i] = 0.0;
+ }
}
int main(int argc, char* argv[]) {
if (argc < 7) {
- printf("Usage: rsem-run-gibbs reference_name imdName statName BURNIN NSAMPLES GAP [-p #Threads] [--var] [-q]\n");
+ printf("Usage: rsem-run-gibbs reference_name imdName statName BURNIN NSAMPLES GAP [-p #Threads] [--seed seed] [-q]\n");
exit(-1);
}
+ strcpy(refName, argv[1]);
strcpy(imdName, argv[2]);
strcpy(statName, argv[3]);
GAP = atoi(argv[6]);
nThreads = 1;
- var_opt = false;
+ hasSeed = false;
quiet = false;
for (int i = 7; i < argc; i++) {
if (!strcmp(argv[i], "-p")) nThreads = atoi(argv[i + 1]);
- if (!strcmp(argv[i], "--var")) var_opt = true;
+ if (!strcmp(argv[i], "--seed")) {
+ hasSeed = true;
+ int len = strlen(argv[i + 1]);
+ seed = 0;
+ for (int k = 0; k < len; k++) seed = seed * 10 + (argv[i + 1][k] - '0');
+ }
if (!strcmp(argv[i], "-q")) quiet = true;
}
verbose = !quiet;
printf("Warning: Number of samples is less than number of threads! Change the number of threads to %d!\n", nThreads);
}
- load_data(argv[1], statName, imdName);
+ load_data(refName, statName, imdName);
+ load_group_info(refName);
+ load_omit_info(imdName);
sprintf(modelF, "%s.model", statName);
FILE *fi = fopen(modelF, "r");
if (verbose) printf("Gibbs finished!\n");
- writeResults(imdName);
-
- if (var_opt) {
- char varF[STRLEN];
-
- sprintf(varF, "%s.var", statName);
- FILE *fo = fopen(varF, "w");
- general_assert(fo != NULL, "Cannot open " + cstrtos(varF) + "!");
- for (int i = 0; i < m; i++) {
- int b = gi.spAt(i), e = gi.spAt(i + 1), number_of_isoforms = e - b;
- for (int j = b; j < e; j++) {
- fprintf(fo, "%s\t%d\t%.15g\t%.15g\n", refs.getRef(j).getName().c_str(), number_of_isoforms, pme_c[j], pve_c[j]);
- }
- }
- fclose(fo);
- }
+ writeResultsGibbs(M, m, m_trans, gi, gt, ta, alleleS, imdName, pme_c, pme_fpkm, pme_tpm, pve_c, pve_c_genes, pve_c_trans);
delete mw; // delete the copy
projects / rsem.git / blobdiff