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

@@ -4,6 +4,7 @@
 
 extern int align(int argc, const char **argv, const Command& command);
 extern int alignall(int argc, const char **argv, const Command& command);
+extern int alignproteome(int argc, const char **argv, const Command& command);
 extern int alignbykmer(int argc, const char **argv, const Command& command);
 extern int appenddbtoindex(int argc, const char **argv, const Command& command);
 extern int apply(int argc, const char **argv, const Command& command);
@@ -43,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"
@@ -632,6 +638,15 @@ std::vector<Command> baseCommands = {
                                                            {"targetDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::sequenceDb },
                                                            {"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,
+                "Proteome clustering and alignment",
+                NULL,
+                "Martin Steinegger <martin.steinegger@snu.ac.kr> & Gyuri Kim <gyuribio@snu.ac.kr>",
+                "<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 },
+                                                           {"alignmentDB", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::alignmentDb },
+                                                           {"alignmentDB", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::alignmentDb }}},
         {"alignall",             alignall,             &par.alignall,             COMMAND_ALIGNMENT,
                 "Within-result all-vs-all gapped local alignment",
                 NULL,
@@ -640,6 +655,7 @@ std::vector<Command> baseCommands = {
                 CITATION_MMSEQS2, {{"sequenceDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::sequenceDb },
                                                            {"resultDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::resultDb },
                                                            {"alignmentDB", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::alignmentDb }}},
+
         {"transitivealign",      transitivealign,      &par.align,                COMMAND_ALIGNMENT,
                 "Transfer alignments via transitivity",
                 //"Infer the alignment A->C via B, B being the center sequence and A,C each pairwise aligned against B",

src/alignment/Matcher.h

@@ -140,6 +140,10 @@ class Matcher{
                 }
             }
         }
+
+        float getSeqId(){
+            return seqId;
+        }
     };
 
     Matcher(int querySeqType, int targetSeqType, int maxSeqLen, BaseMatrix *m,
@@ -217,7 +221,6 @@ class Matcher{
 
 
     static size_t resultToBuffer(char * buffer, const result_t &result, bool addBacktrace, bool compress  = true, 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/DBReader.cpp

@@ -112,6 +112,7 @@ template <typename T> bool DBReader<T>::open(int accessType){
         }
         totalDataSize = 0;
         dataFileCnt = dataFileNames.size();
+
         dataSizeOffset = new size_t[dataFileNames.size() + 1];
         dataFiles = new char*[dataFileNames.size()];
         for(size_t fileIdx = 0; fileIdx < dataFileNames.size(); fileIdx++){
@@ -155,6 +156,22 @@ template <typename T> bool DBReader<T>::open(int accessType){
         }
         lookupData.close();
     }
+    if (dataMode & USE_SOURCE){
+        std::string sourceFilename = (std::string(dataFileName) + ".source");
+        MemoryMapped sourceData(sourceFilename, MemoryMapped::WholeFile, MemoryMapped::SequentialScan);
+        if (sourceData.isValid() == false) {
+            Debug(Debug::ERROR) << "Cannot open source file " << sourceFilename << "!\n";
+            EXIT(EXIT_FAILURE);
+        }
+        char* sourceDataChar = (char *) sourceData.getData();
+        size_t sourceDataSize = sourceData.size();
+        sourceSize = Util::ompCountLines(sourceDataChar, sourceDataSize, threads);
+        source = new(std::nothrow) SourceEntry[this->sourceSize];
+        incrementMemory(sizeof(SourceEntry) * this->sourceSize);
+        readSource(sourceDataChar, sourceDataSize, source);
+
+    }
+
     bool isSortedById = false;
     if (externalData == false) {
         MemoryMapped indexData(indexFileName, MemoryMapped::WholeFile, MemoryMapped::SequentialScan);
@@ -707,6 +724,15 @@ template <typename T> unsigned int DBReader<T>::getLookupFileNumber(size_t id){
     return lookup[id].fileNumber;
 }
 
+template <typename T> std::string DBReader<T>::getSourceFileName (size_t id){
+    if (id >= sourceSize){
+        Debug(Debug::ERROR) << "Invalid database read for id=" << id << ", database index=" << dataFileName << ".source\n";
+        Debug(Debug::ERROR) << "getSourceFileName: local id (" << id << ") >= db size (" << sourceSize << ")\n";
+        EXIT(EXIT_FAILURE);
+    }
+    return source[id].fileName;
+}
+
 template<>
 void DBReader<unsigned int>::lookupEntryToBuffer(std::string& buffer, const LookupEntry& entry) {
     buffer.append(SSTR(entry.id));
@@ -1048,6 +1074,36 @@ void DBReader<T>::readLookup(char *data, size_t dataSize, DBReader::LookupEntry
         currPos = lookupData - (char *) data;
 
         i++;
+
+
+    }
+}
+
+template <typename T>
+void DBReader<T>::readSource(char *data, size_t dataSize, DBReader::SourceEntry *source) {
+    size_t i=0;
+    size_t currPos = 0;
+    char* sourceData = (char *) data;
+    const char * cols[3];
+    while (currPos < dataSize){
+        if (i >= this->sourceSize) {
+            Debug(Debug::ERROR) << "Corrupt memory, too many entries!\n";
+            EXIT(EXIT_FAILURE);
+        }
+        Util::getFieldsOfLine(sourceData, cols, 3);
+        source[i].id = Util::fast_atoi<size_t>(cols[0]);
+        std::string fileName = std::string(cols[1], (cols[2] - cols[1]));
+        size_t lastDotPosition = fileName.rfind('.');
+
+        if (lastDotPosition != std::string::npos) {
+            fileName = fileName.substr(0, lastDotPosition);
+        }
+        source[i].fileName = fileName;
+        sourceData = Util::skipLine(sourceData);
+
+        currPos = sourceData - (char *) data;
+
+        i++;
     }
 }
 

src/commons/DBReader.h

@@ -98,7 +98,10 @@ class DBReader : public MemoryTracker {
             return false;
         }
     };
-
+    struct SourceEntry{
+        T id;
+        std::string fileName;
+    };
     struct LookupEntry {
         T id;
         std::string entryName;
@@ -184,6 +187,8 @@ class DBReader : public MemoryTracker {
 
     size_t getSize() const;
 
+    unsigned int getProteomeTotalLen(size_t id); 
+
     unsigned int getMaxSeqLen(){ 
             return (Parameters::isEqualDbtype(dbtype, Parameters::DBTYPE_HMM_PROFILE ) ) ?
                     (std::max(maxSeqLen, 1u)) / Sequence::PROFILE_READIN_SIZE :
@@ -248,6 +253,7 @@ class DBReader : public MemoryTracker {
     void lookupEntryToBuffer(std::string& buffer, const LookupEntry& entry);
     LookupEntry* getLookup() { return lookup; };
 
+    std::string getSourceFileName(size_t id);
     static const int NOSORT = 0;
     static const int SORT_BY_LENGTH = 1;
     static const int LINEAR_ACCCESS = 2;
@@ -266,6 +272,7 @@ class DBReader : public MemoryTracker {
     static const unsigned int USE_FREAD      = 4;
     static const unsigned int USE_LOOKUP     = 8;
     static const unsigned int USE_LOOKUP_REV = 16;
+    static const unsigned int USE_SOURCE     = 32;
 
 
     // compressed
@@ -309,6 +316,8 @@ class DBReader : public MemoryTracker {
 
     void readLookup(char *data, size_t dataSize, LookupEntry *lookup);
 
+    void readSource(char *data, size_t dataSize, SourceEntry *source);
+
     void readIndexId(T* id, char * line, const char** cols);
 
     unsigned int indexIdToNum(T* id);
@@ -485,7 +494,9 @@ class DBReader : public MemoryTracker {
 
     Index * index;
     size_t lookupSize;
+    size_t sourceSize;
     LookupEntry * lookup;
+    SourceEntry * source;
     bool sortedByOffset;
 
     unsigned int * id2local;

src/commons/Parameters.cpp

@@ -86,6 +86,8 @@ Parameters::Parameters():
         PARAM_ALT_ALIGNMENT(PARAM_ALT_ALIGNMENT_ID, "--alt-ali", "Alternative alignments", "Show up to this many alternative alignments", typeid(int), (void *) &altAlignment, "^[0-9]{1}[0-9]*$", MMseqsParameter::COMMAND_ALIGN),
         PARAM_GAP_OPEN(PARAM_GAP_OPEN_ID, "--gap-open", "Gap open cost", "Gap open cost", typeid(MultiParam<NuclAA<int>>), (void *) &gapOpen, "^[0-9]{1}[0-9]*$", MMseqsParameter::COMMAND_ALIGN | MMseqsParameter::COMMAND_EXPERT),
         PARAM_GAP_EXTEND(PARAM_GAP_EXTEND_ID, "--gap-extend", "Gap extension cost", "Gap extension cost", typeid(MultiParam<NuclAA<int>>), (void *) &gapExtend, "^[0-9]{1}[0-9]*$", MMseqsParameter::COMMAND_ALIGN | MMseqsParameter::COMMAND_EXPERT),
+        PARAM_PROTEOME_SIMILARITY(PARAM_PROTEOME_SIMILARITY_ID, "--proteome-similarity", "Proteome similarity", "Proteome similarity threshold", typeid(float), (void *) &proteomeSimThr, "^0(\\.[0-9]+)?|1(\\.0+)?$", MMseqsParameter::COMMAND_ALIGN),
+
 #ifdef GAP_POS_SCORING
         PARAM_GAP_PSEUDOCOUNT(PARAM_GAP_PSEUDOCOUNT_ID, "--gap-pc", "Gap pseudo count", "Pseudo count for calculating position-specific gap opening penalties", typeid(int), &gapPseudoCount, "^[0-9]+$", MMseqsParameter::COMMAND_ALIGN|MMseqsParameter::COMMAND_EXPERT),
 #endif
@@ -338,6 +340,7 @@ Parameters::Parameters():
     alignall.push_back(&PARAM_SUB_MAT);
     alignall.push_back(&PARAM_ADD_BACKTRACE);
     alignall.push_back(&PARAM_ALIGNMENT_MODE);
+    alignall.push_back(&PARAM_REALIGN_SCORE_BIAS); //gyuri
 //    alignall.push_back(&PARAM_WRAPPED_SCORING);
     alignall.push_back(&PARAM_E);
     alignall.push_back(&PARAM_MIN_SEQ_ID);
@@ -350,7 +353,7 @@ Parameters::Parameters():
     alignall.push_back(&PARAM_NO_COMP_BIAS_CORR);
     alignall.push_back(&PARAM_NO_COMP_BIAS_CORR_SCALE);
 
-//    alignall.push_back(&PARAM_REALIGN);
+    alignall.push_back(&PARAM_REALIGN);
 //    alignall.push_back(&PARAM_MAX_REJECTED);
 //    alignall.push_back(&PARAM_MAX_ACCEPT);
     alignall.push_back(&PARAM_INCLUDE_IDENTITY);
@@ -365,6 +368,41 @@ Parameters::Parameters():
     alignall.push_back(&PARAM_COMPRESSED);
     alignall.push_back(&PARAM_V);
 
+    //alignproteome
+        // alignproteome
+    alignproteome.push_back(&PARAM_SUB_MAT);
+    alignproteome.push_back(&PARAM_ADD_BACKTRACE);
+    alignproteome.push_back(&PARAM_ALIGNMENT_MODE);
+    alignproteome.push_back(&PARAM_REALIGN_SCORE_BIAS); //gyuri
+//    alignproteome.push_back(&PARAM_WRAPPED_SCORING);
+    alignproteome.push_back(&PARAM_E);
+    alignproteome.push_back(&PARAM_MIN_SEQ_ID);
+    alignproteome.push_back(&PARAM_MIN_ALN_LEN);
+    alignproteome.push_back(&PARAM_SEQ_ID_MODE);
+//    alignproteome.push_back(&PARAM_ALT_ALIGNMENT);
+    alignproteome.push_back(&PARAM_C);
+    alignproteome.push_back(&PARAM_COV_MODE);
+    alignproteome.push_back(&PARAM_MAX_SEQ_LEN);
+    alignproteome.push_back(&PARAM_NO_COMP_BIAS_CORR);
+    alignproteome.push_back(&PARAM_NO_COMP_BIAS_CORR_SCALE);
+    alignproteome.push_back(&PARAM_PROTEOME_SIMILARITY);
+
+    alignproteome.push_back(&PARAM_REALIGN);
+//    alignproteome.push_back(&PARAM_MAX_REJECTED);
+//    alignproteome.push_back(&PARAM_MAX_ACCEPT);
+    alignproteome.push_back(&PARAM_INCLUDE_IDENTITY);
+    alignproteome.push_back(&PARAM_PRELOAD_MODE);
+    alignproteome.push_back(&PARAM_PCA);
+    alignproteome.push_back(&PARAM_PCB);
+    alignproteome.push_back(&PARAM_SCORE_BIAS);
+    alignproteome.push_back(&PARAM_GAP_OPEN);
+    alignproteome.push_back(&PARAM_GAP_EXTEND);
+    alignproteome.push_back(&PARAM_ZDROP);
+    alignproteome.push_back(&PARAM_THREADS);
+    alignproteome.push_back(&PARAM_COMPRESSED);
+    alignproteome.push_back(&PARAM_V);
+
+
     // alignment
     align.push_back(&PARAM_SUB_MAT);
     align.push_back(&PARAM_ADD_BACKTRACE);
@@ -1311,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);
@@ -2315,6 +2357,7 @@ void Parameters::setDefaults() {
     maxRejected = INT_MAX;
     maxAccept   = INT_MAX;
     seqIdThr = 0.0;
+    proteomeSimThr = 0.9;
     alnLenThr = 0;
     altAlignment = 0;
     gapOpen = MultiParam<NuclAA<int>>(NuclAA<int>(11, 5));

src/commons/Parameters.h

@@ -434,6 +434,7 @@ class Parameters {
     int    maxAccept;                    // after n accepted sequences stop
     int    altAlignment;                 // show up to this many alternative alignments
     float  seqIdThr;                     // sequence identity threshold for acceptance
+    float  proteomeSimThr;
     int    alnLenThr;                    // min. alignment length
     bool   addBacktrace;                 // store backtrace string (M=Match, D=deletion, I=insertion)
     bool   realign;                      // realign hit with more conservative score
@@ -787,6 +788,7 @@ class Parameters {
     PARAMETER(PARAM_ALT_ALIGNMENT)
     PARAMETER(PARAM_GAP_OPEN)
     PARAMETER(PARAM_GAP_EXTEND)
+    PARAMETER(PARAM_PROTEOME_SIMILARITY)
 #ifdef GAP_POS_SCORING
     PARAMETER(PARAM_GAP_PSEUDOCOUNT)
 #endif
@@ -1078,6 +1080,7 @@ class Parameters {
     std::vector<MMseqsParameter*> threadsandcompression;
 
     std::vector<MMseqsParameter*> alignall;
+    std::vector<MMseqsParameter*> alignproteome;
     std::vector<MMseqsParameter*> align;
     std::vector<MMseqsParameter*> rescorediagonal;
     std::vector<MMseqsParameter*> alignbykmer;
@@ -1115,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/CMakeLists.txt

@@ -80,5 +80,6 @@ set(util_source_files
         util/proteinaln2nucl.cpp
         util/versionstring.cpp
         util/diskspaceavail.cpp
+        util/alignproteome.cpp
         PARENT_SCOPE
         )