Fix SAM output's backtrace for nucl-prot backtraces and show rev-comp…

… seq correctly #845
soedinglab · May 21, 2024 · 5f23f1f · 5f23f1f
1 parent 417f22f
commit 5f23f1f
Showing 1 changed file with 23 additions and 7 deletions.
diff --git a/src/util/convertalignments.cpp b/src/util/convertalignments.cpp
@@ -349,6 +349,9 @@ int convertalignments(int argc, const char **argv, const Command &command) {
         std::string newBacktrace;
         newBacktrace.reserve(1024);
 
+        std::string complementBuffer;
+        complementBuffer.reserve(1024);
+
         const TaxonNode * taxonNode = NULL;
 
 #pragma omp  for schedule(dynamic, 10)
@@ -714,32 +717,45 @@ int convertalignments(int argc, const char **argv, const Command &command) {
                         break;
                     }
                     case Parameters::FORMAT_ALIGNMENT_SAM: {
-                        bool strand = res.qEndPos > res.qStartPos;
+                        bool forward = res.qEndPos > res.qStartPos;
                         int rawScore = static_cast<int>(evaluer->computeRawScoreFromBitScore(res.score) + 0.5);
                         uint32_t mapq = -4.343 * log(exp(static_cast<double>(-rawScore)));
                         mapq = (uint32_t) (mapq + 4.99);
                         mapq = mapq < 254 ? mapq : 254;
-                        int count = snprintf(buffer, sizeof(buffer), "%s\t%d\t%s\t%d\t%d\t",  queryId.c_str(), (strand) ? 16: 0, targetId.c_str(), res.dbStartPos + 1, mapq);
+                        int count = snprintf(buffer, sizeof(buffer), "%s\t%d\t%s\t%d\t%d\t",  queryId.c_str(), (forward) ? 0 : 16, targetId.c_str(), res.dbStartPos + 1, mapq);
                         if (count < 0 || static_cast<size_t>(count) >= sizeof(buffer)) {
                             Debug(Debug::WARNING) << "Truncated line in entry" << i << "!\n";
                             continue;
                         }
                         result.append(buffer, count);
-                        if (isTranslatedSearch == true && targetNucs == true && queryNucs == true) {
+                        if (isTranslatedSearch == true && queryNucs == true) {
                             Matcher::result_t::protein2nucl(res.backtrace, newBacktrace);
                             result.append(newBacktrace);
                             newBacktrace.clear();
-
                         } else {
                             result.append(res.backtrace);
                         }
                         result.append("\t*\t0\t0\t");
                         int start = std::min(res.qStartPos, res.qEndPos);
                         int end   = std::max(res.qStartPos, res.qEndPos);
-                        if (queryProfile) {
-                            result.append(queryProfData.c_str() + start, (end + 1) - start);
+                        if (forward == false && queryNucs == true) {
+                            if (queryProfile) {
+                                complementBuffer.assign(queryProfData.c_str() + start, (end + 1) - start);
+                            } else {
+                                complementBuffer.assign(querySeqData + start, (end + 1) - start);
+                            }
+                            std::reverse(complementBuffer.begin(), complementBuffer.end());
+                            for (size_t i = 0; i < complementBuffer.size(); i++) {
+                                complementBuffer[i] = Orf::complement(complementBuffer[i]);
+                            }
+                            result.append(complementBuffer);
+                            complementBuffer.clear();
                         } else {
-                            result.append(querySeqData + start, (end + 1) - start);
+                            if (queryProfile) {
+                                result.append(queryProfData.c_str() + start, (end + 1) - start);
+                            } else {
+                                result.append(querySeqData + start, (end + 1) - start);
+                            }
                         }
                         count = snprintf(buffer, sizeof(buffer), "\t*\tAS:i:%d\tNM:i:%d\n", rawScore, missMatchCount);
                         if (count < 0 || static_cast<size_t>(count) >= sizeof(buffer)) {