Changeset 8282

Timestamp:

08/12/11 11:31:56 (6 months ago)

Author:

eissler

Message:

PTPan: solved issue with including reverse-complement in ProbeMatch?

Location:

branches/ptpan_back/ptpan

Files:

• ptpan.h (modified) (1 diff)
• ptpan_legacy.cxx (modified) (3 diffs)
• ptpan_tree.cxx (modified) (8 diffs)
• ptpan_tree.h (modified) (1 diff)

branches/ptpan_back/ptpan/ptpan.h

@@ -463 +463 @@
     std::string sq_Query; /* query string */
 
-    bool sq_Reversed; /* query is reversed */
+    bool sq_Reversed; /* check reversed complement as well */
     bool sq_AllowReplace; /* allow replace operations */
     bool sq_AllowInsert; /* allow inserting operations */

branches/ptpan_back/ptpan/ptpan_legacy.cxx

@@ -381 +381 @@
 #endif
     gettimeofday(&tsNow, NULL);
-    printf("... done in %li seconds.\n",
-            (tsNow.tv_sec - tsStart.tv_sec));
+    printf("... done in %li seconds.\n", (tsNow.tv_sec - tsStart.tv_sec));
 #ifndef DEBUG
 }
@@ -1181 +1180 @@
     sq.sq_Query = std::string(locs->pm_sequence);
     sq.sq_MaxErrors = (float) locs->pm_max;
-    sq.sq_Reversed = false;
+    if (locs->pm_reversed) {
+        // include reversed-complement search as well
+        sq.sq_Reversed = true;
+    } else {
+        sq.sq_Reversed = false;
+    }
     if (sq.sq_MaxErrors > 0) {
         sq.sq_AllowReplace = true; // N-mismatches are listed anyways!
@@ -1201 +1205 @@
     ULONG done = create_match_list(sq, results, locs->pm_max_hits);
 
-    if (locs->pm_reversed) {
-        SearchQuery compsq = sq;
-        compsq.sq_Reversed = true;
-        results = pl->pl_pt->searchTree(compsq, true, true);
-        done = create_match_list(compsq, results, locs->pm_max_hits - done);
-    }
+#ifdef DEBUG
+    if (done) {
+        printf("done with search!\n");
+    }
+#endif
 
     // Just to be secure, we clear the PDC structure

branches/ptpan_back/ptpan/ptpan_tree.cxx

@@ -687 +687 @@
     if (!sq.sq_Query.empty()) {
         SearchQuery internal_sq = sq;
+        bool reversed_too = internal_sq.sq_Reversed;
         if (internal_sq.sq_MaxErrors == 0.0) {
             internal_sq.sq_AllowInsert = false;
@@ -692 +693 @@
             internal_sq.sq_AllowReplace = false;
         }
+        internal_sq.sq_Reversed = false;
         SearchQueryHandle sqh = SearchQueryHandle(internal_sq);
 
@@ -735 +737 @@
             }
         }
+
+        if (reversed_too) {
+            STRPTR query_str = const_cast<char *>(internal_sq.sq_Query.c_str());
+            m_as->complementSequence(query_str);
+            m_as->reverseSequence(query_str);
+            internal_sq.sq_Reversed = true;
+
+            SearchQueryHandle sqh_rev = SearchQueryHandle(internal_sq);
+#ifdef DEBUG
+            printf("check reverse complement as well (%s)\n",
+                    internal_sq.sq_Query.c_str());
+#endif
+            // TODO FIXME double code is bad, but there is no time left :-(
+            std::vector<PTPanPartition *>::const_iterator pit;
+            for (pit = m_partitions.begin(); pit != m_partitions.end(); pit++) {
+                search_partition(*pit, internal_sq, sqh_rev);
+            }
+            post_filter_query_hits(sq, sqh_rev);
+
+#ifdef DEBUG
+            printf("after all partitions rev-compl count %lu\n",
+                    sqh_rev.sqh_hits.size());
+#endif
+            if (m_prune_length < sqh_rev.sqh_MaxLength) {
+#ifdef DEBUG
+                printf("may contain UNSAFE! %lu\n", sqh_rev.sqh_MaxLength);
+#endif
+                validate_unsafe_hits(internal_sq, sqh_rev);
+            }
+            if (build_diff) {
+                if (!sqh_rev.sqh_hits.empty()) {
+#ifdef DEBUG
+                    printf(">> CreateDiffAlignment for %lu hits\n",
+                            sqh_rev.sqh_hits.size());
+#endif
+                    if (sqh_rev.sqh_hits.size() < m_num_threads
+                            || m_num_threads <= 1) {
+                        create_diff_alignment(internal_sq, sqh_rev);
+                    } else {
+                        std::size_t i = 0;
+                        ULONG participators = m_num_threads;
+                        for (; i < participators - 1; i++) {
+                            boost::threadpool::schedule(
+                                    *m_threadpool,
+                                    boost::bind(
+                                            &PtpanTree::create_diff_alignment,
+                                            this, boost::ref(internal_sq),
+                                            boost::ref(sqh_rev), i,
+                                            participators));
+                        }
+                        create_diff_alignment(internal_sq, sqh_rev, i,
+                                participators); // main thread can do some work as well!
+                        m_threadpool->wait();
+                    }
+                }
+            }
+            // we want to transfer control over results to sqh
+            // before we sort and return results
+            sqh.sqh_hits.transfer(sqh.sqh_hits.end(), sqh_rev.sqh_hits);
+        }
+
         if (sort) {
             sort_hits_list(internal_sq, sqh);
@@ -901 +964 @@
     UWORD seqcode = 0;
     ULONG length_corr =
-            (length < best_pp->pp_TreePruneLength) ? length :
-                    best_pp->pp_TreePruneLength;
+            (length < best_pp->pp_TreePruneLength) ?
+                    length : best_pp->pp_TreePruneLength;
 
     ULONG buffer_size;
@@ -2908 +2971 @@
         SearchQueryHandle& sqh) const {
     LONG query_len =
-            sq.sq_Query.size() > m_prune_length ? m_prune_length :
-                    (LONG) sq.sq_Query.size();
+            sq.sq_Query.size() > m_prune_length ?
+                    m_prune_length : (LONG) sq.sq_Query.size();
     LONG source_len =
-            sqh.sqh_MaxLength > m_prune_length ? m_prune_length :
-                    (LONG) sqh.sqh_MaxLength;
+            sqh.sqh_MaxLength > m_prune_length ?
+                    m_prune_length : (LONG) sqh.sqh_MaxLength;
 
     LONG max_check =
-            sq.sq_Query.size() < m_prune_length ? m_prune_length
-                    - sq.sq_Query.size() :
-                    0;
+            sq.sq_Query.size() < m_prune_length ?
+                    m_prune_length - sq.sq_Query.size() : 0;
     if (max_check > sq.sq_MaxErrors) {
         max_check = (LONG) sq.sq_MaxErrors;
@@ -3748 +3810 @@
                 qh->qh_sortKey =
                         (LLONG) (
-                                (qh->qh_Flags & QHF_REVERSED) ? (1LL << 62) :
-                                        0LL)
-                                + ((qh->qh_InsertCount | qh->qh_DeleteCount) ? (1LL
-                                        << 61) :
-                                        0LL)
+                                (qh->qh_Flags & QHF_REVERSED) ?
+                                        (1LL << 62) : 0LL)
+                                + ((qh->qh_InsertCount | qh->qh_DeleteCount) ?
+                                        (1LL << 61) : 0LL)
                                 + (((LLONG) (qh->qh_ReplaceCount
                                         + qh->qh_InsertCount
@@ -3780 +3841 @@
                 qh->qh_sortKey =
                         (LLONG) (
-                                (LLONG) (qh->qh_Flags & QHF_REVERSED) ? (1LL
-                                        << 62) :
-                                        0LL)
+                                (LLONG) (qh->qh_Flags & QHF_REVERSED) ?
+                                        (1LL << 62) : 0LL)
                                 + ((LLONG) (qh->qh_InsertCount
-                                        | qh->qh_DeleteCount) ? (1LL << 61) :
-                                        0LL)
+                                        | qh->qh_DeleteCount) ?
+                                        (1LL << 61) : 0LL)
                                 + (((LLONG) round(qh->qh_ErrorCount * 10.0))
                                         << 56)
@@ -4034 +4094 @@
     UWORD seqcode = m_as->wildcard_code + 1;
     ULONG length =
-            (dq.dq_ProbeLength < pp->pp_TreePruneLength) ? dq.dq_ProbeLength :
-                    pp->pp_TreePruneLength;
+            (dq.dq_ProbeLength < pp->pp_TreePruneLength) ?
+                    dq.dq_ProbeLength : pp->pp_TreePruneLength;
     double currtemp = 0.0;
     UWORD currgc = 0;

branches/ptpan_back/ptpan/ptpan_tree.h

@@ -136 +136 @@
             }
             sqh_hits.clear();
+        }
+
+        void clearHitsHash() {
+            if (sqh_HitsHash) {
+                sqh_HitsHash->clear();
+            }
         }