Provide optional strand-specificity flags for kallisto bus

src/MinCollector.cpp

@@ -1,5 +1,6 @@
 #include "MinCollector.h"
 #include <algorithm>
+#include <limits>
 
 // utility functions
 

src/ProcessReads.cpp

@@ -1482,15 +1482,43 @@ void BUSProcessor::processBuffer() {
     u.clear();
 
     // process 2nd read
-    index.match(seq,seqlen, v);
+    index.match(seq, seqlen, v);
 
     // collect the target information
     int ec = -1;
-    int r = tc.intersectKmers(v, v2, false,u);
-    if (!u.empty()) {      
-      ec = tc.findEC(u);
+    int r = tc.intersectKmers(v, v2, false, u);
+
+    // optional strand-specificity
+    if (mp.opt.strand_specific && !u.empty()) {
+      int p = -1;
+      Kmer km;
+      KmerEntry val;
+
+      if (!v.empty()) {
+        vtmp.clear();
+        bool firstStrand = (mp.opt.strand == ProgramOptions::StrandType::FR);
+        p = findFirstMappingKmer(v, val);
+        km = Kmer((seq+p));
+        bool strand = (val.isFw() == (km == km.rep())); // k-mer maps to fw strand?
+        // might need to optimize this
+        const auto &c = index.dbGraph.contigs[val.contig];
+        for (auto tr : u) {
+          for (auto ctx : c.transcripts) {
+            if (tr == ctx.trid) {
+              if ((strand == ctx.sense) == firstStrand) {
+                // swap out
+                vtmp.push_back(tr);
+              }
+              break;
+            }
+          }
+        }
+        if (vtmp.size() < u.size()) {
+          u = vtmp; // copy
+        }
+      }
     }
-
+    
     // find the ec
     if (!u.empty()) {
       BUSData b;      
@@ -1506,7 +1534,7 @@ void BUSProcessor::processBuffer() {
         b.flags = (uint32_t) flags[i / jmax];
       }
 
-      //ec = tc.findEC(u);
+      ec = tc.findEC(u);
 
       // count the pseudoalignment
       if (ec == -1 || ec >= counts.size()) {
@@ -2951,7 +2979,10 @@ FastqSequenceReader::FastqSequenceReader(FastqSequenceReader&& o) :
   o.state = false;
 }
 
-const std::string BamSequenceReader::seq_enc = "=ACMGRSVTWYHKDBN";
+// sequence decoding table
+const std::string BamSequenceReader::seq_enc  = "=ACMGRSVTWYHKDBN";
+// complement sequence decoding table
+const std::string BamSequenceReader::cseq_enc = "=TGKCYSBAWRDMHVN";
 
 BamSequenceReader::~BamSequenceReader() {
   if (fp) {
@@ -3010,16 +3041,31 @@ bool BamSequenceReader::fetchSequences(char *buf, const int limit, std::vector<s
           bufpos += l_bc + l_umi + 1;
 
           pi = buf + bufpos;
-          int len = (l_seq + 1) / 2;
-          if (l_seq % 2) --len;
+          int len = l_seq / 2;
           int j = 0;
-          for (int i = 0; i < len; ++i, ++eseq) {
-            buf[bufpos++] = seq_enc[*eseq >> 4];
-            buf[bufpos++] = seq_enc[*eseq & 0x0F];
-          }
-          if (l_seq % 2) {
-            buf[bufpos++] = seq_enc[*eseq >> 4];
+
+          if (rec->core.flag & BAM_FREVERSE) { // insert reverse complement
+            // start at end
+            bufpos += l_seq;
+            for (int i = 0; i < len; ++i, ++eseq) {
+              buf[--bufpos] = cseq_enc[*eseq >> 4];
+              buf[--bufpos] = cseq_enc[*eseq & 0x0F];
+            }
+            if (l_seq % 2) {
+              buf[--bufpos] = cseq_enc[*eseq >> 4];
+            }
+            // reset position to end
+            bufpos += l_seq;
+          } else { // insert forward sequence
+            for (int i = 0; i < len; ++i, ++eseq) {
+              buf[bufpos++] = seq_enc[*eseq >> 4];
+              buf[bufpos++] = seq_enc[*eseq & 0x0F];
+            }
+            if (l_seq % 2) {
+              buf[bufpos++] = seq_enc[*eseq >> 4];
+            }
           }
+          // insert delimiter
           buf[bufpos++] = '\0';
           seqs.emplace_back(pi, l_seq);
 

src/ProcessReads.h

@@ -165,6 +165,7 @@ class BamSequenceReader : public SequenceReader {
 
 private:
   static const std::string seq_enc;
+  static const std::string cseq_enc;
 };
 
 class MasterProcessor {

src/main.cpp

@@ -546,6 +546,8 @@ void ListSingleCellTechnologies() {
 void ParseOptionsBus(int argc, char **argv, ProgramOptions& opt) {
   int verbose_flag = 0;
   int gbam_flag = 0;
+  int strand_FR_flag = 0;
+  int strand_RF_flag = 0;
 
   const char *opt_string = "i:o:x:t:lbng:c:";
   static struct option long_options[] = {
@@ -560,6 +562,8 @@ void ParseOptionsBus(int argc, char **argv, ProgramOptions& opt) {
     {"genomebam", no_argument, &gbam_flag, 1},
     {"gtf", required_argument, 0, 'g'},
     {"chromosomes", required_argument, 0, 'c'},
+    {"fr-stranded", no_argument, &strand_FR_flag, 1},
+    {"rf-stranded", no_argument, &strand_RF_flag, 1},
     {0,0,0,0}
   };
 
@@ -631,6 +635,15 @@ void ParseOptionsBus(int argc, char **argv, ProgramOptions& opt) {
     opt.genomebam = true;    
   }
 
+  if (strand_FR_flag) {
+    opt.strand_specific = true;
+    opt.strand = ProgramOptions::StrandType::FR;
+  }
+
+  if (strand_RF_flag) {
+    opt.strand_specific = true;
+    opt.strand = ProgramOptions::StrandType::RF;
+  }
 
   // all other arguments are fast[a/q] files to be read
   for (int i = optind; i < argc; i++) {
@@ -1725,6 +1738,8 @@ void usageBus() {
        << "-t, --threads=INT             Number of threads to use (default: 1)" << endl
        << "-b, --bam                     Input file is a BAM file" << endl
        << "-n, --num                     Output number of read in flag column (incompatible with --bam)" << endl
+       << "    --fr-stranded             Strand specific reads, sequence read forward" << endl
+       << "    --rf-stranded             Strand specific reads, sequence read reverse" << endl
        << "    --verbose                 Print out progress information every 1M proccessed reads" << endl;
 }