PMSF might not work correctly in the presence of identical seqs and a potential solution to it

[evolbeginner]

Dear IQ-Tree:

I write to report a potential bug related to PMSF and an easy, likely working solution.

The data are available in the folders iqtree/ and iqtree_no_keepident/. Briefly, for my dataset where there are identical sequences, PMSF seems to works incorrectly if -keep-ident is NOT specified but by -keep-ident the problem can be solved.

Preparation

cd bs_inBV_bestfit_issue/split/0/iqtree/
then
cd ../iqtree_no_keepident/

Purpose

This directory compares several IQ-TREE runs on the same alignment (../combined.phy) and the same fixed tree (ref.tre), under different model / site-frequency settings.

The main point is to show that -keep-ident might affect PMSF-related runs.

---

Shared setup

All runs use:

• ../combined.phy
• ref.tre
• 30 protein sequences
• 100 sites
• fixed-tree evaluation (-te ref.tre), not tree search

So differences in BEST SCORE FOUND come from model/options, not from different searched topologies.

---

Run types

Prefix	Main command idea	Meaning
LG+G	-m LG+G	plain LG + Gamma
LG+C10+G	-m LG+C10+G	LG + C10 mixture + Gamma
LG+C10+G_ft	-m LG+C10+G -ft guide.treefile	infer PMSF/sitefreq from guide.treefile, then evaluate
LG+C10+G_fs	-m LG+C10+G -fs LG+C10+G_ft.sitefreq	reuse inferred sitefreq file
LG+G_fs	-m LG+G -fs LG+C10+G_ft.sitefreq	LG+G with the same imported sitefreq file

Notes:

• *_ft generates LG+C10+G_ft.sitefreq
• *_fs consumes that .sitefreq file
• guide.treefile is only used as input for -ft
• If i understand correctly, In theory, LG+C10+G_ft, LG+C10+G_fs, and LG+G_fs should result in the same lnL (this is what was observed when -keep-ident is specified; see below).

---

Results

Without -keep-ident

In iqtree_no_keepident/:

grep 'BEST SCORE' *log

gives:

LG+C10+G_fs_noki.log:BEST SCORE FOUND : -1634.495
LG+C10+G_ft_noki.log:BEST SCORE FOUND : -1634.495
LG+C10+G_noki.log:BEST SCORE FOUND : -1634.495
LG+G_fs_noki.log:BEST SCORE FOUND : -1666.411
LG+G_noki.log:BEST SCORE FOUND : -1666.411

This means:

• LG+C10+G, LG+C10+G_ft, and LG+C10+G_fs all collapse to the same lnL.
• • LG+G and LG+G_fs also collapse to the same lnL
• i.e. PMSF / imported site frequencies appear to have no effect without -keep-ident

---

With -keep-ident

In iqtree/:

grep 'BEST SCORE' *log

gives:

LG+C10+G_fs.log:BEST SCORE FOUND : -1475.677
LG+C10+G_ft.log:BEST SCORE FOUND : -1475.677
LG+C10+G.log:BEST SCORE FOUND : -1634.495
LG+G_fs.log:BEST SCORE FOUND : -1475.677
LG+G.log:BEST SCORE FOUND : -1666.412

This means:

• LG+C10+G_ft, LG+C10+G_fs, and LG+G_fs all give the improved PMSF-like lnL: -1475.677
• while the ordinary non-PMSF runs remain:
  • LG+C10+G: -1634.495
  • LG+G: -1666.412

So with -keep-ident, the PMSF-related runs behave as expected.

---

Interpretation

The issue appears to be tied to identical-sequence handling:

• without -keep-ident, the PMSF/sitefreq-based runs do not produce their distinct improved likelihoods
• with -keep-ident, they do

In other words, the problematic behavior seems to disappear once -keep-ident is enabled.

---

bs_inBV_bestfit_issue.tar.gz

cheers,
sishuo

[StefanFlaumberg]

Hi Sishuo @evolbeginner,

The issue really exists, I already spotted it before. In fact, I remember making a fix to it during the summer, but I haven't published the code.

In theory, LG+C10+G_ft, LG+C10+G_fs, and LG+G_fs should result in the same lnL

Yes, they should. In fact, LG+C10+G_fs here is a misleading construction, which should be banned and warned about, since this model never really uses the +C10 mixture at any step -- it is a site-specific model with given site frequencies, not a mixture one.

The technical cause of the issue of not using site-specific parameters in case the --keep-ident option is not enabled is that while removing identical sequences requires removing and reinitializing the Alignment object, the functions performing such operations do not copy some critical parameters of the original Alignment object, such as site-specific frequencies. This leaves the resultant Alignment without site-specific frequencies, which does not allow the ModelSet model to be initialized (the PMSF implementation) and leads instead to initializing a simple DNA or Protein model.
This particular issue stems from the problem in the Alignment::extractSubAlignment function. But there are quite a few other functions that copy alignments, all suffering the same problem of not copying the site-specific parameters.

I guess I will add the comprehensive fix either as part of my PMSFR model PR or as a separate PR.

Best,
Stefan

[bqminh]

Thanks for spotting this, Stefan. May I suggest that you make another pull request, separately from the PMSF feature. That PR would be much simpler, and easier for me to review it.

[StefanFlaumberg]

@bqminh,

In fact, the fixes that I made have a number of subtleties depending on the PMSFR implementation, such as new variable names and copying site-specific rates (instead of only the frequencies) to the reinitialized alignments.
But generally, I agree that a lighter PR would be a better solution here. And the issue itself is quite serious, since it appears that in a standard analysis users can end up running a simple model while thinking they're running a PMSF model.
Maybe, I could reapply the main changes to the current master branch now for a small PR and later fix the subtleties for the PMSFR PR. I'll see what I can do!

[evolbeginner]

thx much stefan! actually -m LG+Cxx+G -fs sitefreq_file is used in the official tutorial which, however, as you pointed out, doesn't make too much sense but potentially confuses users. this is quite a lot of work, and i believe you're very right that it could potentially affect many analyses. thx so much!

[evolbeginner]

I guess there might be other things that work correctly only if -keep-ident is specified. do you also know any, @StefanFlaumberg ? one might be the hessian matrix generation. anyway, -keep-ident seems the safest way.

[StefanFlaumberg]

@evolbeginner,

actually -m LG+Cxx+G -fs sitefreq_file is used in the official tutorial

It is, yet it shouldn't :-)
I remember that having read the tutorial for the first time, I was very confused about it myself. Like, "what info does the +Cxx specifier convey, given we've already provided a .sitefreq profile?". I believe this tutorial example should actually be fixed.

I guess there might be other things that work correctly only if -keep-ident is specified

I haven't checked this. Regarding the branch length Hessian matrix, I see that it's computed after removing identical sequences and before inserting them back, so using --keep-ident could make a difference, yet there is nothing a priori suspicious here.
But you can check it as with the PMSF case: if the results differ significantly, then there'll be a reason to look closer to it. Here in PMSF, the problem stems not from identical sequence management itself, but from shortcomings of alignment copying functions.

[StefanFlaumberg]

@bqminh,
I've looked through the bug and my previous fix attempts once again. The problem's turned out to be somewhat more complicated than I thought, but yes, it can be solved within a separate PR.
I hope to drop the bug-fixing PR during the coming week or so...