Merge branch 'master' of https://github.com/steineggerlab/foldseek

steineggerlab · Dec 26, 2023 · 7b68363 · 7b68363
2 parents 886021d + 592ffa8
commit 7b68363
Show file tree

Hide file tree

Showing 5 changed files with 26 additions and 47 deletions.
diff --git a/data/easycomplexsearch.sh b/data/easycomplexsearch.sh
@@ -26,30 +26,10 @@ if notExists "${TARGET}.dbtype"; then
     TARGET="${TMP_PATH}/target"
 fi
 
-if notExists "${TMP_PATH}/result.dbtype"; then
-    # shellcheck disable=SC2086
-    "$MMSEQS" search "${QUERY}" "${TARGET}" "${TMP_PATH}/result" "${TMP_PATH}/search_tmp" ${SEARCH_PAR} \
-        || fail "Search died"
-fi
-
-RESULT="${TMP_PATH}/result"
-if [ "$PREFMODE" != "EXHAUSTIVE" ]; then
-    if notExists "${TMP_PATH}/result_expand_pref.dbtype"; then
-        # shellcheck disable=SC2086
-        "$MMSEQS" expandcomplex "${QUERY}" "${TARGET}" "${RESULT}" "${TMP_PATH}/result_expand_pref" ${THREADS_PAR} \
-            || fail "Expandcomplex died"
-    fi
-    if notExists "${TMP_PATH}/result_expand_aligned.dbtype"; then
-        # shellcheck disable=SC2086
-        "$MMSEQS" $COMPLEX_ALIGNMENT_ALGO "${QUERY}" "${TARGET}" "${TMP_PATH}/result_expand_pref" "${TMP_PATH}/result_expand_aligned" ${COMPLEX_ALIGN_PAR} \
-            || fail "something died"
-    fi
-    RESULT="${TMP_PATH}/result_expand_aligned"
-fi
 if notExists "${TMP_PATH}/complex_result.dbtype"; then
     # shellcheck disable=SC2086
-    $MMSEQS scorecomplex "${QUERY}" "${TARGET}" "${RESULT}" "${TMP_PATH}/complex_result" ${SCORECOMPLEX_PAR} \
-        || fail "ScoreComplex died"
+    "$MMSEQS" complexsearch "${QUERY}" "${TARGET}" "${TMP_PATH}/complex_result" "${TMP_PATH}/complexsearch_tmp" ${COMPLEXSEARCH_PAR} \
+    || fail "ComplexSearch died"
 fi
 
 # shellcheck disable=SC2086
@@ -91,6 +71,6 @@ if [ -n "${REMOVE_TMP}" ]; then
         # shellcheck disable=SC2086
         "$MMSEQS" rmdb "${TMP_PATH}/query_ss" ${VERBOSITY}
     fi
-    rm -rf "${TMP_PATH}/search_tmp"
-    rm -f "${TMP_PATH}/easyscorecomplex.sh"
-fi
+    rm -rf "${TMP_PATH}/complexsearch_tmp"
+    rm -f "${TMP_PATH}/easycomplexsearch.sh"
+fi
diff --git a/src/FoldseekBase.cpp b/src/FoldseekBase.cpp
@@ -268,13 +268,13 @@ std::vector<Command> foldseekCommands = {
                 "# Search a single/multiple PDB file against a set of PDB files and get complex level alignments\n"
                 "foldseek complexsearch queryDB targetDB result tmp\n"
                 "# Format output differently\n"
-                "foldseek easy-complexsearch queryDB targetDB result tmp --format-output query,target,qstart,tstart,cigar\n"
+                "foldseek complexsearch queryDB targetDB result tmp --format-output query,target,qstart,tstart,cigar\n"
                 "# Align with TMalign (global)\n"
                 "foldseek complexsearch queryDB targetDB result tmp --alignment-type 1\n"
                 "# Skip prefilter and perform an exhaustive alignment (slower but more sensitive)\n"
                 "foldseek complexsearch queryDB targetDB result tmp --exhaustive-search 1\n\n",
                 "Woosub Kim <[email protected]>",
-                "<i:queryDB> <i:targetDB> <o:outputFileName> <tmpDir>",
+                "<i:queryDB> <i:targetDB> <o:alignmentDB> <tmpDir>",
                 CITATION_FOLDSEEK, {
                                            {"queryDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA | DbType::NEED_HEADER, &DbValidator::sequenceDb},
                                            {"targetDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA | DbType::NEED_HEADER, &DbValidator::sequenceDb},
@@ -302,7 +302,7 @@ std::vector<Command> foldseekCommands = {
                                    }
         },
         {"createcomplexreport", createcomplexreport, &localPar.createcomplexreport, COMMAND_FORMAT_CONVERSION,
-                "Convert complex DB to tsv format\"",
+                "Convert complexDB to tsv format",
                 "# Create output in tsv format (9 columns):  qComplexName.c_str(), tComplexName.c_str(), qChainString.c_str(), tChainString.c_str(), qTMScore, tTMScore, u, t, assId\n"
                 "#  (1,2) identifiers for query and target complex,\n"
                 "#  (3,4) chains of query complex and target complex,\n"
@@ -319,8 +319,8 @@ std::vector<Command> foldseekCommands = {
                                            {"complexFile", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::flatfile}
                                    }
         },
-        {"expandcomplex", expandcomplex, &localPar.expandcomplex, COMMAND_ALIGNMENT,
-        "expand complex",
+        {"expandcomplex", expandcomplex, &localPar.expandcomplex, COMMAND_PREFILTER,
+        "Re-prefilter to ensure complete alignment between complexes",
         NULL,
         "Woosub Kim <[email protected]>",
         "<i:queryDB> <i:targetDB> <i:alignmentDB> <o:prefilterDB>",

diff --git a/src/strucclustutils/expandcomplex.cpp b/src/strucclustutils/expandcomplex.cpp
@@ -6,7 +6,7 @@
 #include "LocalParameters.h"
 #include "MemoryMapped.h"
 #include "createcomplexreport.h"
-
+#include <set>
 #ifdef OPENMP
 #include <omp.h>
 #endif
@@ -53,7 +53,7 @@ int expandcomplex(int argc, const char **argv, const Command &command) {
         thread_idx = static_cast<unsigned int>(omp_get_thread_num());
 #endif
         resultToWrite_t result;
-        std::vector<unsigned int> dbFoundIndices;
+        std::set<unsigned int> dbFoundIndices;
         std::vector<ChainKeyPair_t>  chainKeyPairs;
 #pragma omp for schedule(dynamic, 1)
         // for each q complex
@@ -72,20 +72,23 @@ int expandcomplex(int argc, const char **argv, const Command &command) {
                     const auto dbChainKey = (unsigned int) strtoul(dbKeyBuffer, NULL, 10);
                     const unsigned int dbComplexId = dbChainKeyToComplexIdMap.at(dbChainKey);
                     // find all db complex aligned to the query complex.
-                    if (std::find(dbFoundIndices.begin(), dbFoundIndices.end(), dbComplexId) == dbFoundIndices.end())
-                        dbFoundIndices.emplace_back(dbComplexId);
+                    dbFoundIndices.insert(dbComplexId);
                     data = Util::skipLine(data);
                 }
             }
-            if (dbFoundIndices.empty())
+            if (dbFoundIndices.empty()) {
+                for (size_t qChainIdx=0; qChainIdx<qChainKeys.size(); qChainIdx++) {
+                    resultWriter.writeData(result.c_str(),result.length(),qChainKeys[qChainIdx],thread_idx);
+                }
                 continue;
+            }
             // Among all db complexes aligned to query complex
-            for (size_t dbIdx=0; dbIdx<dbFoundIndices.size(); dbIdx++) {
-                std::vector<unsigned int> &dbChainKeys = dbComplexIdToChainKeysMap.at(dbFoundIndices[dbIdx]);
+            for (auto dbIter = dbFoundIndices.begin(); dbIter != dbFoundIndices.end(); dbIter++) {
+                std::vector<unsigned int> &dbChainKeys = dbComplexIdToChainKeysMap.at(*dbIter);
                 // for all query chains
                 for (size_t qChainIdx=0; qChainIdx<qChainKeys.size(); qChainIdx++) {
                     // and target chains
-                    for (size_t dbChainIdx=0; dbChainIdx<dbChainKeys.size(); dbChainIdx++) {
+                    for (size_t dbChainIdx = 0; dbChainIdx < dbChainKeys.size(); dbChainIdx++) {
                         // get all possible alignments
                         chainKeyPairs.emplace_back(qChainKeys[qChainIdx], dbChainKeys[dbChainIdx]);
                     }

diff --git a/src/strucclustutils/scorecomplex.cpp b/src/strucclustutils/scorecomplex.cpp
@@ -342,15 +342,12 @@ class DBSCANCluster {
             }
         }
 
-//        if (maxClusterSize == idealClusterSize) {
-//            bestClusters = currClusters;
-//            prevMaxClusterSize = maxClusterSize;
-//            return finishDBSCAN();
-//        } else
         if (maxClusterSize < prevMaxClusterSize)
             return finishDBSCAN();
-        else if (maxClusterSize == prevMaxClusterSize && currClusters.size() < bestClusters.size())
+
+        if (maxClusterSize == prevMaxClusterSize && currClusters.size() < bestClusters.size())
             return finishDBSCAN();
+
         bestClusters = currClusters;
         prevMaxClusterSize = maxClusterSize;
         eps += LEARNING_RATE;

diff --git a/src/workflow/EasyComplexSearch.cpp b/src/workflow/EasyComplexSearch.cpp
@@ -118,8 +118,7 @@ int easycomplexsearch(int argc, const char **argv, const Command &command) {
     cmd.addVariable("LEAVE_INPUT", par.dbOut ? "TRUE" : NULL);
     par.filenames.pop_back();
     cmd.addVariable("CREATEDB_PAR", par.createParameterString(par.structurecreatedb).c_str());
-    cmd.addVariable("SEARCH_PAR", par.createParameterString(par.structuresearchworkflow, true).c_str());
-    cmd.addVariable("SCORECOMPLEX_PAR", par.createParameterString(par.scorecomplex).c_str());
+    cmd.addVariable("COMPLEXSEARCH_PAR", par.createParameterString(par.complexsearchworkflow).c_str());
     cmd.addVariable("CONVERT_PAR", par.createParameterString(par.convertalignments).c_str());
     cmd.addVariable("REPORT_PAR", par.createParameterString(par.createcomplexreport).c_str());
     cmd.addVariable("THREADS_PAR", par.createParameterString(par.onlythreads).c_str());
@@ -131,4 +130,4 @@ int easycomplexsearch(int argc, const char **argv, const Command &command) {
     // Should never get here
     assert(false);
     return EXIT_FAILURE;
-}
+}