Allowed for a subset of reference sequences to be declared in an input SAM/BAM file...

[rsem.git] / Gibbs.cpp

diff --git a/Gibbs.cpp b/Gibbs.cpp
index d5e0b9c49634e29c9a4f2b5ed2a5c74026f5642f..7563dd49f508129992d6a52aa70ed3a8c3130b40 100644 (file)
--- a/Gibbs.cpp
+++ b/Gibbs.cpp
@@ -64,6 +64,7 @@ vector<Item> hits;
 vector<double> eel;
 double *mw;
 
 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;
 
 vector<double> pme_c, pve_c; //global posterior mean and variance vectors on counts
 vector<double> pme_tpm, pme_fpkm;
 


@@ -86,18 +87,6 @@ int m_trans;
 GroupInfo gt, ta;
 vector<double> pve_c_genes, pve_c_trans;
 
 GroupInfo gt, ta;
 vector<double> pve_c_genes, pve_c_trans;
 

-void load_group_info(char* refName) {
-  // Load group info
-  sprintf(groupF, "%s.grp", refName);
-  gi.load(groupF);
-  m = gi.getm();
-  
-  alleleS = isAlleleSpecific(refName, &gt, &ta); // if allele-specific 
-  m_trans = (alleleS ? ta.getm() : 0);
-
-  if (verbose) { printf("Loading group information is finished!\n"); }
-}
-

 void load_data(char* refName, char* statName, char* imdName) {
        ifstream fin;
        string line;
 void load_data(char* refName, char* statName, char* imdName) {
        ifstream fin;
        string line;


@@ -138,6 +127,29 @@ void load_data(char* refName, char* statName, char* imdName) {
        if (verbose) { printf("Loading data is finished!\n"); }
 }
 
        if (verbose) { printf("Loading data is finished!\n"); }
 }
 

+void load_group_info(char* refName) {
+  // Load group info
+  sprintf(groupF, "%s.grp", refName);
+  gi.load(groupF);
+  m = gi.getm();
+  
+  alleleS = isAlleleSpecific(refName, &gt, &ta); // if allele-specific 
+  m_trans = (alleleS ? ta.getm() : 0);
+
+  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>
 void init_model_related(char* modelF) {
        ModelType model;
 template<class ModelType>
 void init_model_related(char* modelF) {
        ModelType model;


@@ -206,7 +218,7 @@ void sampleTheta(engine_type& engine, vector<double>& theta) {
        theta.assign(M + 1, 0);
        denom = 0.0;
        for (int i = 0; i <= M; i++) {
        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);
                denom += theta[i];
        }
        assert(denom > EPSILON);


@@ -226,7 +238,7 @@ void* Gibbs(void* arg) {
        Params *params = (Params*)arg;
 
        vector<double> theta, tpm, fpkm;
        Params *params = (Params*)arg;
 
        vector<double> theta, tpm, fpkm;

-       vector<int> z, counts;
+       vector<int> z, counts(pseudo_counts);

        vector<double> arr;
 
        uniform_01_generator rg(*params->engine, uniform_01_dist());
        vector<double> arr;
 
        uniform_01_generator rg(*params->engine, uniform_01_dist());


@@ -235,8 +247,6 @@ void* Gibbs(void* arg) {
        sampleTheta(*params->engine, theta);
 
        z.assign(N1, 0);
        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++) {
        counts[0] += N0;
 
        for (READ_INT_TYPE i = 0; i < N1; i++) {


@@ -274,8 +284,8 @@ void* Gibbs(void* arg) {
                                polishTheta(M, theta, eel, mw);
                                calcExpressionValues(M, theta, eel, tpm, fpkm);
                                for (int i = 0; i <= M; i++) {
                                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] += double(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];
                                }
                                        params->pme_tpm[i] += tpm[i];
                                        params->pme_fpkm[i] += fpkm[i];
                                }


@@ -283,7 +293,7 @@ void* Gibbs(void* arg) {
                                for (int i = 0; i < m; i++) {
                                  int b = gi.spAt(i), e = gi.spAt(i + 1);
                                  double count = 0.0;
                                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] - 1;
+                                 for (int j = b; j < e; j++) count += counts[j] - pseudo_counts[j];

                                  params->pve_c_genes[i] += count * count;
                                }
 
                                  params->pve_c_genes[i] += count * count;
                                }
 


@@ -291,7 +301,7 @@ void* Gibbs(void* arg) {
                                  for (int i = 0; i < m_trans; i++) {
                                    int b = ta.spAt(i), e = ta.spAt(i + 1);
                                    double count = 0.0;
                                  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] - 1;
+                                   for (int j = b; j < e; j++) count += counts[j] - pseudo_counts[j];

                                    params->pve_c_trans[i] += count * count;
                                  }
                        }
                                    params->pve_c_trans[i] += count * count;
                                  }
                        }


@@ -410,8 +420,9 @@ int main(int argc, char* argv[]) {
                printf("Warning: Number of samples is less than number of threads! Change the number of threads to %d!\n", nThreads);
        }
 
                printf("Warning: Number of samples is less than number of threads! Change the number of threads to %d!\n", nThreads);
        }
 

-       load_group_info(refName);

        load_data(refName, 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");
 
        sprintf(modelF, "%s.model", statName);
        FILE *fi = fopen(modelF, "r");