dataesr/scientific-paragraphs-categorization huggingface dataset:
In all above analyses, we considered AS rates measured per intron, and not per gene. Yet, the trait under selection is the per-gene error rate, which depends not only on the error rate per intron, but also on the number of introns per gene. Given that intron density varies widely across clades (from 2. 8 introns per gene in diptera to 8. 4 introns per gene in vertebrates; Supplementary Tab. 1), the correlations reported above between AS rates and N e may undervalue the predictive power of the drift-barrier model. The RNA-seq datasets that we analyzed consist of short-read sequences, which do not allow a direct quantification of the per-gene AS rate. We therefore indirectly estimated the per-gene AS rate in each species, based on the per-intron AS rate and on the number of introns per gene (Materials & Methods).
Genomic and transcriptomic data collection To analyze AS rate variation across metazoans, three types of information are required: transcriptome sequencing (RNA-seq) datasets, genome assemblies, and gene annotations. To obtain this data, we first queried the Short Read Archive database  (Leinonen et al. , 2011)  to extract publicly available RNA-seq datasets. We also queried the NCBI Genomes database (NCBI Resource Coordinators, 2018) to retrieve genomic sequences and annotations. When this project was initiated, the vast majority of metazoans represented in this database corresponded to vertebrates or insects. We therefore decided to focus our analyses on these two clades (n=69 species).
Life history traits We used various life history traits to approximate the effective population size of each species. For vertebrates species we considered the maximum lifespan (i. e. from birth to death) and body length referenced. For insects we took the maximum lifespan and body length of the imago. For eusocial insects and the eusocial mammal Heterocephalus glaber, the selected values correspond to the queens. The sources from which the lifespan and the body length information was taken are listed in data/Data9 supp. pdf in the Zenodo repository (see Data and code availability).
, 2015). This dataset included 80,868,061 SNPs that were genotyped in 2,504 individuals. For each minor intron sharing one boundary with a major intron, we computed the number of SNPs that occur at their respective splice sites: at their shared boundary, and at the major intron and minor introns specific boundaries. We focused our study on minor introns that have their specific boundary folding in the exons adjacent to the major intron or in the major intron. As a control, for each minor intron, we searched for one GT and one AG dinucleotides in the interval between 20 and 60 bp with respect to the major splice site, in the neighboring exon and in the intron, and computed the number of SNPs that occur on these sites. We searched for control AG dinucleotides in the 5' of the intron and for GT dinucleotides in the 3' of the intron, to avoid studying sites that might correspond to minor splice sites.
Acknowledgements We thank  Loïc Guille  for his contribution to an initial pilot study,  Tristan Lefébure  for insightful discussions and  Laurent Guéguen  for his help on dN/dS analyses. Computational analyses were performed using the computing facilities of the CC LBBE/PRABI and the Core Cluster of the  Institut Français de Bioinformatique (IFB)  ( ANR-11-INBS-0013 ).
The observed number of amino-acid changes is estimated using ancestral character reconstruction, and their expected number in a null model using computer simulations. In the following sections, we describe this approach that is implemented in a software and web service named ConDor (Convergence Detector; condor. pasteur. cloud). Its performance is assessed on HIV-like simulated datasets, on a real HIV reverse transcriptase dataset and on a fish rhodopsin dataset. We notably compare its performance in retrieving drug resistance mutations in HIV, with the results of a standard positive selection-based approach.
Simulation-based approach Our method identifies amino-acid mutations that emerged several times in independent lineages and occurred significantly more frequently than expected under a null model of evolution. The workflow of the method is presented in Figure  1. It is made up of four main steps: (1) estimate the parameters of the null model from the real data (phylogenetic tree, parameters of the substitution model, site-wise evolutionary rates, etc. ); (2) infer ancestral amino acids and count the number of observed emergence events of mutations (EEMs) for every position and amino acid of interest in the real data; (3) simulate new datasets under the inferred null model and count simulated EEMs;  (4)  compare the observed and the simulated numbers of EEMs, and then determine which mutations occurred significantly more often than expected by chance, assuming the null model.
Taking only the amino acid with the highest posterior tends to bias the simulations and yields poorer results (not shown), especially if the reconstruction is uncertain (e. g. , two amino acids with posteriors of 0. 55 and 0. 45). Simulations are carried out along the inferred tree, and then we count the simulated numbers of EEMs (10,000 values per position and per studied AA) using the algorithm presented below. For example, let us consider the mutation M41L from our real HIV dataset, where at position 41, a Methionine (M) is substituted by a Leucine (L) in 211 sequences. The observed number of EEMs towards L is 47, which is smaller than 211 as in some subtrees all tips have L, corresponding to only 1 EEM. This number is compared to the distribution of the number of EEMs towards L, starting from an M at the tree root every time (no ambiguity in ACR), among 10,000 simulations in the null model; this The method takes as input an amino-acid alignment, which is used for inference of the null model (phylogenetic tree, substitution model and its parameters, site-wise evolutionary rates, etc.
Creation of a synthetic HIV-like dataset To our knowledge, there is no convergent evolution model allowing simulating thousands of sequences without prior knowledge of the phenotype or environmental constraints. We therefore created our own convergent dataset inspired from a real case of convergence. Drug resistance mutations (DRMs) occur independently in patients receiving a drug treatment and thus are a perfect example of evolutionary convergence. In HIV, they are well characterized and studied since their emergence can lead to treatment failure and transmission of resistant viruses. They are mostly found in proteins targeted by the antiretroviral treatment: the protease, the reverse transcriptase and the integrase. The list of known DRMs affecting these proteins is publicly available at https://hivdb.
stanford. edu/. DRMs are written as "XposY", with X the ancestral amino acid, "pos" the position of the substitution in the protein alignment, with numbering based on the reference sequence HXB2, and Y the mutated amino acid. We extracted positions and sequences with DRMs from a real HIV polymerase amino-acid dataset (described in Material and Methods) and replaced the corresponding residues with gaps in the multiple sequence alignment (MSA). DRMs were retrieved from the Stanford University Drug resistance database "Essential DRM Data" lists (https://hivdb. stanford. edu/pages/poc. html). We then reconstructed a tree, estimated the substitution model parameters, and inferred the rates of evolution per position (see Material and Methods for details). This tree represents the evolutionary relationships between the sequences in the real data and its structure is not affected by the DRMs.
Such procedure is standard when inferring trees from HIV sequences, to avoid convergence perturbations of the inferred tree. We then simulated the evolution of HXB2 (GenBank TM accession number K03455) reference sequence (reverse transcriptase only) along this tree. We performed the simulations 5 times for robustness purposes, resulting in 5 MSAs without convergence. DRMs were then manually added in the sequences and positions where they were found in the real MSA. For example, mutation M41L was found in the real data in 211 sequences, so we replaced in the corresponding sequences of our simulated alignment, any amino acid found at position 41 (which turned out to all be methionine (M), the same as HXB2 at that position) by a leucine (L). Since we used the real data tree to create the synthetic dataset, we did not randomly place DRMs, which would make the detection task easy.
Methionine and leucine are closely related, with a genetic barrier of 1 (M and L can exchange through a single nucleotide substitution), and the detection is much more difficult than between highly different amino acids (e. g. , D and W with a genetic barrier of 3). As M41L is at the same position and in the same tip sequences in our synthetic dataset as in real data, detecting this convergence in the synthetic and real datasets should be of similar difficulty. This guarantees a certain realism of our simulated data. We implemented this insertion procedure for the 37 most common DRMs of our real dataset (i. e. , present in ≥ 12 sequences and > 2 independent clades). The fiveresulting synthetic MSAs thus have no convergent events, except the "realistic" added DRMs.
Results We applied our approach on three datasets for which we knew a priori mutations due to convergent evolution: (1) an HIV-like synthetic dataset with "realistic" added DRMs; (2) a real HIV dataset of reverse transcriptase with 20% sequences with DRMs; and (3) a real dataset of fish rhodopsin, a light-sensitive receptor protein that is highly conserved but known to vary at certain positions among species depending on their environment. On simulated data, we know exactly which mutations are true convergent events or not. With real data, even if we know certain convergent events, some other mutations are likely convergent, but are unknown. In other words, we will be able to assess the rate of false positives with simulated data, but not with real data, where the method sensitivity in detecting the known convergence events will be the main criterion.
Synthetic HIV-like alignment The dataset consists in 5 MSAs of 3,551 sequences and 250 amino acids each, mimicking HIV reverse transcriptase and simulated with HIVb model of evolution  (Nickle et al. 2007). In total, 37 DRMs were placed in each of the MSAs at 27 distinct positions (see above). These DRMs are found in at least 12 sequences. The most common one, M184V, is found in 273 (8%) sequences. However, 19 DRMs are found in less than 1% of the sequences (i. e. , in 12 to 33 sequences) so they are expected to be difficult to detect.
Performance of the method on the detection of true convergence The model inferred from the datasets by ModelFinder  (Kalyaanamoorthy et al. 2017)  was HIVb, which is the model we used for generating them. Thus, the whole analysis (tree reconstruction, ACR and simulations) was first achieved with HIVb, which is the true model of evolution. In a second stage, we also used JTT  (Jones et al. 1992) , to study the impact of model misspecification. We tested on average 441 mutations per dataset, the ones present in at least 12 sequences and with more than 2 EEMs. Using HIVb, on average 27. 4 mutations are found to be convergent by our method, 26 of which are true DRMs (out of 37 added DRMs). The numbers of EEMs for the DRMs range from 9 (K101P) to 225 (M184V). We detect DRMs with a higher number of EEMs better, and especially those with more than 30 EEMs that we detect 90% as convergent, as illustrated in Figure  3a.
If there are several DRMs at one position, we often only detect the most frequent one(s). For example, at position 219 (Fig3a, bold characters) we detect mutations towards Q and E but not N. Similarly, at position 215, we do not detect mutations towards S and D. However, we detect mutation T215C although there are fewer EEMs towards C. This is explained by the low substitution rate between T and C (BLOSUM62 score = -1), and thus few changes are expected from T to C; 2 EEMs between T and C at position 215 are observed on average over the 10,000 null-model simulations, while 17 are found in the synthetic, convergent MSAs. If we focus on detecting positions with convergence (e. g. , position 219) rather than DRMs (e. g. , K219Q, K219E, K219N, etc. ) accuracy increases and we detect on average 22 of the 27 convergent positions for all datasets, while the number of false positives remains equal to 1.
4 on average (Tab. 1). Table  1 : Method accuracy with synthetic data. We display the number of mutations and positions (1) detected and (2) tested but not detected on the 5 synthetic HIV-like MSAs, analyzed with HIVb and JTT substitution models. We report the average for the 5 datasets and the standard deviation between parentheses. True positives are at the intersection between detected and DRMs, and false positives at the intersection between detected and non-DRMs (i. e. , mutations resulting from the evolution of the root sequence under the null model). In these experiments, we tested all mutations (regardless of their DRM status) exhibiting more than 2 EEMs and found in at least 12 sequences. Mutations in bold characters are discussed in text. The evolutionary rate (normalized between 0 and 1) of each mutation position is reported on the y-axis.
The dashed vertical line represents the limit of 30 EEMs, all events on the left have less than 30 EEMS. In the MSA represented in Figure  3a , we find two false positives, one of which exhibits a very high evolutionary rate. The same results are observed in the other MSAs: as the evolutionary rate increases, more changes are observed at the given position and thus more variability and uncertainty in the simulations. Thus, very fast positions can bias convergence detection and lead to the detection of false positives. As expected, we observe very few false positives when analyzing our synthetic datasets with the true model of evolution. Since we never have the true model of evolution with real data, we tested the effect of model violation on the synthetic dataset. We run the whole analysis with JTT, for tree reconstruction, ACR and simulations, instead of letting the workflow infer and use the best model of evolution (here HIVb).
The strongest effect can be seen on the number of false positives, which increases from 1. 4 to 17. 2 (Tab. 1). Compared to the number of negatives, this remains very low, with 4% of falsely detected random mutations (i. e. , mutations resulting from the evolution under the null model) among more than 400. Moreover, as we see Figure  3b , the 6 events with highest evolutionary rate are false positives. The same trend can be observed in the other MSAs, where the mutations at positions with highest rates are always false positives. Since the substitution model does not correspond exactly to the synthetic datasets being analyzed, we tend to detect more mutations as convergent, but this does not impact the detection of DRMs. DRM detection remains sensitive and is robust to model violation, as we still detect more than 25 (among 37) DRMs, and more than 22 DRM positions (among 27), that is, nearly the same true positive fractions as with the true HIVb substitution model.
Real HIV dataset This dataset consists in a MSA of truncated polymerase from HIV-1 subtype B. It was retrieved from the paper by  Lemoine et al (2018). After removal of recombinant sequences, it contains 3,546 sequences of 1,043 nucleotide positions that were translated into 347 amino acids. Among these 347 amino acids, 250 are on the reverse transcriptase and are analyzed here. Slightly more than 20% of the sequences have at least one known DRM (https://hivdb. stanford. edu/pages/poc. html) and, on average, the DRMs are found in ~11 sequences. The most common one, M184V, is found in 273 sequences. There are 37 DRMs present in at least 12 sequences, corresponding to those used to generate the synthetic datasets. They are distributed on 27 positions. We focused on these 37 DRMs to assess the performance of our approach, but, as already explained, we expect to detect other mutations, some being truly convergent but unknown, and some corresponding to false positives, likely located on fast positions, due to model misspecification.
The evolutionary model selected by ModelFinder  (Kalyaanamoorthy et al. 2017 ) on this dataset is HIVb, with 'freerates' rates-across-site model and 9 rate categories. We tested 255 mutations in total: those present in at least 12 sequences and showing more than 2 EEMS. Among these, we detected 74 convergent events, after applying the Benjamini-Hochberg correction  (Benjamini and Hochberg 1995)  for multiple testing (non-corrected p-value threshold of 4e-4, corresponding to a corrected alpha level of 5%). Among these detections, 20 are DRMs, which represents 54% of true positives and is a higher proportion than what is expected by chance (Fisher's exact test p-value = 4. 8e-4). The non-DRM detected events correspond to 11 mutations on fast evolving positions (Fig. 4a ), which are likely false positives, plus other events for which we cannot conclude using simple arguments (but see below).
Among non-fast positions, only a small fraction (9/36) of the non-DRM events detected with ConDor were thus found at positions with positive selection. However, while reviewing the proposed convergent mutations found by ConDor, we could identify 20 mutations on 18 (among 36) positions previously described in the literature as potential accessory mutations, correlated with virological failure (see Supplementary Table  S2 ). Moreover, as we previously noticed with synthetic datasets, our method tends to detect events at fast positions. The same holds with positive selection, thus leading us to consider the 5% of the fastest positions (calculated on the whole dataset) as false positives. This is comforted by several studies indicating that the main polymorphisms in HIV-1 are found on the fast positions  (Bao et al.
Fish Rhodopsin Rhodopsin is a photosensitive protein pigment responsible for the eye's sensitivity to light. It is found in many vertebrates and has been shown to be under positive selection among species that evolve in different environments  (Spady et al. 2005). Depending on the habitat and the amount of light available, different amino acids are observed at the same positions, which result in variations in rhodopsin structure and different maximum wavelength absorption. Certain substitutions corresponding to these amino-acid changes have been described as resulting from convergent evolution. In particular D83N, E122Q, F261Y and A292S (using similar substitution encoding as with HIV) occurred several times independently  (Yokoyama 2008). The dataset we used comes from a study in which the authors characterized substitution F261Y as convergent in fish rhodopsin, as a possible result from a transition from marine to brackish or freshwater environments  (Hill et al.
2019). It contains an alignment of 2,047 sequences with 308 amino-acid positions. The sequences have been classified by the authors in two groups: species found only in marine water and species that can evolve (exclusively or not) in brackish or fresh water. Species annotated within the habitats brackish or fresh water can therefore also be found in marine water. The neutral model inferred by ModelFinder on this dataset was 'MtZoa' and 'freerates' with 8 rate categories. The reconstructed tree is well supported with three quarters of the bootstrap supports above 70%. We tested 355 substitutions (present in at least 12 sequences and showing more than 2 EEMs) with our method and 55 were retained as possibly convergent (15%). They are distributed over 49 positions out of the 136 tested positions, which means that 36% of the tested positions are detected as convergent (16% considering all alignment positions).
These 9 reversions most likely are true convergence events. We tested the correlation between the 55 predicted convergent events and the habitat and found a significant dependence for 32 of them (9 on fast positions), meaning that, again, a substantial fraction of these detections most likely corresponds to true convergence. This was confirmed by a literature review (see Supplementary Table  S4  for details) showing that ConDor detected several interesting mutations such as T166S responsible for blue-shifting absorption  (Luk et al. 2016 ), S299A found mainly in bathypelagic species  (Dungan and Chang 2017) , or A124G causing a spectral shift towards red wavelengths  (Van Nynatten et al. 2021 ). To confirm our results and compare our method, we also ran PCOC  (Rey et al 2018)  on this dataset, considering that the convergent clades were those annotated with the brackish or freshwater environment.
Implementation and availability Our method is implemented in a workflow, named ConDor, which is accessible through a web service: http://condor. pasteur. cloud. ConDor inputs consist of: (1) a protein alignment in fasta format; (2) a file containing outgroup sequence identifiers; and (3) a newick file with a tree (rooted or not) whose tips are the same sequences. The workflow is executed on the Institut Pasteur cluster and takes ~160 minutes for a dataset of ~2,050 sequences and ~300 positions. This time corresponds to ~140 minutes for the model inference and ~20 minutes for the convergence detection, using 10,000 simulations per position and target amino acid. The output consists of all the tested mutations with several statistics such as their p-value and if they passed or not the threshold to be considered as convergent mutations.
Discussion In this work, we have developed the ConDor approach, which makes it possible to detect evolutionary convergence at the resolution of a mutation, without prior knowledge of constraints and phenotype. We could retrieve more than 70% of the positions with DRMs on a real HIV dataset, and ~55% of all DRMs present in the sequences. On synthetic datasets mimicking the evolution of HIV, we showed that the detection power depended on the number of emergence-events of mutations and the exchangeability between amino acids. Even though ConDor was primarily developed for the analysis of viral datasets, it could detect several convergent mutations involved in the change in absorption wavelength in a fish rhodopsin dataset. These results confirm that our method detects realistic convergent evolution signal and could be applied to a broad range of organisms.
We tested the robustness of ConDor to model violation by using JTT  (Jones et al. 1992 ) instead of HIVb  (Nickle et al. 2007 ) as a neutral model of evolution for the study of synthetic HIV-like datasets. In doing so, the sensitivity remained high, and we still detected ~70% of positions with DRMs. However, the number of false positives increased, and thus we expect false positives with real data. The number of false positives can be partly explained as ConDor tends to be biased towards mutations found at positions with very high evolutionary rates. Indeed, our method relies on how realistic simulations are as a null model. Our results show that for most positions and most mutations, we are close to what is observed in real data, and our simulations represent a satisfactory null model. However, at certain fast positions we observed that simulations tend to differ from the real data, which resulted in an increased rate of false positives.