Unable to reproduce the same FunDi results with the same command and the same seed

[novigit]

IQTREE version:

IQ-TREE multicore version 2.2.5 COVID-edition for Linux 64-bit built Sep 15 2023

Command line:

TAXA="xd-EG_Diplonema_sp,xd-EG_Bodo_saltans,am-AM_Dictyostelium_discoideum,di-CH_Micromonas_commoda,am-AM_Echinamoeba_exudans,di-CH_Dunaliella_salina,di-AL_Neospora_caninum,xd-HL_Acrasis_kona,am-FU_Catenaria_anguillulae,xd-EG_Rhabdomonas_costata,di-CH_Chlamydomonas_reinhardtii,di-PL_Selaginella_moellendorffii,xd-HL_Naegleria_lovaniensis,di-CH_Volvox_carteri,di-PL_Phalaenopsis_equestris,di-PL_Chara_braunii,am-FU_Ustilago_maydis,xm-PX_Paratrimastix_eleionoma,di-RP_Rhodosorus_marinus,xd-EG_Trypanosoma_sp,di-ST_Aplanochytrium_sp,xm-PX_Streblomastix_strix,di-PL_Oryza_brachyantha,xd-EG_Leptomonas_pyrrhocoris,di-AL_Amphidinium_massartii,di-ST_Ochromonas_sp,am-AM_Polysphondylium_pallidum,am-FU_Mortierella_elongata,am-AM_Cavostelium_apophysatum,am-FU_Aspergillus_aculeatus,xd-EG_Lacrimia_lanifica,xm-FR_Giardia_muris,di-RZ_Lotharella_sp,xm-PX_Paratrimastix_pyriformis,di-RP_Rhodella_sp,xd-JA_Andalucia_godoyi,di-CH_Monoraphidium_neglectum,xm-PB_Trichomonas_tenax,am-AM_Vermamoeba_vermiformis,am-FU_Allomyces_macrogynus,xd-EG_Hemistasia_phaeocysticola,xd-HL_Pharyngomonas_kirbyi,xd-EG_Leishmania_braziliensis,am-HO_Priapulus_caudatus,di-AL_Tetrahymena_thermophila,xm-PX_Trimastix_marina,xm-PB_Trichomitus_batrachorum,di-PL_Klebsormidium_nitens,di-RP_Gracilariopsis_chorda,am-AM_Dictyostelium_purpureum,di-RP_Porphyra_umbilicalis,am-FU_Aspergillus_nidulans,xd-HL_Percolomonas_cosmoAE1,xd-EG_Rhynchopus_humris,di-AL_Perkinsus_marinus,am-HO_Trichoplax_adhaerens,di-PL_Amborella_trichopoda,xd-HL_Naegleria_gruberi,di-ST_Phytophthora_sp,am-HO_Exaiptasia_pallida,xd-EG_Distigma_spP2,di-RZ_Partenskyella_glossopodia,am-FU_Rhizoclosmatium_globosum,xm-PB_Tritrichomonas_foetus,di-RZ_Bigelowiella_natans,xm-PB_Dientamoeba_fragilis,xm-FR_Carpediemonas_membranifera,di-CH_Bathycoccus_prasinos,di-ST_Ectocarpus_siliculosus,di-RP_Chondrus_crispus,xd-HL_Willaertia_magna,di-PL_Setaria_italica,xm-FR_Chilomastix_sp,xd-EG_Neobodo_designis,xd-EG_Ploeotia_vitrea,am-FU_Batrachochytrium_dendrobatidis,am-AM_Cryptodifflugia_operculata,di-CH_Coccomyxa_subellipsoidea,am-HO_Nematostella_vectensis,am-FU_Moniliophthora_roreri,xd-EG_Artemidia_motanka,xd-JA_Stygiella_incarcerata,am-FU_Gonapodya_prolifera,di-CH_Micractinium_conductrix,di-AL_Scrippsiella_Hangoei,xm-PB_Histomonas_meleagridis,am-HO_Monosiga_brevicollis,xm-PB_Trichomonas_gallinae,di-RZ_Plasmodiophora_brassicae,di-PL_Dioscorea_cayenensis,am-HO_Acanthoecalike_sp,am-HO_Danio_rerio,am-HO_Branchiostoma_floridae,di-CH_Auxenochlorella_protothecoides,xm-FR_Retortamonas_caviae,di-AL_Paramecium_tetraurelia,am-AM_Clastostelium_recurvatum,am-HO_Apis_mellifera,di-AL_Vitrella_brassicaformis,di-AL_Stentor_coeruleus,am-HO_Capsaspora_owczarzaki,am-HO_Hydra_vulgaris,di-AL_Polarella_glacialis,di-AL_Besnoitia_besnoiti,xd-HL_Percolomonas_cosmoWS,di-ST_Fistulifera_solaris,di-AL_Hammondia_hammondi,di-ST_Nannochloropsis_gaditana,di-AL_Toxoplasma_gondii,di-CH_Chlorella_variabilis,xd-EG_Eutreptiella_gymnastica,xd-EG_Namystynia_karyoxenos,di-AL_Stylonychia_lemnae,xm-FR_Dysnectes_brevis,am-FU_Spizellomyces_punctatus,am-HO_Stegodyphus_mimosarum,xd-JA_Seculamonas_ecuadoriensis,am-HO_Salpingoeca_rosetta,xm-FR_Kipferlia_bialata,di-AL_Togula_jolla,am-AM_Filamoeba_nolandi,di-ST_Aphanomyces_astaci,di-CH_Edaphochlamys_debaryana,am-AM_Copromyxa_protea,xd-EG_Sulcionema_specki,di-RZ_Chlorarachnion_reptans,xm-PB_Pentatrichomonas_hominis,di-PL_Ananas_comosus,am-FU_Laccaria_bicolor,am-FU_Rhizophagus_irregularis,di-PL_Triticum_urartu,xd-EG_Rhynchopus_euleeides,di-ST_Albugo_candida,di-RP_Porphyridium_purpureum,di-ST_Schizochytrium_aggregatum,di-ST_Thalassiosira_pseudonana,di-AL_Chromera_velia,am-AM_Acytostelium_subglobosum,xm-PX_Monocercomonoides_sp,xd-HL_Naegleria_fowleri,xd-EG_Strigomonas_culicis,xm-FR_Carpediemonas_likeNY0171,di-PL_Physcomitrella_patens,di-PL_Marchantia_polymorpha,am-AM_Protostelium_nocturnum"

iqtree2 -s supermatrix_loci_183.fasta -te c60.treefile -seed 137888 -m LG+fundi_C10+G4 -mdef re-optimized-model.nexus -a 0.6250003068091398 --fundi $TAXA,estimate --prefix supermatrix_loci_183.fullTree.C10.out/fundi_full_tree -nt 20 -prec 10 --quiet

The model fundi_C10 is specified in the re-optimized-model.nexus file, which states the per-class frequencies and weights.

I ran this exact run twice, and it gave me two markedly different results with regards to FunDi results:

Run1:

NOTE: 51844 MB RAM (50 GB) is required!
Estimate model parameters (epsilon = 0.010)
1. Initial log-likelihood: -6590864.099
2. Current log-likelihood: -6589798.467
3. Current log-likelihood: -6589429.568
4. Current log-likelihood: -6589373.746
5. Current log-likelihood: -6589364.260
6. Current log-likelihood: -6589362.464
7. Current log-likelihood: -6589362.104
8. Current log-likelihood: -6589362.030
Optimal log-likelihood: -6589362.018
Mixture weights: 0.099 0.308 0.067 0.040 0.097 0.089 0.096 0.029 0.097 0.078
Gamma shape alpha: 0.625
Parameters optimization took 8 rounds (503.782 sec)
Wrote distance file to...
BEST SCORE FOUND : -6589362.018
Optimizing FunDi model...
rho = 0.000
Central branch length: 1.593
Optimizing FunDi model parameters...
WARNING: Numerical underflow for non-rev lh-branch Noname
Best FunDi parameter rho: 0.494
Best FunDi central branch length: 0.450
FunDi log-likelihood: -6492752.884
Time taken to optimize FunDi model: 4.158 sec
Total tree length: 90.597

Run2:

NOTE: 51844 MB RAM (50 GB) is required!
Estimate model parameters (epsilon = 0.010)
1. Initial log-likelihood: -6590864.099
2. Current log-likelihood: -6589798.467
3. Current log-likelihood: -6589429.568
4. Current log-likelihood: -6589373.746
5. Current log-likelihood: -6589364.260
6. Current log-likelihood: -6589362.464
7. Current log-likelihood: -6589362.104
8. Current log-likelihood: -6589362.030
Optimal log-likelihood: -6589362.018
Mixture weights: 0.099 0.308 0.067 0.040 0.097 0.089 0.096 0.029 0.097 0.078
Gamma shape alpha: 0.625
Parameters optimization took 8 rounds (358.080 sec)
Wrote distance file to...
BEST SCORE FOUND : -6589362.018
Optimizing FunDi model...
rho = 0.000
Central branch length: 1.593
Optimizing FunDi model parameters...
WARNING: Numerical underflow for non-rev lh-branch Noname
Best FunDi parameter rho: 0.491
Best FunDi central branch length: 2.627
FunDi log-likelihood: -6495746.928
Time taken to optimize FunDi model: 2.660 sec
Total tree length: 92.774

Note that these FunDi results are different, particularly the inferred central branch length and log likelihood:

Run1:

Best FunDi parameter rho: 0.494
Best FunDi central branch length: 0.450
FunDi log-likelihood: -6492752.884

Run2:

Best FunDi parameter rho: 0.491
Best FunDi central branch length: 2.627
FunDi log-likelihood: -6495746.928

The only difference between these runs is that they were run on different nodes of the same cluster.

Run1:

Host:    perun17 (AVX2, FMA3, 755 GB RAM)
Command: iqtree2 -s supermatrix_loci_183.fasta -te c60.treefile -m LG+fundi_C10+G4 ....
Seed:    137888 (Using SPRNG - Scalable Parallel Random Number Generator)
Time:    Wed Dec  4 15:22:23 2024
Kernel:  AVX+FMA - 40 threads (64 CPU cores detected)

Run2:

Host:    perun23 (AVX2, FMA3, 2003 GB RAM)
Command: iqtree2 -s supermatrix_loci_183.fasta -te c60.treefile -m LG+fundi_C10+G4 ....
Seed:    137888 (Using SPRNG - Scalable Parallel Random Number Generator)
Time:    Tue Nov 26 11:12:53 2024
Kernel:  AVX+FMA - 20 threads (256 CPU cores detected)

I suspect that perhaps the starting seed ensures reproducibility up to optimization of the "normal" tree inference, but then when FunDi optimization starts, some other random element causes two different runs to diverge?

Also, there is the numerical underflow warning, could that be related somehow?

[bqminh]

We're currently still testing the FunDi model implementation with Andrew Roger's group. I can't really guarantee it works in all cases at the moment... I will update here once we make some progress.

[bqminh]

Thanks for the productive meeting!
Apart from --fundi-epsilon option, I have now introduced two new options:

• --fundi-init-rho to set the initial value for rho before the optimisation. The default is 0.5 instead of 0 like previous version.
• --fundi-init-branch to set the initial branch length before optimisation. The default is the branch length from the tree file.

This is now implemented in iqtree3 repo (not this iqtree2 repo): https://github.com/iqtree/iqtree3
You will need to compile the code there yourself. Hope that it's not difficult to compile.

Some more thoughts:

• I believe the rho parameter should not reach 1.0, but instead bounded between 0 and the number of variable sites in the alignment. So a more proper initial value can be half of the number of variable sites.
• If branch length is inferred by maximum likelihood under standard non-FunDi model, then the initial branch length of FunDi model should be smaller than that value (as we know that without FunDi, we tend to inflate the branch length). Therefore, a more sensible starting branch length should be smaller than that value, and we can bound the branch length between 0 and the branch length under non-FunDi model. Starting point can be half of this interval.

But please check this version and we can follow up later when we have some testing results....

[novigit]

We were able to compile this version and run it.

Setting --fundi-init-rho to 0.5 and 0.01 both led to much more reproducible, and sensible results! Also, decreasing the --fundi-epsilon from 1e-4 to 1e-6 also greatly enhanced to reproducibilty. For my dataset, 1e-6 seems to be a good enough value for all replicate runs to get the same results up to several digits after the dot. I didn't set the --fundi-init-branch, so it seems that setting a non-zero --fundi-init-rho is sufficient. I also tried to reproduce the "old" results, by setting the initial rho to 0.00, but the software wouldn't let me haha. I guess that is a good thing.

I think your thoughts make sense!

[bqminh]

awesome. Thanks for the confirmation!
Just for completeness (we also discussed): Another factor (but less impactful for your setting) for non-reproducibility is due to using 20 CPU cores -nt 20. I discussed this at length in #233.

[evolbeginner]

Hi Minh,

I have a relevant Q regarding the fundi model. I tried to perform analysis on a small, toy DNA alignment given a fixed topology. I used GTR or GTR+fundi with $\rho=0.5$. I found that both analyses gave exactly the same estimates of the subs model params and branch lengths (except for the internal branch that separate the indicated two subtrees). I wondered whether this looks correct. As i understand, under the fundi model, and under the FD hypothesis, the branch lengths of all subtrees as well as the params in the subs models would be estimated separately. I was using IQ-Tree v3.1.0. Did i understand it correctly? Any help would be highly appreciated.

GTR vs fundi_GTR

.iqtree summary

• Both runs used 1.aln with user tree 1.tre under GTR+F.
• GTR rates: A-C 2.1013, A-G 1.6951, A-T 1.0702, C-G 1.2558, C-T 1.0285, G-T 1.0000.
• Base frequencies: A 0.2375, C 0.245, G 0.235, T 0.2825.
• difference in bl of the .treefile file: internal branch length 0.5739104576 vs 0.0339178381 (see below, last number in the nwk formatted tree).

.treefile comparison

• iqtree-GTR.treefile: (t1:0.1726716596,t2:0.1615619893,(t3:0.7670966642,t4:0.48230 [iqtree-fundi.tar.gz](https://github.com/user-attachments/files/29283344/iqtree-fundi.tar.gz) 63580):0.5739104576);
• iqtree-fundi.treefile: (t1:0.1726716596,t2:0.1615619893,(t3:0.7670966642,t4:0.4823063580):0.0339178381);

cmd

• GTR iqtree -s 1.aln -m GTR -pre iqtree-GTR -te 1.tre -redo -wsl
• fundi: iqtree -s 1.aln -m GTR --fundi t1,t2,0.5 -pre iqtree-fundi -te 1.tre --fundi-init-rho 0.3 -redo --fundi-epsilon 1e-8 --fundi-init-branch 0.1 -wsl

iqtree-GTR.tar.gz
iqtree-fundi.tar.gz
1.aln.txt

best,
sishuo

[bqminh]

Right now the Fundi model implementation is based on the parameter estimates of the "main" model (GTR in your case). What iqtree does, is that at the end, when all parameters (incl. branch lengths, subst. model parameters) are estimated, it will switch to the Fundi model, and only either calculate the Fundi likelihood of a fixed rho (0.05 in your case) or estimate rho (which can be switched on by replacing 0.5 with the word "estimate").

There were some discussions of iteratively then reestimating other parameters, but it's not easy in terms of implementation. So this is not done yet.

In short, the behaviour you reported is expected.