From 79569b0447d0818ebcb16c501c81a683bba9caec Mon Sep 17 00:00:00 2001
From: Mervin Fansler <mmfansler@gmail.com>
Date: Thu, 25 Jun 2020 23:59:54 -0400
Subject: [PATCH] adds strand-preserving bam seqs

---
 src/ProcessReads.cpp | 36 +++++++++++++++++++++++++++---------
 src/ProcessReads.h   |  1 +
 2 files changed, 28 insertions(+), 9 deletions(-)

diff --git a/src/ProcessReads.cpp b/src/ProcessReads.cpp
index 8c9f0353..b1878cba 100644
--- a/src/ProcessReads.cpp
+++ b/src/ProcessReads.cpp
@@ -2979,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) {
@@ -3038,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);
 
diff --git a/src/ProcessReads.h b/src/ProcessReads.h
index 32501194..0b3ad68c 100644
--- a/src/ProcessReads.h
+++ b/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 {