Unable to reproduce Precision-Recall plot Supp. Fig. 5

[rcedgar]

I have tried reproducing Supp. Fig. 5 using the scripts in this repository and also my own code but I get quite different results. For example, as shown in the figure above I find TM-align is better than DALI on Superfamily over the entire range, while your figure shows TM-align to be substantially worse (Fig. S5 on the left, my results on the right). My plot on the right was generated using your data and scripts as follows.

Hits downloaded from https://wwwuser.gwdguser.de/~compbiol/foldseek/scop.benchmark.result.tar.gz

sort -rgk3 ../alns/TMalign.txt > TMalign.sorted.txt
sort -rgk3 ../alns/dali.txt > dali.sorted.txt
bench.fdr.noselfhit.awk TMalign.sorted.txt scop_lookup.fix.tsv <(cat $TMalign.sorted.txt) > TM-align.rocx
bench.fdr.noselfhit.awk dali.sorted.txt scop_lookup.fix.tsv <(cat dali.sorted.txt) > dali.rocx

Plot column PREC_SFAM (y axis) vs. RECALL_SFAM (x axis).

As you can see, the plot for DALI looks right but TM-align looks very wrong.

Any help in resolving this discrepancy will be much appreciated.

Also, the calculation of precision and recall in bench.fdr.noselfhit.awk appears to be making corrections compared to the standard formulas, but I don't understand how it works. Can you clarify? In particular, what is the variable norm is doing? Thanks!

[rcedgar]

After looking more closely at the source code, as best I can tell the calculations of both precision and recall are not correct.

Line 5 in bench.fdr.noselfhit.awk outputs 6 column headings:

                    1           2           3            4             5             6
      print "PREC_FAM","PREC_SFAM","PREC_FOLD","RECALL_FAM","RECALL_SFAM","RECALL_FOLD";

Line 37 in bench.fdr.noselfhit.awk outputs 7 columns:

                                     1                     2                     3                 4                  5                  6       7
NR %1000 == 0{ print tp_fam/(tp_fam+fp), tp_sfam/(tp_sfam+fp), tp_fold/(tp_fold+fp), tp_fam / queries, tp_sfam / queries, tp_fold / queries, tp_fam; }

Matching column headings to expressions in line 37:

1. PREC_FAM     = tp_fam/(tp_fam+fp)
2. PREC_SFAM    = tp_sfam/(tp_sfam+fp)
3. PREC_FOLD    = tp_fold/(tp_fold+fp)

4. RECALL_FAM   = tp_fam / queries # queries = constant = 3566
5. RECALL_SFAM  = tp_sfam / queries
6. RECALL_FOLD  = tp_fold / queries

Definitions of precision and recall are (see e.g. https://scikit-learn.org/1.5/auto_examples/model_selection/plot_precision_recall.html):

# TP = number of true positive hits above the threshold
# FP = number of false positive hits above the threshold
# FN = number of false negatives at the threshold

Precision = TP/(TP + FP)
Recall = TP/(TP + FN)

The calculation of precision appears to be wrong for family and superfamily because fp is always calculated for fold (see line 38, it is increased if id2fold[$1] != id2fold[$2]).

The calculation of recall appears to be wrong because the divisor should be (TP+FN) but is constant (always 3556).

[martin-steinegger]

Thanks for flagging this. We will have a look at it. But since the postdoc left the lab it might need some time.
Regarding queries, we are computing a weighted ROC, meaning each query can contribute in total 1 if it recalls all possible TPs (TP=1/foldSize(norm)). Therefore, the number of queries is calculated as TP + FN. Precision is computed as expected: Precision = TP / (TP + FP) but we also apply weighting.

As for the TM-align results, we tried various methods to sort the hits. A reviewer recommended using the average of qTM and tTM scores, which indeed worked best. It’s possible that the uploaded file does not reflect this averaging.

[QY6616]

Do you emersion DALI and TM - align index, is to use the https://wwwuser.gwdguser.de/~compbiol/foldseek/scop.benchmark.result.tar.gz related script? Because I didn't find in the https://github.com/steineggerlab/foldseek-analysis DALI and TM - align related scripts @rcedgar

[rcedgar]

@QY6616 sorry, I don't understand the question

[QY6616]

@QY6616抱歉，我不明白这个问题

@rcedgar What I mean is: How did you reproduce the PR curve of DALI and TM-align? I also want to reproduce the PR curve of DALI and TM-align on the SCOPE dataset, but I haven't found the relevant scripts.

[rcedgar]

foldseek-analysis/scopbenchmark/scripts/bench.fdr.noselfhit.awk
foldseek-analysis/scopbenchmark/scripts/bench.noselfhit.awk
foldseek-analysis/training/roc1.awk

[rcedgar]

fyi the corresponding scripts & figures for the reseek paper are here https://github.com/rcedgar/reseek_bench

[QY6616]

fyi the corresponding scripts & figures for the reseek paper are here https://github.com/rcedgar/reseek_bench

@rcedgar Thank you very much for your help, but the link you provided cannot be opened. Could you please check it?

[rcedgar]

oops! that repo was still private, it should have been made public when the paper was published -- I fixed that, you should be able to reach it now.

[QY6616]

@rcedgar It is accessible now! Actually, what I want to ask is: as shown in the figure, where can I download the ../data/scop.file.full.name.list and ./TM-align/TMalign files in foldseek-analysis? I see that you provided a download link: https://wwwuser.gwdguser.de/~compbiol/foldseek/scop.benchmark.result.tar.gz. Does this file contain the relevant files for TMalign?

[rcedgar]

Script foldseek_rocx.bash in rcedgar/reseek_bench runs the awk script from steineggerlab/foldsee-analysis , also reseek_bench/awk has a port of their awk code to python to make sure I understood the logic by ensuring that it generates the same results. I don't remember all the details, IIRC I had to make some changes and reconstruct some helper files, I don't think the code in steineggerlab/foldsee-analysis runs as-is. The https://wwwuser.gwdguser.de/~compbiol/foldseek/scop.benchmark.result.tar.gz download has alignments for DALI, TM-align and the other methods, I don't remember what else it has offhand.

[QY6616]

Thank you for your help, I will give it a try.

[QY6616]

@rcedgar Hello, I obtained the same PR curves for TMalign and DALI using the same method as you, but when I tried to reproduce the MMseq2 results, I encountered the data shown in the image. Did I overlook any details, or is there an issue with the data itself?

[rcedgar]

The image seems to show a table of numbers, I don't really understand what you are asking. If you can't reproduce a figure from the Foldseek paper, then this seems to confirm the issue here. If so, I suggest it would be useful to show the MMseqs curve you obtained side-by-side to clarify the discrepancy, similar to the image I showed at the start of this issue. FWIW, I believe that some of the curves in the Foldseek paper are wrong because some curves are reproduced by their code but others are not; this could possibly be explained by something like a copy/paste error. The TM-align example is striking because the published figure implies that TM-align is much worse than DALI, while I find that TM-align is better than DALI using both the code here and also my own code.