Add sampling option to alignment-based salmon; changed projection limit

COMBINE-lab · Jun 17, 2015 · 2b57c2e · 2b57c2e
1 parent c49dccd
commit 2b57c2e
Show file tree

Hide file tree

Showing 3 changed files with 21 additions and 9 deletions.
diff --git a/include/TranscriptCluster.hpp b/include/TranscriptCluster.hpp
@@ -93,7 +93,7 @@ class TranscriptCluster {
             if (round % 100 == 0) {
                 std::cerr << "\r\rproject to polytope: " << round;
             }
-            if (round > 50000) {
+            if (round > 5000) {
                 return;
             }
 

diff --git a/src/SalmonQuantify.cpp b/src/SalmonQuantify.cpp
@@ -3114,14 +3114,6 @@ transcript abundance from RNA-seq reads
         optimizer.optimize(experiment, sopt, 0.01, 10000);
         jointLog->info("Finished optimizer");
 
-        if (sopt.useGSOpt) {
-            jointLog->info("Starting Gibbs Sampler");
-            CollapsedGibbsSampler sampler;
-            sampler.sample(experiment, sopt, sopt.numGibbsSamples);
-            jointLog->info("Finished Gibbs Sampler");
-        }
-
-
         free(memOptions);
         size_t tnum{0};
 
@@ -3135,6 +3127,14 @@ transcript abundance from RNA-seq reads
         salmon::utils::writeAbundancesFromCollapsed(
                 sopt, experiment, estFilePath, commentString);
 
+        if (sopt.useGSOpt) {
+            jointLog->info("Starting Gibbs Sampler");
+            CollapsedGibbsSampler sampler;
+            sampler.sample(experiment, sopt, sopt.numGibbsSamples);
+            jointLog->info("Finished Gibbs Sampler");
+        }
+
+
         // Now create a subdirectory for any parameters of interest
         bfs::path paramsDir = outputDirectory / "libParams";
         if (!boost::filesystem::exists(paramsDir)) {

diff --git a/src/SalmonQuantifyAlignments.cpp b/src/SalmonQuantifyAlignments.cpp
@@ -60,6 +60,7 @@ extern "C" {
 #include "spdlog/spdlog.h"
 #include "EquivalenceClassBuilder.hpp"
 #include "CollapsedEMOptimizer.hpp"
+#include "CollapsedGibbsSampler.hpp"
 
 namespace bfs = boost::filesystem;
 using salmon::math::LOG_0;
@@ -904,6 +905,10 @@ int salmonAlignmentQuantify(int argc, char* argv[]) {
                         "Use mass \"banking\" in subsequent epoch of inference.  Rather than re-observing uniquely "
                         "mapped reads, simply remember the ratio of uniquely to ambiguously mapped reads for each "
                         "transcript and distribute the unique mass uniformly throughout the epoch.")
+    ("useGSOpt", po::bool_switch(&(sopt.useGSOpt))->default_value(false), "[*super*-experimental]: After the initial optimization has finished, "
+                "use collapsed Gibbs sampling to refine estimates even further (and obtain variance)")
+    ("numGibbsSamples", po::value<uint32_t>(&(sopt.numGibbsSamples))->default_value(500), "[*super*-experimental]: Number of Gibbs sampling rounds to "
+     		"perform.")
     ("useVBOpt,v", po::bool_switch(&(sopt.useVBOpt))->default_value(false), "Use the Variational Bayesian EM rather than the "
      			"traditional EM algorithm for optimization in the batch passes.");
 
@@ -1167,6 +1172,13 @@ int salmonAlignmentQuantify(int argc, char* argv[]) {
                     salmon::utils::writeAbundancesFromCollapsed(
                         sopt, alnLib, outputFile, commentString);
 
+                    if (sopt.useGSOpt) {
+                        jointLog->info("Starting Gibbs Sampler");
+                        CollapsedGibbsSampler sampler;
+                        sampler.sample(alnLib, sopt, sopt.numGibbsSamples);
+                        jointLog->info("Finished Gibbs Sampler");
+                    }
+
                     /*
                     fmt::print(stderr, "\n\nwriting output \n");
                     salmon::utils::writeAbundances(sopt, alnLib, outputFile, commentString);