alignproteome added

data/workflow/CMakeLists.txt

@@ -3,6 +3,7 @@ set(GENERATED_WORKFLOWS
         workflow/easycluster.sh
         workflow/easytaxonomy.sh
         workflow/easyrbh.sh
+        workflow/easyalignproteome.sh
         workflow/blastp.sh
         workflow/blastpgp.sh
         workflow/map.sh

data/workflow/easyalignproteome.sh

@@ -0,0 +1,53 @@
+#!/bin/sh -e
+fail() {
+    echo "Error: $1"
+    exit 1
+}
+
+notExists() {
+   [ ! -f "$1" ]
+}
+
+
+if notExists "${TMP_PATH}/input.dbtype"; then
+    # shellcheck disable=SC2086
+    "$MMSEQS" createdb "$@" "${TMP_PATH}/input" ${CREATEDB_PAR} \
+        || fail "query createdb died"
+fi
+
+if notExists "${TMP_PATH}/clu.dbtype"; then
+    # shellcheck disable=SC2086
+    "$MMSEQS" linclust "${TMP_PATH}/input" "${TMP_PATH}/clu" "${TMP_PATH}/clu_tmp" ${CLUSTER_PAR} \
+        || fail "Search died"
+fi
+
+if notExists "${TMP_PATH}/aln.dbtype"; then
+    # shellcheck disable=SC2086
+    "$MMSEQS" alignproteome "${TMP_PATH}/input" "${TMP_PATH}/clu" "${TMP_PATH}/aln_protein" "${TMP_PATH}/aln_proteome" ${ALIGNPROTEOME_PAR} \
+        || fail "Convert Alignments died"
+fi
+
+if notExists "${RESULTS}_protein_cluster.tsv"; then
+    # shellcheck disable=SC2086
+    "$MMSEQS" createtsv "${TMP_PATH}/input" "${TMP_PATH}/clu" "${RESULTS}_protein_cluster.tsv" ${THREADS_PAR} \
+            || fail "createtsv protein died"
+fi
+
+if notExists "${RESULTS}_proteome_cluster.tsv"; then
+    # shellcheck disable=SC2086
+    "$MMSEQS" createtsv "${TMP_PATH}/input" "${TMP_PATH}/aln_proteome" "${RESULTS}_proteome_cluster.tsv" ${THREADS_PAR} \
+            || fail "createtsv proteome died"
+fi
+
+if [ -n "${REMOVE_TMP}" ]; then
+    # shellcheck disable=SC2086
+    "$MMSEQS" rmdb "${TMP_PATH}/input" ${VERBOSITY_PAR}
+    # shellcheck disable=SC2086
+    "$MMSEQS" rmdb "${TMP_PATH}/input_h" ${VERBOSITY_PAR}
+    # shellcheck disable=SC2086
+    "$MMSEQS" rmdb "${TMP_PATH}/clu" ${VERBOSITY_PAR}
+    # shellcheck disable=SC2086
+    "$MMSEQS" rmdb "${TMP_PATH}/aln" ${VERBOSITY_PAR}
+    rm -rf "${TMP_PATH}/clu_tmp"
+    rm -f "${TMP_PATH}/easyproteomecluster.sh"
+fi

src/CommandDeclarations.h

@@ -44,6 +44,7 @@ extern int diffseqdbs(int argc, const char **argv, const Command& command);
 extern int easycluster(int argc, const char **argv, const Command& command);
 extern int easyrbh(int argc, const char **argv, const Command& command);
 extern int easylinclust(int argc, const char **argv, const Command& command);
+extern int easyalignproteome(int argc, const char **argv, const Command& command);
 extern int easysearch(int argc, const char **argv, const Command& command);
 extern int easylinsearch(int argc, const char **argv, const Command& command);
 extern int tsv2exprofiledb(int argc, const char **argv, const Command& command);

src/MMseqsBase.cpp

@@ -79,6 +79,12 @@ std::vector<Command> baseCommands = {
                 CITATION_MMSEQS2|CITATION_LINCLUST, {{"fastaFile[.gz|.bz2]", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA|DbType::VARIADIC, &DbValidator::flatfileAndStdin },
                                                             {"clusterPrefix", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::flatfile },
                                                             {"tmpDir", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::directory }}},
+        {"easy-alignproteome", easyalignproteome, &par.easyalignproteome, COMMAND_EASY,
+                "Proteome clustering and alignemnt",
+                "<i:fastaFile1[.gz|.bz2]> ... <i:fastaFileN[.gz|.bz2]>"
+                // CITATION_MMSEQS2,
+
+        },
         {"easy-taxonomy",        easytaxonomy,         &par.easytaxonomy,         COMMAND_EASY,
                 "Taxonomic classification",
                 "# Assign taxonomic labels to FASTA sequences\n"
@@ -633,9 +639,9 @@ std::vector<Command> baseCommands = {
                                                            {"resultDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::resultDb },
                                                            {"alignmentDB", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::alignmentDb }}},
         {"alignproteome",             alignproteome,             &par.alignproteome,             COMMAND_ALIGNMENT,
-                "Within-result all-vs-all gapped local alignment",
+                "Proteome clustering and alignment",
                 NULL,
-                "Martin Steinegger <[email protected]>",
+                "Martin Steinegger <[email protected]> & Gyuri Kim <[email protected]>",
                 "<i:sequenceDB> <i:resultDB> <o:alignmentDB>",
                 CITATION_MMSEQS2, {{"sequenceDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::sequenceDb },
                                                            {"resultDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::resultDb },

src/alignment/Matcher.h

@@ -221,7 +221,6 @@ class Matcher{
 
 
     static size_t resultToBuffer(char * buffer, const result_t &result, bool addBacktrace, bool compress  = true, bool addOrfPosition = false);
-    static size_t resultToBuffer_str(char * buffer, const result_t &result, bool addBacktrace, bool compress, bool addOrfPosition=false);
     static int computeAlnLength(int anEnd, int start, int dbEnd, int dbStart);
 
     static void updateResultByRescoringBacktrace(const char *querySeq, const char *targetSeq, const char **subMat, EvalueComputation &evaluer,

src/commons/CMakeLists.txt

@@ -8,7 +8,6 @@ set(commons_header_files
         commons/Command.h
         commons/CommandCaller.h
         commons/Concat.h
-        # commons/ClusterSpecies.h
         commons/DBConcat.h
         commons/DBReader.h
         commons/DBWriter.h
@@ -53,7 +52,6 @@ set(commons_source_files
         commons/BaseMatrix.cpp
         commons/Command.cpp
         commons/CommandCaller.cpp
-        # commons/ClusterSpecies.cpp
         commons/DBConcat.cpp
         commons/DBReader.cpp
         commons/DBWriter.cpp

src/commons/DBReader.cpp

@@ -137,8 +137,6 @@ template <typename T> bool DBReader<T>::open(int accessType){
         }
     }
     if (dataMode & USE_LOOKUP || dataMode & USE_LOOKUP_REV) {
-        Debug(Debug::INFO) << "ReadLookup file: " << dataFileName << "\n"; //gyuri
-
         std::string lookupFilename = (std::string(dataFileName) + ".lookup");
         MemoryMapped lookupData(lookupFilename, MemoryMapped::WholeFile, MemoryMapped::SequentialScan);
         if (lookupData.isValid() == false) {
@@ -147,9 +145,7 @@ template <typename T> bool DBReader<T>::open(int accessType){
         }
         char* lookupDataChar = (char *) lookupData.getData();
         size_t lookupDataSize = lookupData.size();
-        Debug(Debug::INFO) << "Lookup Data size is " << lookupDataSize << "\n"; //gyuri
         lookupSize = Util::ompCountLines(lookupDataChar, lookupDataSize, threads);
-        Debug(Debug::INFO) << "Lookup size is " << lookupSize << "\n"; //gyuri
         lookup = new(std::nothrow) LookupEntry[this->lookupSize];
         incrementMemory(sizeof(LookupEntry) * this->lookupSize);
         readLookup(lookupDataChar, lookupDataSize, lookup);
@@ -1025,19 +1021,6 @@ size_t DBReader<T>::getOffset(size_t id) {
     return index[id].offset;
 }
 
-template<typename T>
-size_t DBReader<T>::getIndexLen(size_t id) {
-    if (id >= size){
-        Debug(Debug::ERROR) << "Invalid database read for id=" << id << ", database index=" << indexFileName << "\n";
-        Debug(Debug::ERROR) << "getOffset: local id (" << id << ") >= db size (" << size << ")\n";
-        EXIT(EXIT_FAILURE);
-    }
-    if (local2id != NULL) {
-        id = local2id[id];
-    }
-    return index[id].length;
-}
-
 template<typename T>
 size_t DBReader<T>::findNextOffsetid(size_t id) {
     size_t idOffset = getOffset(id);

src/commons/DBReader.h

@@ -187,7 +187,7 @@ class DBReader : public MemoryTracker {
 
     size_t getSize() const;
 
-    unsigned int getProteomeTotalLen(size_t id); //gyuri
+    unsigned int getProteomeTotalLen(size_t id); 
 
     unsigned int getMaxSeqLen(){ 
             return (Parameters::isEqualDbtype(dbtype, Parameters::DBTYPE_HMM_PROFILE ) ) ?
@@ -446,8 +446,6 @@ class DBReader : public MemoryTracker {
 
     size_t findNextOffsetid(size_t id);
 
-    size_t getIndexLen(size_t id);
-
     int isCompressed(){
         return isCompressed(dbtype);
     }

src/commons/DBWriter.cpp

@@ -283,7 +283,6 @@ size_t DBWriter::writeAdd(const char* data, size_t dataSize, unsigned int thrIdx
     }
     size_t totalWriten = 0;
     if(isCompressedDB && (state[thrIdx] == INIT_STATE || state[thrIdx] == COMPRESSED) ) {
-        Debug(Debug::INFO) << "Compressing data in thread " << thrIdx << "\n"; //gyuri
         state[thrIdx] = COMPRESSED;
         // zstd seems to have a hard time with elements < 60
         ZSTD_inBuffer input = { data, dataSize, 0 };

src/commons/Parameters.cpp

@@ -1349,6 +1349,10 @@ Parameters::Parameters():
     // easylinclustworkflow
     easylinclustworkflow = combineList(linclustworkflow, createdb);
 
+    // easyalignproteomeworkflow
+    easyalignproteome = combineList(easylinclustworkflow, alignproteome);
+
+
     // clustering workflow
     clusterworkflow = combineList(prefilter, align);
     clusterworkflow = combineList(clusterworkflow, rescorediagonal);
@@ -2353,7 +2357,7 @@ void Parameters::setDefaults() {
     maxRejected = INT_MAX;
     maxAccept   = INT_MAX;
     seqIdThr = 0.0;
-    proteomeSimThr = 0.0;
+    proteomeSimThr = 0.9;
     alnLenThr = 0;
     altAlignment = 0;
     gapOpen = MultiParam<NuclAA<int>>(NuclAA<int>(11, 5));

src/commons/Parameters.h

@@ -1118,6 +1118,7 @@ class Parameters {
     std::vector<MMseqsParameter*> kmersearch;
     std::vector<MMseqsParameter*> countkmer;
     std::vector<MMseqsParameter*> easylinclustworkflow;
+    std::vector<MMseqsParameter*> easyalignproteome;
     std::vector<MMseqsParameter*> linclustworkflow;
     std::vector<MMseqsParameter*> easysearchworkflow;
     std::vector<MMseqsParameter*> searchworkflow;

src/util/alignall.cpp

@@ -27,7 +27,6 @@ int alignall(int argc, const char **argv, const Command &command) {
     }
     unsigned int swMode = Alignment::initSWMode(par.alignmentMode, par.covThr, par.seqIdThr);
 
-    // DBReader<unsigned int> tdbr(par.db1.c_str(), par.db1Index.c_str(), par.threads, DBReader<unsigned int>::USE_DATA|DBReader<unsigned int>::USE_INDEX);
     DBReader<unsigned int> tdbr(par.db1.c_str(), par.db1Index.c_str(), par.threads, DBReader<unsigned int>::USE_DATA|DBReader<unsigned int>::USE_INDEX|DBReader<unsigned int>::USE_LOOKUP);
 
     tdbr.open(DBReader<unsigned int>::NOSORT);
@@ -58,7 +57,6 @@ int alignall(int argc, const char **argv, const Command &command) {
     EvalueComputation evaluer(tdbr.getAminoAcidDBSize(), subMat, gapOpen, gapExtend);
     const size_t flushSize = 100000000;
     size_t iterations = static_cast<int>(ceil(static_cast<double>(dbr_res.getSize()) / static_cast<double>(flushSize)));
-    Debug(Debug::INFO) << "Number of iterations: " << iterations << " Size of linclust dbr_res : " <<  dbr_res.getSize() <<  '\n';
 
     for (size_t i = 0; i < iterations; i++) {
         size_t start = (i * flushSize);

src/util/alignproteome.cpp

@@ -50,7 +50,7 @@ struct __attribute__((__packed__)) ProteomeEntry{
         totalSeqId += seqId;
     }
 
-    void addSeqLength(unsigned int seqLength) {
+    void addSeqLen(unsigned int seqLength) {
         proteomeAALen += seqLength;
     }
 
@@ -64,14 +64,14 @@ struct __attribute__((__packed__)) ProteomeEntry{
     }
 };
 
-struct __attribute__((__packed__)) memberProteinEntry{
+struct __attribute__((__packed__)) MemberProteinEntry{
     unsigned int proteomeKey;
     unsigned int proteinKey;
 };
 
 struct ClusterEntry {
     bool isAvailable;
-    std::vector<memberProteinEntry> memberProteins;
+    std::vector<MemberProteinEntry> memberProteins;
 
     ClusterEntry() : isAvailable(false) {}
 
@@ -89,7 +89,7 @@ void calculateProteomeLength(std::vector<ProteomeEntry>& ProteomeList, DBReader<
     for (size_t i = 0; i < lookupSize; i++) {
         const unsigned int ProteomeId = lookup[i].fileNumber;
         const unsigned int ProteinId = lookup[i].id;
-        ProteomeList[ProteomeId].addSeqLength(tProteinDB.getIndexLen(ProteinId) - 2);
+        ProteomeList[ProteomeId].addSeqLen(tProteinDB.getSeqLen(ProteinId));
     }
 }
 
@@ -110,7 +110,7 @@ void initLocalClusterReps(size_t& id, std::vector<ClusterEntry>& localClusterRep
         for (auto& eachMemberKey : memberKeys){
             const unsigned int proteinId = tProteinDB.getId(eachMemberKey);
             const unsigned int proteomeKey = tProteinDB.getLookupFileNumber(proteinId);
-            memberProteinEntry mem;
+            MemberProteinEntry mem;
             mem.proteomeKey = proteomeKey;
             mem.proteinKey = proteinId;
             eachClusterRep.memberProteins.push_back(mem);

src/workflow/CMakeLists.txt

@@ -14,6 +14,7 @@ set(workflow_source_files
         workflow/Search.cpp
         workflow/Taxonomy.cpp
         workflow/EasyTaxonomy.cpp
+        workflow/EasyAlignproteome.cpp
         workflow/CreateIndex.cpp
         PARENT_SCOPE
         )

src/workflow/EasyAlignproteome.cpp

@@ -0,0 +1,41 @@
+#include <cassert>
+
+#include "FileUtil.h"
+#include "CommandCaller.h"
+#include "Util.h"
+#include "Debug.h"
+#include "Parameters.h"
+
+#include "easyalignproteome.sh.h"
+
+void setEasyAlignproteomeDefaults(Parameters *p) {
+    p->proteomeSimThr = 0.9;
+}
+
+// void setEasyAlignproteomeMustPassAlong(Parameters *p){
+
+// }
+
+int easyalignproteome(int argc, const char **argv, const Command &command) {
+    Parameters &par = Parameters::getInstance();
+    setEasyAlignproteomeDefaults(&par);
+    par.parseParameters(argc, argv, command, true, Parameters::PARSE_VARIADIC, 0);
+    // setEasyAlignproteomeMustPassAlong(&par);
+    std::string tmpDir = par.filenames.back();
+    std::string hash = SSTR(par.hashParameter(command.databases, par.filenames, *command.params));
+    if (par.reuseLatest) {
+        hash = FileUtil::getHashFromSymLink(tmpDir + "/latest");
+    }
+    tmpDir = FileUtil::createTemporaryDirectory(tmpDir, hash);
+    par.filenames.pop_back();
+
+    CommandCaller cmd;
+    cmd.addVariable("ALIGNPROTEOME_PAR", par.createParameterString(par.alignproteome,true).c_str()); // what?
+    std::string program = tmpDir + "/easyalignproteome.sh";
+    FileUtil::writeFile(program, easyalignproteome_sh, easyalignproteome_sh_len);
+    cmd.execProgram(program.c_str(), par.filenames);
+
+    // Should never get here
+    assert(false);
+    return 0;
+}