change gap penalties, add citation

gamcil · gamcil · commit 36eb84ca4db8 · 2025-04-03T11:36:01.000+09:00
diff --git a/src/commons/FoldmasonParameters.cpp b/src/commons/FoldmasonParameters.cpp
@@ -13,7 +13,8 @@ FoldmasonParameters::FoldmasonParameters() :
         PARAM_PAIR_THRESHOLD(PARAM_PAIR_THRESHOLD_ID, "--pair-threshold", "LDDT pair threshold", "% of pair subalignments with LDDT information [0.0,1.0]",typeid(float), (void *) &pairThreshold, "^0(\\.[0-9]+)?|1(\\.0+)?$"),
         PARAM_REPORT_COMMAND(PARAM_REPORT_COMMAND_ID, "--report-command", "", "", typeid(std::string), (void *) &reportCommand, ""),
         PARAM_REPORT_PATHS(PARAM_REPORT_PATHS_ID, "--report-paths", "", "", typeid(bool), (void *) &reportPaths, ""),
-        PARAM_REFINE_SEED(PARAM_REFINE_SEED_ID, "--refine-seed", "Random number generator seed", "Random number generator seed", typeid(int), (void *) &refinementSeed, "^([-]?[0-9]*)$")
+        PARAM_REFINE_SEED(PARAM_REFINE_SEED_ID, "--refine-seed", "Random number generator seed", "Random number generator seed", typeid(int), (void *) &refinementSeed, "^([-]?[0-9]*)$"),
+        PARAM_ONLY_SCORING_COLS(PARAM_ONLY_SCORING_COLS_ID, "--only-scoring-cols", "Normalise LDDT by no. scoring columns", "Normalise LDDT by no. scoring columns", typeid(bool), (void *) &onlyScoringCols, "")
 {
     // structuremsa
     structuremsa.push_back(&PARAM_WG);
@@ -40,6 +41,7 @@ FoldmasonParameters::FoldmasonParameters() :
     structuremsa.push_back(&PARAM_NO_COMP_BIAS_CORR);
     structuremsa.push_back(&PARAM_V);
     structuremsa.push_back(&PARAM_REFINE_SEED);
+    structuremsa.push_back(&PARAM_ONLY_SCORING_COLS);
 
     structuremsacluster = combineList(structuremsacluster, structuremsa);
 
@@ -50,6 +52,7 @@ FoldmasonParameters::FoldmasonParameters() :
     msa2lddt.push_back(&PARAM_V);
     msa2lddt.push_back(&PARAM_REPORT_COMMAND);
     msa2lddt.push_back(&PARAM_REPORT_PATHS);
+    msa2lddt.push_back(&PARAM_ONLY_SCORING_COLS);
 
     // refinemsa
     refinemsa = combineList(refinemsa, structuremsa);
@@ -78,8 +81,9 @@ FoldmasonParameters::FoldmasonParameters() :
     filterMsa = 0;
     compBiasCorrection = 0;
     refinementSeed = -1;
-    gapOpen = 25;
-    gapExtend = 3;
+    gapOpen = 20;
+    gapExtend = 2;
+    onlyScoringCols = false;
 
-    citations.emplace(CITATION_FOLDMASON, " << TODO >> ");
+    citations.emplace(CITATION_FOLDMASON, "Gilchrist, C. L. M., Mirdita, M. & Steinegger, M. Multiple Protein Structure Alignment at Scale with FoldMason. bioRxiv, doi.org/10.1101/2024.08.01.606130 (2024)");
 }
diff --git a/src/commons/FoldmasonParameters.h b/src/commons/FoldmasonParameters.h
@@ -41,6 +41,7 @@ class FoldmasonParameters : public LocalParameters {
     PARAMETER(PARAM_REPORT_COMMAND)
     PARAMETER(PARAM_REPORT_PATHS)
     PARAMETER(PARAM_REFINE_SEED)
+    PARAMETER(PARAM_ONLY_SCORING_COLS)
     
     std::string guideTree;
     bool recomputeScores;
@@ -52,5 +53,6 @@ class FoldmasonParameters : public LocalParameters {
     bool reportPaths;
     float pairThreshold;
     int refinementSeed;
+    bool onlyScoringCols;
 };
 #endif
diff --git a/src/strucclustutils/msa2lddt.cpp b/src/strucclustutils/msa2lddt.cpp
@@ -245,7 +245,8 @@ std::tuple<std::vector<float>, std::vector<int>, float, int> calculate_lddt(
     std::vector<size_t> &subset,
     std::vector<size_t> &keys,
     DBReader<unsigned int> * seqDbrCA,
-    float pairThreshold
+    float pairThreshold,
+    bool onlyScoringCols
 ) {
     // mapping:
     // cigar vectors can be any size
@@ -375,6 +376,7 @@ std::tuple<std::vector<float>, std::vector<int>, float, int> calculate_lddt(
     std::vector<int> colCounts = countColumns(cigars, subset, alnLength);
     float scaledSum = 0.0;
     int numCols = 0;
+    int numColsScoring = 0;
     for (size_t i = 0; i < perColumnCount.size(); i++) {
         // float pairSupport = perColumnCount[i] / static_cast<float>(numPairs);
         float residuesInColumn = colCounts[i];
@@ -392,13 +394,16 @@ std::tuple<std::vector<float>, std::vector<int>, float, int> calculate_lddt(
             //perColumnScore[i] /= perColumnCount[i];  // get mean LDDT for this column
             perColumnScore[i] = (2 * score * occupancy) / (score + occupancy);
             scaledSum += perColumnScore[i];
+            numColsScoring++;
         } else {
             perColumnScore[i] = 0.0;
         }
         // scaledSum += colscore;
         numCols++;
     }
-    float lddtScore = (numCols > 0) ? scaledSum / static_cast<float>(numCols) : 0.0;
+    float lddtScore = onlyScoringCols
+        ? ((numColsScoring > 0) ? scaledSum / static_cast<float>(numColsScoring) : 0.0f)
+        : ((numCols > 0) ? scaledSum / static_cast<float>(numCols) : 0.0f);
     return std::make_tuple(perColumnScore, perColumnCount, lddtScore, numCols);
 }
 
@@ -464,7 +469,8 @@ float getLDDTScore(
     DBReader<unsigned int> &seqDbr3Di,
     DBReader<unsigned int> &seqDbrCA,
     std::string msa,
-    float pairThreshold
+    float pairThreshold,
+    bool onlyScoringCols
 ) {
     KSeqWrapper* kseq = new KSeqBuffer(msa.c_str(), msa.length());
     int alnLength = 0;
@@ -479,7 +485,7 @@ float getLDDTScore(
     std::vector<int>   perColumnCount;
     float lddtScore;
     int numCols;
-    std::tie(perColumnScore, perColumnCount, lddtScore, numCols) = calculate_lddt(cigars_aa, inds, inds, &seqDbrCA, pairThreshold);
+    std::tie(perColumnScore, perColumnCount, lddtScore, numCols) = calculate_lddt(cigars_aa, inds, inds, &seqDbrCA, pairThreshold, onlyScoringCols);
 
     return lddtScore;
 }
@@ -531,7 +537,7 @@ int msa2lddt(int argc, const char **argv, const Command& command, int makeReport
     std::iota(subset.begin(), subset.end(), 0);
     
     if (caExist) {
-        std::tie(perColumnScore, perColumnCount, lddtScore, numCols) = calculate_lddt(cigars_aa, subset, indices, seqDbrCA, par.pairThreshold);
+        std::tie(perColumnScore, perColumnCount, lddtScore, numCols) = calculate_lddt(cigars_aa, subset, indices, seqDbrCA, par.pairThreshold, par.onlyScoringCols);
         std::string scores;
         for (float score : perColumnScore) {
             if (scores.length() > 0) scores += ",";
diff --git a/src/strucclustutils/msa2lddt.h b/src/strucclustutils/msa2lddt.h
@@ -24,7 +24,8 @@ std::tuple<std::vector<float>, std::vector<int>, float, int> calculate_lddt(
     std::vector<size_t> &subset,
     std::vector<size_t> &keys,
     DBReader<unsigned int> * seqDbrCA,
-    float pairThreshold
+    float pairThreshold,
+    bool scoringColsOnly
 );
 
 double calculate_lddt_pair(
@@ -49,7 +50,8 @@ float getLDDTScore(
     DBReader<unsigned int> &seqDbr3Di,
     DBReader<unsigned int> &seqDbrCA,
     std::string msa,
-    float pairThreshold
+    float pairThreshold,
+    bool scoringColsOnly
 );
 
 #endif
diff --git a/src/strucclustutils/refinemsa.cpp b/src/strucclustutils/refinemsa.cpp
@@ -259,7 +259,8 @@ void refineMany(
     std::string qid,
     float pairThreshold,
     std::vector<size_t> indices,
-    int seed
+    int seed,
+    bool onlyScoringCols
 ) {
     std::cout << "Running " << iterations << " refinement iterations\n";
 
@@ -268,7 +269,7 @@ void refineMany(
         subset[i] = i;
     }
 
-    float prevLDDT = std::get<2>(calculate_lddt(cigars_aa, subset, indices, seqDbrCA, pairThreshold));
+    float prevLDDT = std::get<2>(calculate_lddt(cigars_aa, subset, indices, seqDbrCA, pairThreshold, onlyScoringCols));
     float initLDDT = prevLDDT;
     std::cout << "Initial LDDT: " << prevLDDT << '\n';
 
@@ -304,7 +305,7 @@ void refineMany(
             sequences_aa, sequences_ss,
             rng
         );
-        float lddtScore = std::get<2>(calculate_lddt(cigars_new_aa, subset, indices, seqDbrCA, pairThreshold));
+        float lddtScore = std::get<2>(calculate_lddt(cigars_new_aa, subset, indices, seqDbrCA, pairThreshold, onlyScoringCols));
         // std::cout << std::fixed << std::setprecision(4) << "New LDDT: " << lddtScore << '\t' << "(" << i + 1 << ")\n";
         // for (std::vector<Instruction> &ins : cigars_new_aa) {
         //     std::cout << expand(ins) << '\n';
@@ -407,7 +408,7 @@ int refinemsa(int argc, const char **argv, const Command& command) {
         structureSmithWaterman, par.refineIters, par.compBiasCorrection, par.wg, par.filterMaxSeqId,
         par.qsc, par.Ndiff, par.covMSAThr,
         par.filterMinEnable, par.filterMsa, par.gapExtend.values.aminoacid(), par.gapOpen.values.aminoacid(),
-        par.maxSeqLen, par.qid, par.pairThreshold, indices, par.refinementSeed
+        par.maxSeqLen, par.qid, par.pairThreshold, indices, par.refinementSeed, par.onlyScoringCols
     );
     
     // Write final MSA to file
diff --git a/src/strucclustutils/refinemsa.h b/src/strucclustutils/refinemsa.h
@@ -30,7 +30,8 @@ void refineMany(
     std::string qid,
     float pairThreshold,
     std::vector<size_t> indices,
-    int seed
+    int seed,
+    bool onlyScoringCols
 );
 void deleteGapCols(std::vector<std::string> &sequences);
 void buildSubMSA(std::vector<std::string> &headers, std::vector<std::string> &sequences, std::string &subMSA);
diff --git a/src/strucclustutils/structuremsa.cpp b/src/strucclustutils/structuremsa.cpp
@@ -1642,7 +1642,7 @@ int structuremsa(int argc, const char **argv, const Command& command, bool preCl
             par.refineIters, par.compBiasCorrection, par.wg, par.filterMaxSeqId, par.qsc,
             par.Ndiff, par.covMSAThr, par.filterMinEnable, par.filterMsa, par.gapExtend.values.aminoacid(),
             par.gapOpen.values.aminoacid(), par.maxSeqLen, par.qid, par.pairThreshold, msa.dbKeys,
-            par.refinementSeed
+            par.refinementSeed, par.onlyScoringCols
         );
     }
 }

Original file line number	Diff line number	Diff line change
`@@ -1642,7 +1642,7 @@ int structuremsa(int argc, const char **argv, const Command& command, bool preCl`
`1642`	`1642`	`par.refineIters, par.compBiasCorrection, par.wg, par.filterMaxSeqId, par.qsc,`
`1643`	`1643`	`par.Ndiff, par.covMSAThr, par.filterMinEnable, par.filterMsa, par.gapExtend.values.aminoacid(),`
`1644`	`1644`	`par.gapOpen.values.aminoacid(), par.maxSeqLen, par.qid, par.pairThreshold, msa.dbKeys,`
`1645`		`- par.refinementSeed`
	`1645`	`+ par.refinementSeed, par.onlyScoringCols`
`1646`	`1646`	`);`
`1647`	`1647`	`}`
`1648`	`1648`	`}`