alignproteome added by martin-steinegger · Pull Request #875 · soedinglab/MMseqs2

martin-steinegger · 2024-08-18T08:51:22Z

No description provided.

martin-steinegger · 2024-08-18T08:52:26Z

src/MMseqsBase.cpp

                                                           {"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",


I would improve this text.

martin-steinegger · 2024-08-18T08:52:45Z

src/MMseqsBase.cpp

+        {"alignproteome",             alignproteome,             &par.alignproteome,             COMMAND_ALIGNMENT,
+                "Within-result all-vs-all gapped local alignment",
+                NULL,
+                "Martin Steinegger <[email protected]>",


Please add you @Gyuuul2

martin-steinegger · 2024-08-18T08:53:35Z

src/commons/ClusterSpecies.cpp

This file looks incomplete

martin-steinegger · 2024-08-18T08:53:55Z

src/commons/ClusterSpecies.h

martin-steinegger · 2024-08-18T08:54:28Z

src/commons/DBReader.cpp

        }
    }
    if (dataMode & USE_LOOKUP || dataMode & USE_LOOKUP_REV) {
+        Debug(Debug::INFO) << "ReadLookup file: " << dataFileName << "\n"; //gyuri


I would remove this print out.

martin-steinegger · 2024-08-18T08:55:03Z

src/commons/DBReader.cpp

        }
        char* lookupDataChar = (char *) lookupData.getData();
        size_t lookupDataSize = lookupData.size();
+        Debug(Debug::INFO) << "Lookup Data size is " << lookupDataSize << "\n"; //gyuri


I would remove this print out.

martin-steinegger · 2024-08-18T08:55:08Z

src/commons/DBReader.cpp

        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


I would remove this print out.

martin-steinegger · 2024-08-18T09:02:00Z

src/commons/CMakeLists.txt

        commons/BaseMatrix.cpp
        commons/Command.cpp
        commons/CommandCaller.cpp
+        # commons/ClusterSpecies.cpp


Please remove it.

martin-steinegger · 2024-08-18T09:02:59Z

src/alignment/Matcher.h


    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);


Why is this needed?

martin-steinegger · 2024-08-18T09:03:40Z

src/util/alignall.cpp


-    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> 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);


I assume this is not needed anymore

martin-steinegger · 2024-08-18T09:03:46Z

src/util/alignall.cpp

    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';


martin-steinegger · 2024-08-18T09:12:43Z

src/util/alignproteome.cpp

+    }
+};
+
+struct __attribute__((__packed__)) memberProteinEntry{


MemberProteinEntry

milot-mirdita · 2024-08-19T08:54:09Z

src/commons/DBReader.h

    size_t getSize() const;

-    unsigned int getProteomeTotalLen(size_t id); //gyuri
+    unsigned int getProteomeTotalLen(size_t id); 


I would name this getSourceTotalLen or getSetTotalLen to be consistent with the .source naming

milot-mirdita · 2024-08-19T08:55:37Z

src/workflow/EasyAlignproteome.cpp

+}
+
+// void setEasyAlignproteomeMustPassAlong(Parameters *p){
+


You need to implement this whenever you use the createParameterString(..., true) parameter. Or it won't get passed along.

milot-mirdita · 2024-08-19T08:56:16Z

src/workflow/EasyAlignproteome.cpp

+    par.filenames.pop_back();
+
+    CommandCaller cmd;
+    cmd.addVariable("ALIGNPROTEOME_PAR", par.createParameterString(par.alignproteome,true).c_str()); // what?


CLUSTER_PAR, THREADS_PAR etc are missing

milot-mirdita · 2024-08-19T08:57:35Z

src/workflow/EasyAlignproteome.cpp

+        hash = FileUtil::getHashFromSymLink(tmpDir + "/latest");
+    }
+    tmpDir = FileUtil::createTemporaryDirectory(tmpDir, hash);
+    par.filenames.pop_back();


You need to pop_back() the output path too and pass it as an environment variable. Or else it will also be passed to createdb.

…of different fasta files

…usage to createdb and easy-search in MMseqsBase.cpp.

milot-mirdita · 2024-08-27T15:57:40Z

src/commons/IndexReader.h

    ) : sequenceReader(NULL), index(NULL) {
        int targetDbtype = FileUtil::parseDbType(dataName.c_str());
        if (Parameters::isEqualDbtype(targetDbtype, Parameters::DBTYPE_INDEX_DB)) {
-            index = new DBReader<unsigned int>(dataName.c_str(), (dataName + ".index").c_str(), 1, DBReader<unsigned int>::USE_DATA|DBReader<unsigned int>::USE_INDEX);


what happened here?

milot-mirdita · 2024-08-27T16:01:48Z

src/util/proteomecluster.cpp

+        const unsigned int ProteomeId = lookup[i].fileNumber;
+        const unsigned int ProteinId = lookup[i].id;
+        ProteomeList[ProteomeId].addSeqLen(tProteinDB.getSeqLen(ProteinId));
+        if (ProteomeList[ProteomeId].proteomeKey == -1){


style: please add spaces around if:
if[space](...)[space]{

milot-mirdita · 2024-08-27T16:03:23Z

src/util/proteomecluster.cpp

+    std::vector<unsigned int> memberKeys;
+    memberKeys.reserve(50); // store key for every protein in a cluster
+
+    std::vector<bool> isProteomeInCluster(localMemCount.size(), false); ;


double ;

also careful about bool vectors. They are not multi-threading friendly. Here its fine, as you don't do any multi-threading

milot-mirdita · 2024-08-27T16:03:36Z

src/util/proteomecluster.cpp

+    if (memberKeys.size() > 1) { //If not a singleton cluster
+        ClusterEntry eachClusterRep(memberKeys.size());
+        //init MemberProteinEntry and add it to memberProteins vector
+        for (auto& eachMemberKey : memberKeys){


style: missing space after the for

milot-mirdita · 2024-08-27T16:03:52Z

src/util/proteomecluster.cpp

+    for (size_t idx = 0; idx < ProteomeList.size(); idx++){
+        if (ProteomeList[idx].isCovered()) {
+            continue;
+        }else{


style: missing spaces around else

milot-mirdita · 2024-08-27T16:05:58Z

src/util/proteomecluster.cpp

+    DBReader<unsigned int>::LookupEntry* tLookup = tProteinDB.getLookup();
+    const size_t tLookupSize = tProteinDB.getLookupSize();
+    unsigned int totalProteomeNumber = tLookup[tLookupSize - 1].fileNumber;
+    std::vector<ProteomeEntry> ProteomeList(totalProteomeNumber + 1);


style: should be lowercase: proteomeList

milot-mirdita · 2024-08-27T16:07:02Z

src/util/proteomecluster.cpp

+
+    int gapOpen, gapExtend;
+    BaseMatrix *subMat;
+    subMat = new SubstitutionMatrix(par.scoringMatrixFile.values.aminoacid().c_str(), 2.0, par.scoreBias);


You can avoid the explicit allocation here:

SubstitutionMatrix subMat(par.scoringMatrixFile.values.aminoacid().c_str(), 2.0, par.scoreBias); ...

milot-mirdita · 2024-08-27T16:08:16Z

src/util/proteomecluster.cpp

+        #pragma omp critical
+        {
+            for (size_t idx=0; idx < localMemCount.size(); idx++){
+                ProteomeList[idx].incrementMemCount(localMemCount[idx]);


if you use __sync_fetch_and_add inside incrementMemCount you don't need the critical

milot-mirdita · 2024-08-27T16:09:38Z

src/util/proteomecluster.cpp

+            #pragma omp critical
+            {
+                for (size_t idx = 0; idx < ProteomeList.size(); idx++) {
+                    ProteomeList[idx].addSeqId(localSeqIds[idx]);


here you might need the critical as atomic floats can be tricky to do correctly

milot-mirdita · 2024-08-27T16:10:08Z

src/util/proteomecluster.cpp

+    proteomeClustWriter.close();
+    proteinClustWriter.close();
+    tProteinDB.close();
+    delete subMat;


if you rewrite the subMat init above, you can remove this

…mber proteome is detected)

…teome handling update

alignproteome added

1cb2074

martin-steinegger commented Aug 18, 2024

View reviewed changes

martin-steinegger commented Aug 19, 2024

View reviewed changes

add workflow, revise based on feedback

12c0c81

milot-mirdita reviewed Aug 19, 2024

View reviewed changes

elpis51613 and others added 8 commits August 25, 2024 18:49

Change createdb.cpp so that it takes in ".txt" file containing paths …

df80f07

…of different fasta files

Merge remote-tracking branch 'refs/remotes/origin/master'

1073ec9

Change filepath format from txt to tsv. Also add explanations of its …

3617ea7

…usage to createdb and easy-search in MMseqsBase.cpp.

createtsv - handle source identifier

3fc7fcd

proteomcluster - change writer / way to select repProteome

0487a1f

easy proteomecluster workflow

3f4fc91

indexreader - handle source identifier

eeabad3

update mmseqsbase and change module name into proteomecluster

9996821

milot-mirdita reviewed Aug 27, 2024

View reviewed changes

Gyuuul2 added 12 commits August 28, 2024 14:16

set alignment mode default

27d9f68

latest

e2ef229

original index

fd1f5d6

use sync_fetch_and_add when increment memcount

0d7db0f

change submat

1fe489e

add timer

86f6b15

Merge branch 'pr-879'

8344e5b

Remove __packed__ to resolve sync_fetch_and_add error

be5a53d

change chmod in createdb

8d3a38f

update proteome-similarity threshold parameter

bfa1911

Sort redundant proteomes by similarity score in proteome_cluster.tsv

7e912d2

delete redundant code - checking proteome is singleton(no redudant me…

618e2bf

…mber proteome is detected)

Gyuuul2 added 11 commits September 9, 2024 11:12

stop

e990773

jaccardindexTrial

e9ffe22

add normalized scoring

d14217f

add proteome relative(normalized) similarity threshold parameter

5e6ee6d

relative similarity threshold

743370d

latest version-bugfix,relativeProteomeSimilarity added, singleton pro…

44a6fca

…teome handling update

Change Relative Score formula

910bb7f

Remove Debug printout

e24b9e3

style change

ad1acdf

Latest Version

e026260

Add cluster count report

d9213ef


		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);

Conversation

martin-steinegger commented Aug 18, 2024

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants