Suggestion for selection analysis using multiple-hits + alignment errors and troubleshooting GARD

[desmodus1984]

Hi,

I am interested in detecting genes under selection and I recently found out about the alignment error masking and the multi-hit options.
My goal is to identify potential genes under selection in a couple of sister species. I would appreciate if you could help with comments and advice on how to conduct the analysis.

First, I identified single-copy orthogroups with Orthofinder from a set of 13 genomes, then my analysis steps are:

1. Retrieve the transcripts from a database, and align them with MACSE.
2. Filter the alignment with Clipkit using the codon model in default mode, and then trim further using a restrictive gap threshold ( -g 0.2).
3. Remove the stop codons with hyphy cln
4. Lastly, run HYPHYMPI absrel with --srv Yes

My question is where to do the alignment filtering (BUSTED-E) , and whether the multiple-hit should be done before or after alignment filtering because BUSTED-E can also be done with the --multiple-hits option.

For instance, I ran the multi-hit on some aligments, and I got for instance:

### Summary of rate estimates and significance testing
|                Model                 |   Log-L    |   omega    | 2-hit rate |              p-value               | 3-hit rate |              p-value               |
|:------------------------------------:|:----------:|:----------:|:----------:|:----------------------------------:|:----------:|:----------------------------------:|
|            Standard MG94             |-10977.33...|    0.2316  |    N/A     |                N/A                 |    N/A     |                N/A                 |
|        Standard MG94 + 2 hits        |-10952.29...|    0.1801  |    0.4034  |       0.0000 (2-hit rate = 0)      |    N/A     |                N/A                 |
|     Standard MG94 + 2 or 3 hits      |-10923.71...|    0.1723  |    0.1899  |      0.0000 (2&3-hit rates = 0)    |    1.8212  |      0.0000 (3-hit rate(s) = 0)    |

### 5 individual sites which showed sufficiently strong preference for multiple-hit models

|   Site   |       ER (2 vs 1)       |       ER (3 vs 2)       |                       Substitutions                        |
|:--------:|:-----------------------:|:-----------------------:|:----------------------------------------------------------:|
|   695    |          63.5193        |         227.6457        |                     TCA->AGC(1)CTC(1)                      |
|   698    |         124.6144        |         156.6477        |                     AGC->CAT(1)TCT(1)                      |
|   1096   |           4.6150        |          24.3391        |                  TCG->AGC(1)CAA(1)TCT(1)                   |
|   1097   |           5.3539        |           9.7650        |                     CAA->GCC(1)TTG(1)                      |
|   1106   |           3.1503        |           7.5585        |                     TCT->AAG(1)CCC(1)                      |

and the alignments stats are:

>Loaded a multiple sequence alignment with **13** sequences, **1112** codons, and **1** partitions

According to the fmm analysis, the model with the highest logL is the one with either 2 or 3 hits, but out of 1112 codons, only 5 sites have preference for multiple hits. Is this enough evidence (enough sites) to use BUSTED-E or absrel with the 3-hits-multiple-hits model?

For instance, another alignment:

>Loaded a multiple sequence alignment with **13** sequences, **268** codons, and **1** partitions

### Summary of rate estimates and significance testing
|                Model                 |   Log-L    |   omega    | 2-hit rate |              p-value               | 3-hit rate |              p-value               |
|:------------------------------------:|:----------:|:----------:|:----------:|:----------------------------------:|:----------:|:----------------------------------:|
|            Standard MG94             |  -3175.17  |    0.4972  |    N/A     |                N/A                 |    N/A     |                N/A                 |
|        Standard MG94 + 2 hits        |  -3089.86  |    0.3681  |    1.1937  |       0.0000 (2-hit rate = 0)      |    N/A     |                N/A                 |
|     Standard MG94 + 2 or 3 hits      |  -3051.41  |    0.3651  |    0.7279  |      0.0000 (2&3-hit rates = 0)    |    4.2756  |      0.0000 (3-hit rate(s) = 0)    |

### 19 individual sites which showed sufficiently strong preference for multiple-hit models

|   Site   |       ER (2 vs 1)       |       ER (3 vs 2)       |                       Substitutions                        |
|:--------:|:-----------------------:|:-----------------------:|:----------------------------------------------------------:|
|    4     |          15.0814        |           6.8236        |                     CAA->AGG(1)CTG(1)                      |
|    18    |          20.5916        |           0.3047        |               ACT->AGG(1)GCA(1), GCT->ACT(1)               |
|    21    |          72.7447        |           7.7902        |                     ACA->CGA(1)GAG(1)                      |
|   105    |          53.3908        |           0.4668        |               CCG->CCC(1)GGG(1), GGG->CCG(1)               |
|   214    |          10.0952        |           0.5638        |                  CCC->CCT(1), CCT->CAG(1)                  |
|   215    |          23.3547        |          10.7724        |            CTC->AGC(1), TCT->CCT(1)CTC(1)TCG(1)            |
|   216    |          45.6818        |           1.7204        |                  CCA->TGG(1), CTG->CCA(1)                  |
|   217    |          44.1538        |           7.6918        |               GCA->TGC(1), GTG->GCA(1)GCG(1)               |
|   221    |           9.5131        |           2.3019        |                     GCA->CGG(1)GCG(1)                      |
|   223    |          10.7141        |          16.6627        |                     AGT->GGT(1)GTA(1)                      |
|   248    |          13.0184        |          13.1496        |                        GCC->AGA(1)                         |
|   259    |           8.8395        |           0.5507        |                        GAA->CAC(1)                         |
|   263    |           9.7945        |           4.6021        |        GCC->CAC(1)GCG(1), GCG->GTG(1), GTG->AGC(1)         |
|   268    |          26.8636        |          20.4134        |                        CTC->ACT(1)                         |
|    20    |           4.6980        |          13.0971        |                  GAG->AAG(1)AGA(1)GAA(1)                   |
|   182    |           5.2119        |           8.2645        |                        TGC->CAG(1)                         |
|   251    |           6.3157        |          13.5233        |                  ACC->ACA(1)AGC(1)GAG(1)                   |
|   258    |           4.2727        |           8.3152        |            CCG->CCA(1)GAA(1), CCT->CCC(1)CCG(1)            |
|   267    |           4.9757        |           7.8712        |                  CTC->AGG(1)CAC(1)CTT(1)                   |

The results were more convincing about using multi-hits, 19 sites out of 268 codons.

My concern is if multiple-hits sites are found, but then BUSTED-E marks them as erroneous, making the usage of multiple-hit model unnecessary or wrong.

I'd appreciate if you could help me decide whether I should use BUSTED-E + multi-hit detection.

Also, I have Hyphy installed in a server, and I am trying to run GARD to detect recombination and I get an error.
I first ran:

for f in `cat test_ortho.tsv`; do echo -e "\n$f"; hyphy gard --type nucleotide --alignment "${f}.CDS.NT.trimmed.g0.2.STOP.fa" --model JTT --mode Normal --rv None --rate-classes 4 --output ${f}.gard.json ; echo -e "Done $f\n" ; done

and I got the error:

### Fitting the baseline (single-partition; no breakpoints) model
* Log(L) = -8770.42, AIC-c = 17603.40 (31 estimated parameters)

### Performing an exhaustive single breakpoint analysis
[GARD] Breakpoint          6 of 576. Best cAIC =   17603.3967 with no breakpoints.Error:
Internal error
Internal error in ComputeBranchCache (branch gard.tree.part_0.Musnig, eval #3993 ) reversible model cached likelihood =   -89.57536821397055, directly computed likelihood =   -89.57536805683237. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.

Function call stack
1 :  [namespace = LYWdsGte] Optimize(mles, ^lf_id);

        Keyword arguments:
                {
                 "model":"JTT",
                 "rate-classes":"4"
                }
-------
2 :  [namespace = ZDRcMlPV] res=estimators.FitExistingLF(&likelihoodFunction,modelObjects);

        Keyword arguments:
                {
                 "model":"JTT",
                 "rate-classes":"4"
                }
-------
3 :  [namespace = fE_Rixoy] fit=gard.fitPartitionedModel(breakPoints,model,initialValues,null,FALSE);

        Keyword arguments:
                {
                 "model":"JTT",
                 "rate-classes":"4"
                }
-------

then, I tried

for f in `cat test_ortho.tsv`; do echo -e "\n$f"; ENV="TOLERATE_NUMERICAL_ERRORS=1;" hyphy gard --type nucleotide --alignment "${f}.CDS.NT.trimmed.g0.2.STOP.fa" --model JTT --mode Normal --rv None --rate-classes 4 --output ${f}.gard.json ; echo -e "Done $f\n" ; done

and I still got the error.
For example, I ran the multiple-hits with this script without any problem:

for f in `cat test_ortho.tsv`; do echo -e "\n$f"; hyphy fmm --alignment ${f}.CDS.NT.trimmed.g0.2.STOP.fa --tree tree.nwk --srv Yes ; echo -e "Done $f\n" ; done

I would be grateful of any suggestion on how or where to add the code to correct for numerical instability.

Thanks

[spond]

Dear @desmodus1984,

hyphy ENV argument passing is a bit funny. It does not read those from the shell environment, but rather as a a command line argument.

So instead of

ENV="TOLERATE_NUMERICAL_ERRORS=1;" hyphy gard ....

do

hyphy gard ... ENV="TOLERATE_NUMERICAL_ERRORS=1;"

and the error will become a warning.

Best,
Sergei

[spond]

Dear @desmodus1984,

Yes, you should use BUSTED-E if unsure.

hyphy busted --error-sink Yes --rates 2 --syn-rates 2 --starting-points 5 --multuple-hits Double+Triple ...

1. You don't have to prescreen the alignment before running BUSTED-E -- it's a one stop shop.
2. If you want to check the influence of "errors" run busted with and without the --error-sink Yes flag. If you get selection without, and no selection with, then you have a good indication that errors are to blame.
3. For 13 sequences uses fewer than the 3 standard rates (2 here)
4. The tell-tale sign of errors are multi-nucleotide substitution in consecutive sites (local mis-alignment). In your second example even fmm (which does not allow branch-by-branch variation) picks up 214-217 (and 221, 223). My suspicion is one of your sequences at those positions will "stick out" even if you look at it by eye.

Best,
Sergei

[desmodus1984]

Good morning,

Thanks for the reply.
I tried gard as you mentioned, it ran and I got several "warnings".

>Loaded a nucleotide multiple sequence alignment with **13** sequences, **1041** sites (207 of which are variable) from `/data/common/data/genomes/LW_13_sel/OG0007202.CDS.NT.trimmed.g0.2.STOP.fa`
>Minimum size of a partition is set to be 23 sites


### Fitting the baseline (single-partition; no breakpoints) model
* Log(L) = -2909.38, AIC-c =  5882.74 (31 estimated parameters)

### Performing an exhaustive single breakpoint analysis
[GARD] Breakpoint          9 of 206. Best cAIC =    5882.7356 with no breakpoints.
***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_0.Halgry, eval #8102 ) reversible model cached likelihood =   -147.5932649632766, directly computed likelihood =   -147.5932621242337. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.
[GARD] Breakpoint         10 of 206. Best cAIC =    5882.7356 with no breakpoints.
***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_0.Hydlep, eval #9552 ) reversible model cached likelihood =   -184.6199094887444, directly computed likelihood =   -184.6199066495592. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.
[GARD] Breakpoint         15 of 206. Best cAIC =    5882.7356 with no breakpoints.
***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_0.Lepwed, eval #15171 ) reversible model cached likelihood =   -261.1710520493303, directly computed likelihood =   -261.1710492096663. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.
[GARD] Breakpoint        174 of 206. Best cAIC =    5882.7356 with no breakpoints.
***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_1.Mirang, eval #202498 ) reversible model cached likelihood =   -450.9695115380054, directly computed likelihood =   -450.9695086986148. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.
[GARD] Breakpoint        178 of 206. Best cAIC =    5882.7356 with no breakpoints.
***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_1.Mirang, eval #207307 ) reversible model cached likelihood =    -471.499913787015, directly computed likelihood =    -471.499910946857. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.
[GARD] Breakpoint        182 of 206. Best cAIC =    5882.7356 with no breakpoints.
***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_1.Mirleo, eval #212323 ) reversible model cached likelihood =   -418.9256432612639, directly computed likelihood =   -418.9256404226348. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.
[GARD] Breakpoint        186 of 206. Best cAIC =    5882.7356 with no breakpoints.
...
...
...
### Performing multi breakpoint analysis using a genetic algorithm

***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_2.Hydlep, eval #275115 ) reversible model cached likelihood =   -207.2163729833418, directly computed likelihood =   -207.2163701439508. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.

***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_2.Mirleo, eval #348713 ) reversible model cached likelihood =   -394.0054783254287, directly computed likelihood =   -394.0054754860377. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.

***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_1.Eumjub, eval #351466 ) reversible model cached likelihood =   -1552.567301233025, directly computed likelihood =   -1552.533169080638. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.

***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_1.Mirang, eval #354406 ) reversible model cached likelihood =   -389.4131826151919, directly computed likelihood =   -389.4131797758009. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.
[GARD]   2 breakpoints [ generation      1, total models      242,    0% converged, 34.57/sec]. Min (c-AIC) =    5910.4783 [no improvement], breakpoints at 368, 875
***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_2.Mirleo, eval #425769 ) reversible model cached likelihood =   -425.7408048167185, directly computed likelihood =   -425.7408019773276. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.
...
...
***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_1.Neosch, eval #13220224 ) reversible model cached likelihood =   -464.9345911600223, directly computed likelihood =   -464.9020619052868. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.
[GARD]   2 breakpoints [ generation    114, total models     5157,   99% converged, 22.72/sec]. Min (c-AIC) =    5890.5055 [no improvement], breakpoints at 380, 591
***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_1.Odoros, eval #13311428 ) reversible model cached likelihood =   -353.7671375160788, directly computed likelihood =   -353.7671346766879. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.
...
...
***WARNING***
Internal error in ComputeBranchCache (branch gard.tree.part_2.Mirang, eval #13336920 ) reversible model cached likelihood =   -259.2339845778273, directly computed likelihood =   -259.2339817379635. This is most likely because a non-reversible model was incorrectly auto-detected (or specified by the model file in environment variables; for smaller errors, this could be due to numerical instability of calculations for larger alignments). To treat numerical errors as warnings, please specify "ENV=TOLERATE_NUMERICAL_ERRORS=1;" as the command line argument. This often resolves the issue, which is indicative of numerical instability.
Done with 2 breakpoint analysis.
    Best break point locations: 380, 591
    c-AIC = 5890.50549596655

Done OG0007202

Are these "warning/errors" "signal" of bad alignments?
I extracted and aligned the longest isoform per gene with AGAT (agat_sp_keep_longest_isoform.pl) and MACSE, respectively. Then, I inspected a few alignments and some sequences were either shorter or longer than the others, which I am worried it might affect the selection analysis due to gaps.

You said to go and try BUSTED, BUSTED with error-sink, and GART, and I wanted to inquire if alignments like those might need perhaps a "quick trimming" before doing any analysis.
Also, after doing GARD on the trimmed alignments, and I checked the output, and the trees for each partition were different compared to the original with an inverted branching pattern, with the most derived species at the base of the tree.

Original tree:

Can I still specify the analysis to use my tree despite the GARD output including different ones?

Finally, I would like to ask about the multi-hit model selection.
Let's say that I ran GARD on an alignment, and it splitted in two blocks; will the multi-hit test work?
Maybe I am overthinking; - it fmm indicated double-hits just for only the first block, I would have to then parse and split the blocks in the GARD output and separately indicate the multiple-hits parameter for each block for BUSTED/aBSREL analysis?

I would like to run all the recommended tests (GARD/multi-hits), but I feel overwhelmed - new to bioinformatics, since I have 8906 single-copy orthologs alignments.

Thank you very much;