fix automated strand detection in single-end alignment-based moded

COMBINE-lab · Aug 24, 2018 · a02c0e3 · a02c0e3
1 parent 2bd7786
commit a02c0e3
Show file tree

Hide file tree

Showing 4 changed files with 15 additions and 12 deletions.
diff --git a/include/BAMQueue.tpp b/include/BAMQueue.tpp
@@ -553,6 +553,9 @@ inline bool BAMQueue<FragT>::getFrag_(UnpairedRead& sread, FilterT filt) {
             !(bam_flag(sread.read) & BAM_FUNMAP) and
             sread.transcriptID() >= 0) {
             haveValidRead = true;
+            bool isFwd = !(bam_strand(sread.read));
+            int32_t startPos = bam_pos(sread.read); 
+            sread.libFmt = salmon::utils::hitType(startPos, isFwd);
         }
 
         if (!haveValidRead) { 

diff --git a/include/UnpairedRead.hpp b/include/UnpairedRead.hpp
@@ -20,7 +20,7 @@ extern "C" {
 struct UnpairedRead {
   bam_seq_t* read = nullptr;
   double logProb;
-  LibraryFormat libFmt{ReadType::PAIRED_END, ReadOrientation::NONE,
+  LibraryFormat libFmt{ReadType::SINGLE_END, ReadOrientation::NONE,
                        ReadStrandedness::U};
 
   UnpairedRead()

diff --git a/src/SalmonQuantify.cpp b/src/SalmonQuantify.cpp
@@ -219,7 +219,6 @@ void processMiniBatch(ReadExperimentT& readExp, ForgettingMassCalculator& fmCalc
     incompatPrior = salmon::math::LOG_1;
   }
 
-  auto expectedLibraryFormat = readLib.format();
   uint64_t zeroProbFrags{0};
 
   // EQClass
@@ -237,8 +236,8 @@ void processMiniBatch(ReadExperimentT& readExp, ForgettingMassCalculator& fmCalc
   double startingCumulativeMass =
       fmCalc.cumulativeLogMassAt(firstTimestepOfRound);
 
-  auto isUnexpectedOrphan = [expectedLibraryFormat](AlnT& aln) -> bool {
-    return (expectedLibraryFormat.type == ReadType::PAIRED_END and
+  auto isUnexpectedOrphan = [](AlnT& aln, LibraryFormat expectedLibFormat) -> bool {
+    return (expectedLibFormat.type == ReadType::PAIRED_END and
             aln.mateStatus != rapmap::utils::MateStatus::PAIRED_END_PAIRED);
   };
 
@@ -253,6 +252,7 @@ void processMiniBatch(ReadExperimentT& readExp, ForgettingMassCalculator& fmCalc
       if (alnGroup.size() == 0) {
         continue;
       }
+      LibraryFormat expectedLibraryFormat = readLib.format();
       std::fill(libTypeCountsPerFrag.begin(), libTypeCountsPerFrag.end(), 0);
 
       // We start out with probability 0
@@ -340,7 +340,7 @@ void processMiniBatch(ReadExperimentT& readExp, ForgettingMassCalculator& fmCalc
           double logFragProb = LOG_1;
           // If we are expecting a paired-end library, and this is an orphan,
           // then logFragProb should be small
-          if (isUnexpectedOrphan(aln)) {
+          if (isUnexpectedOrphan(aln, expectedLibraryFormat)) {
             logFragProb = LOG_EPSILON;
           }
 

diff --git a/src/SalmonQuantifyAlignments.cpp b/src/SalmonQuantifyAlignments.cpp
@@ -185,7 +185,7 @@ void processMiniBatch(AlignmentLibraryT<FragT>& alnLib,
 
   double startingCumulativeMass =
       fmCalc.cumulativeLogMassAt(firstTimestepOfRound);
-  auto expectedLibraryFormat = alnLib.format();
+  LibraryFormat expectedLibraryFormat = alnLib.format();
   uint32_t numBurninFrags{salmonOpts.numBurninFrags};
   bool noLengthCorrection{salmonOpts.noLengthCorrection};
 
@@ -200,8 +200,8 @@ void processMiniBatch(AlignmentLibraryT<FragT>& alnLib,
 
   double maxZeroFrac{0.0};
 
-  auto isUnexpectedOrphan = [expectedLibraryFormat](FragT* aln) -> bool {
-    return (expectedLibraryFormat.type == ReadType::PAIRED_END and
+  auto isUnexpectedOrphan = [](FragT* aln, LibraryFormat expectedLibFormat) -> bool {
+    return (expectedLibFormat.type == ReadType::PAIRED_END and
             !aln->isPaired());
   };
 
@@ -226,7 +226,7 @@ void processMiniBatch(AlignmentLibraryT<FragT>& alnLib,
 
     // If we actually got some work
     if (miniBatch != nullptr) {
-
+      expectedLibraryFormat = alnLib.format();
       useAuxParams = (processedReads >= salmonOpts.numPreBurninFrags);
       bool considerCondProb = (useAuxParams or burnedIn);
       ++activeBatches;
@@ -288,7 +288,7 @@ void processMiniBatch(AlignmentLibraryT<FragT>& alnLib,
             double logFragProb = LOG_1;
             // If we are expecting a paired-end library, and this is an orphan,
             // then logFragProb should be small
-            if (isUnexpectedOrphan(aln)) {
+            if (isUnexpectedOrphan(aln, expectedLibraryFormat)) {
               logFragProb = LOG_EPSILON;
             }
 
@@ -335,11 +335,11 @@ void processMiniBatch(AlignmentLibraryT<FragT>& alnLib,
               detector->addSample(aln->libFormat());
               if (detector->canGuess()) {
                 detector->mostLikelyType(alnLib.getFormat());
-                expectedLibraryFormat = alnLib.getFormat();
+                expectedLibraryFormat = alnLib.format();
                 incompatPrior = salmonOpts.incompatPrior;
                 autoDetect = false;
               } else if (!detector->isActive()) {
-                expectedLibraryFormat = alnLib.getFormat();
+                expectedLibraryFormat = alnLib.format();
                 incompatPrior = salmonOpts.incompatPrior;
                 autoDetect = false;
               }