dataesr/scientific-paragraphs-categorization huggingface dataset:
The study showed that some of the mutations (RT:D67N, RT:K70R, RT:M184V, RT:K219Q) decreased the fitness, one (PR:L90M) seemed to increase the fitness, while the others did not have a significant effect. The models above, and more generally, the family of multi-type birth-death models  [7]  with a Bayesian birth-death skyline plot (allowing the parameters to change in a piecewise constant manner)  [12]  they belong to, define ODEs for fine-tuned parameter estimation. However, their complexity prevents resolving them analytically. Numerical solutions of the ODEs, on the other hand, have long computational times, which prevents the application of these models to larger datasets (dozens of thousands of sequences), while larger datasets are desirable for more accurate parameter estimation. A compromise between the model complexity and computational speed when applied to large datasets needs to be found.
Moreover, we account for the changes in treatment policies by allowing for separate ACRs for different time intervals (e. g. , before and after the first DRM-provoking ARV introduction). Once the reconstruction is performed, we visualize the results with PastML and calculate various statistics for transmission patterns. We apply our approach to analyze the patterns of DRM emergence, transmission, and loss in HIV-1-infected individuals in the UK, using sequences and metadata from the UK HIV Drug Resistance Database  [14].
Materials and Methods The UK HIV Drug Resistance Database provides HIV protease (PR) and reverse transcriptase (RT) sequences extracted during resistance tests and their corresponding metadata (e. g. , treatment status of the patient before the test: treatment-experienced, -naive, or unknown; date of the test). These sequences were obtained by Sanger sequencing, which provides the consensus sequence of the predominant virus variants in a patient at the time of sampling. These sequences are subject to sampling uncertainty and potential sequencing errors, but the most discordant sequences were detected as outliers and removed from our analyses (see Transmission tree reconstruction). Furthermore, our dataset is so large (∼80,000 sequences) that the few remaining errors, which are unavoidable, most likely have a small impact on the overall results.
Polymorphic mutations can appear spontaneously in treatment-naive individuals and hence are an important source of phylogenetic signal and were kept. It should also be noted that our trees were reconstructed from DNA sequences, not proteins, which reduces the impact of convergent mutations at the amino acid level. We then dated each tree with LSD2  [21]  (v2. 3: github. com/tothuhien/lsd2/tree/v1. 4. 2. 2 accessed on 1 April 2023) using tip sampling dates, under strict molecular clock with outlier removal (see Appendix A for detailed settings). Sequences whose mutation rates were different from the mean rate for more than three standard deviations were considered as outliers.
Ancestral Character Reconstruction For each DRM (surveillance, polymorphic or accessory) listed in the Stanford HIV Drug Resistance Database [22] (hivdb. stanford. edu/pages/surveillance. html accessed on 10 May 2023), we extracted its presence/absence in the (unaligned) sequences of our datasets and the ARVs that can provoke it with Sierra, the Stanford Algorithm  [23]  web service. We then analyzed the DRMs that were found in at least 0. 5% of sequences (after filtering by patient and temporal outlier removal) of our dataset (either B or C, analyzed separately). Each DRM name (e. g. , RT:T215D) contains 2 pieces of information: the DRM position (e. g. , RT:T215) and its mutated amino acid, associated with resistance (e. g. , D). The DRM position contains the protein name: RT (reverse transcriptase) or PR (protease), the reference position of the amino acid (e.
g. , 215), and its wild-type amino acid (e. g. , T). We analyzed DRMs with prevalence > 0. 5% for each dataset. Each DRM position (e. g. , RT:V179 in the C dataset) was analyzed independently by reconstructing its states in the ancestral nodes, based on the tip states. The DRMs with prevalence > 0. 5% found in this position (e. g. , RT:V179D and RT:V179E in C) were analyzed together. For the majority of the DRM positions, only one DRM (i. e. , mutated amino acid, associated with resistance) with prevalence > 0. 5% was found (e. g. , PR:L90M for the position PR:L90). In the B dataset two positions contained several DRMs (> 0. 5%): RT:T215 (RT:T215D, RT:T215F, RT:T215S, and RT:T215Y) and RT:K219 (RT:K219E, RT:K219N, and RT:K219Q). In the C dataset, one position contained several DRMs: RT:V179 (RT:V179D and RT:V179E).
Possible states for ancestral character reconstruction (ACR) corresponded to DRM presence (i. e. , the resistant state) or absence (sensitive state) for DRM positions with only one DRM. For instance, for PR:L90M, the resistant state corresponds to the amino acid M, and the sensitive state to any other amino acid; in practice, the sensitive state is almost uniquely L. In the B dataset, 97. 79% of sequences have L at the position PR:90; 1. 93% have M; less than 0. 01% have W or F, and 0. 23% have an ambiguity at this position (so their initial state for ACR is unresolved between sensitive and resistant). For positions with several DRMs, the resistant state was split into all the possibilities (e. g. , D, F, S, or Y for RT:T215). For polymorphic mutations (e. g. , RT:S68G), ACR was performed on the corresponding (B or C) time-scaled tree with PastML (v1.
9. 40, MAP (maximum a posteriori) decision rule), without taking into account the year of ARV acceptance, as the these mutations could be present independently of ARVs. To reconstruct the ancestral character states for non-polymorphic DRMs, we used the procedure visualized in Figure  1a  (which we first proposed and applied to study HIV resistance patterns in Cuba in  [24] ). For each ARV, we extracted the dates of their acceptance with the Wikipedia python package (https://github. com/goldsmith/Wikipedia accessed on 1 April 2023). We cut the time-scaled tree at the earliest of the dates of acceptance of ARVs that can provoke the DRM (e. g. , for PR:L90M, saquinavir (SQV) was accepted in 1995). We hence obtained the pre-treatment-introduction tree and a forest of post-treatmentintroduction subtrees.
We count all the identified DRM loss events, all the identified (observed and hidden) TDR events, and only those ADR events that are not at the root of naive-only TDR clusters, as the latter happened in unsampled treatment-experienced patients (see Figure  2b, d ). Note that N resistant tips = N ADR + N TDR -N loss. For example, in Figure  2c , all the events correspond to observed tips, so we count one ADR, three TDR, and one DRM loss events: N resistant tips = 3 = 1 + 3 -1. Figure  2d  represents a more complex case: we count one hidden TDR event (as it led to the resistance status of one of the observed tips) and three observed TDR events (leading to resistance statuses of other observed tips). We do not count the ADR event (as it corresponds to an unobserved patient, whose virus is not in our dataset).
Times of DRM Loss To estimate the DRM loss times, we used survival analysis with an exponential (constant hazard) model (Weibull model with β = 1), implemented in Python3 package SurPyval (github. com/derrynknife/SurPyval accessed on 1 April 2023, v0. 10. 10). This model takes as input observations about event durations and estimates the rate at which the event occurs. The input data might be left-, right-, or interval-censored. Left-censored data represent times that are longer than the event occurrences, e. g. , if the DRM loss occurred in exactly 2 years, but the observation was only made after 3 years, the 3-year duration represents a left-censored data point. Right-censored data represent times that are shorter than the event occurrences, e. g. , if for the same DRM loss the only observation was made after 1 year (and observed no DRM loss), the 1-year duration represents a right-censored data point.
Interval-censored data represent cases when both a left-and a right-censored data point are available, e. g. , for the same DRM loss an interval-censored datum might state that it occurred sometime between 1 and 3 years. For each individual represented in our dataset, we extracted at most one data point for the loss survival analysis, as described below. We estimated the loss times for nonpolymorphic DRMs, assuming each DRM can be acquired only in the presence of selective pressure of the ARVs that could provoke it and lost in the absence of such ARVs. However, the treatment composition of individuals was not recorded in our datasets; we only had information on whether they were on treatment or were treatment-naive. Therefore, if a DRM loss did happen (i. e. , a sample with a DRM followed by a subsequent sample without the DRM) and the corresponding individual's state was treatment-experienced, we assumed the treatment was such that it could not provoke this DRM and counted this as a DRM loss, just as in the case of an ARV interruption.
A right-censored data point represents the maximal observed duration during which a mutation loss did not occur. We extracted such points for the individuals who had several consecutive treatment-naive samples with the DRM of interest (and of the subtype of interest) in our metadata: we took the difference in sampling times of the last such sample and the first one. A left-censored data point represents a duration that is longer than the mutation loss time. To estimate such a duration, we needed to know not only (1) the time by which the individual's virus lost the DRM, but also (2) the earliest time by which it could have acquired it (the difference making an upper limit on the loss duration). For (1), we used the time of the earliest sample without the DRM, provided it was preceded by samples with the DRM.
For (2), we used either (2a) the time of the latest sample without the DRM preceding the aforementioned samples (where the DRM was present and then lost), if such sample existed in the metadata, or (2b) if the earliest metadata sample already had the DRM (which implies it corresponded to a resistant tip in the tree), the time of the tip's most recent ancestral node whose status was sensitive (with marginal probability > 0. 95). For individuals for whom both a left-and a right-censored data point was present, we converted them to an interval-censored one. We reported the resulting DRM loss time estimates (i. e. , inverse of the loss rates) for non-polymorphic DRMs with at least 5 left-censored and 5 right-censored data points (interval-censored data points counted as both). We estimated confidence intervals (CIs) as the 2.
If resistance is the reason for a treatment failure, the treatment is changed. Therefore, in the case of treatment failure, there is a window of opportunity for the virus to be transmitted: between the time the viral load increases to transmittable levels and the time when the clinician realizes it has increased and changes treatment. The probability of transmission varies across patients and depends on various factors  [5]. The information collected from the HIV drug resistance tests carried out in the UK since 1996 is available in the UK HIV Drug Resistance Database. The database stores protease (PR) and reverse transcriptase (RT) sequences for about 50% of infected individuals in the UK.
UK HIV Dataset We used the data from the UK HIV Drug Resistance Database containing samples from 1996 to 2016 to estimate transmission mechanisms for different common DRMs. Out of 88,009 initial sequences obtained from the database, the majority were of subtypes B (58,569 sequences, 66. 5%) and C (27,151 sequences, 30. 1%); we also detected 8 D, 1 F, 2 G, and 3 K (<0. 0001%) sequences and 2276 potentially recombinant sequences (2. 6%; in particular 494 A, B, G and 446 B, K recombinants (0. 5%)). We report these and other dataset statistics in Table  3. Table  3. Statistics on the B and C datasets. The "with DRM(s)" statistics count samples with at least one unambiguous resistant amino acid at any DRM position. Samples that contained either non-resistant or ambiguous amino acids at all DRM positions were considered as "without DRMs".
We hence obtained a 40,055-sequence dataset for B and a 19,139-sequence dataset for C. We further filtered these datasets by removing temporal outliers (<2% of sequences), as they could correspond to erroneous dates or poorly sequenced samples (e. g. , due to viral load-mediated stochasticity, base-calling software errors, or issues with sample transport and receipt). The final datasets contained 39,159 sequences for B and 18,809 for C. We detected 161 DRMs found in at least one sequence of the B dataset, and 146 DRMs for C. 31. 4% of B and 27. 4% of C sequences had at least one of these DRMs present, 18. 5% of B and 16. 9% of C sequences had only one mutation, while the others had multiple DRMs present. While the subtypes B and C are different, as are the locations where these subtypes are most prevalent (African countries for C versus the UK and other European countries for B), we did not detect major differences in DRM distribution in the B and C datasets.
Hence, while more C than B sequences correspond to imported cases, the UK health policies must play an important role on their DRM patterns, independently of the subtype. In a recent study Blassel et al. [30]  compared DRMs in a UK and an African datasets and reported that the median number of DRMs in resistant sequences differed between the two datasets (three in the African sequences versus one in the UK sequences). In our case, there was no difference between B and C datasets if all DRMs were considered (median number of one DRM for both B and C datasets in resistant sequences); if we considered only non-polymorphic DRMs, a slight difference appeared (one for B vs. two for C). Detailed statistics on DRM number distributions are shown in Table  A1. There was, however, a significant difference in the TDR distribution: more TDR could be suspected among the B samples (12% of treatment-naive sequences had non-polymorphic DRMs present, while in the C samples, there were only 7% of such sequences).
Drug Resistance Analyses We reconstructed time-scaled phylogenetic trees for B and C datasets and performed ancestral character reconstruction for each of the selected DRMs and positions to look at their transmission patterns. Consistently with what was previously reported in HIV-1 group M studies (of the pol gene  [31]  and of the full-genome  [32] ), we estimated a faster mutation rate (1. 9 × 10 -3 (mutations per site per year)) and a more recent root date  (1965)  for subtype B than for subtype C (1. 4 × 10 -3 ; 1944). More details on B and C datasets can be found in Table  3. On the time-scaled trees we analyzed the transmission patterns of the DRMs found in at least 0. 5% of sequences: 31 DRMs (on 26 different positions) for B and 21 (on 20 different positions) for C. The threshold of 0.
5% permitted us to analyze all the major DRMs while having enough sequences representing these DRMs in the dataset. The major drug resistance patterns found in the B and C datasets are visualized in Figures  3  and 4. The statistics on these DRMs and their loss times are shown in Tables  1 2 3 4. While some of the DRMs (e. g. , RT:M184V) are comparably prevalent in B and C (4. 8% of resistant cases vs. 5. 4%), others are very subtype-specific. For instance, the nonpolymorphic mutation PR:L90M is present in 2. 2% of B resistant cases and only in 0. 6% of C. Another example is the mutations in position RT:106. In the B dataset, the polymorphic DRM RT:V106I is present in 3. 9% of resistant cases and is 30 times more prevalent than the non-polymorphic DRM RT:V106M, which was not selected for our analysis due to its low prevalence; for the C dataset, we have the opposite distribution: RT:V106M is present in 2% of resistant cases and is 16 times more prevalent than RT:V106I, which was not selected for our analyses.
The RT:T215D/S are reversions often developed in patients primarily infected with strains with RT:T215Y/F and hence have a higher fitness  [33]. This is further confirmed by our estimation of the loss times: the loss times of RT:T215D/S are long (9. 3 and 6. 8 years versus 1. 1 and 1. 8 years for RT:T215Y/F, see Table  4 ), which may explain their prevalence in treatment-naive patients. Similarly, we estimated a rather long loss time for RT:K219N (3. 7 years). We did not have enough data to estimate the loss time of RT:V179E. While this mutation is generally considered as non-polymorphic  [34] , its natural presence in treatment-naive patients has been reported for the HIV-1 common recombinant form CRF55_01B  [35]. Using the information from the tree, we refined the mutation statistics further, classifying resistant mutations sources into TDR vs.
ADR and detecting DRM loss events (see Tables  1  and 2 ). The C dataset featured smaller TDR clusters (apart from the polymorphic mutation RT:E138A) than B. This could be explained by multiple introductions of subtype C into different regions of the UK and different risk groups, particularly from Africa via immigration, which is consistent with the higher diversity of the C strains observed in our data (C: 0. 019 vs. B: 0. 014 (mutations per site per branch), Table  3 ) and with the dates of origin of the two UK sub-epidemics (C: 1944 vs. B: 1965, Table  3 ). A large size (e. g. , 78 individuals in the B dataset for RT:K103N) of some of the TDR clusters and a rather high proportion of TDR cases among the resistant ones (see Tables  1  and 2 ) is clinically problematic, as it means a high level of resistant strain transmission, leading to a decrease in treatment choice at the population level.
9% to 2. 2% in the B dataset and from 8. 3% to 11. 6% in the C dataset (Tables  S8  and S31 ). Similarly, the proportion of resistant cases with polymorphic RT:S68G has been increasing from 4. 6% in 2001 to 8. 1% in 2016 in B and from 0. 3% to 0. 9% in C (Tables  S21  and S41 ). The proportion of resistant cases with polymorphic RT:V106I has been increasing in B from 2% in 2001 to 2. 7% in 2016, while the proportion of non-polymorphic RT:V106M (similar to RT:V106I) in C seems to have stabilized at 2% over the last five sampling years (2011-2016, Tables  S23  and S43 ). The proportion of resistant cases with polymorphic RT:V179D has been increasing in B from 1. 3% in 2001 to 2% in 2016, as did the proportion of non-polymorphic RT:V179E (similar to RT:V179D) in C, from 0. 1% in 2006 to 0. 6% in 2016.
DRM Loss Times We estimated the times of DRM loss for non-polymorphic DRMs in our datasets and compared them to the estimates previously reported by  Castro et al. [5]. Castro et al. analyzed 313 patients from the UK Drug Resistance database who were treatment-naive and had a DRM present in their first resistance test (performed between 1997 and 2009), mixing all the subtypes and using survival analysis. We also used survival analysis, but we had a larger dataset that included information not only from the metadata, but also from the tree, and thus analyzed the subtypes separately. Our results and the comparison are shown in Table  4. Overall, out estimates are compatible with those by Castro et al. The visualization using ACR for the DRMs at the position RT:T215 (Figure  5 ) is consistent with the estimated loss patterns.
Discussion We proposed fast maximum-likelihood ACR methods for the investigation of drug resistance patterns in large sequence datasets. Their application to ∼40,000 subtype B and ∼20,000 subtype C sequences from the UK HIV Drug Resistance Database allowed us to investigate the trends in drug resistance patterns between 1996 and 2016 and to estimate the loss times for 25 common non-polymorphic DRMs. An important advantage of our methods is their applicability to very large datasets (dozens of thousands of sequences). Previous studies had to face an uncomfortable choice between using more complex models on filtered data  [9]  or using less accurate (e. g. , parsimony) approaches on full datasets  [4]. Our approach uses a robust maximum likelihood framework and permits the extraction of global drug-resistant patterns from all the available data.
[5]: the CIs of the two studies intersect for all the DRMs except for RT:K103N in the C dataset. The difference for RT:K103N could be explained by the fact that Castro et al. analyzed different subtypes together (though the majority of samples used were from B), while we performed a subtype-specific analysis: our estimate for RT:K103N on the B dataset (2. 0-2. 6 years) is compatible with the one by Castro et al. [5] (2. 0-6. 8 years). Our CIs are systematically narrower than those of Castro et al. While the CIs estimated in our study and the ones by Castro et al. intersect, our point estimates are systematically lower than those of Castro et al. This could be explained by the fact that Castro et al. used only left-censored intervals (i. e. , longer than the event time) in their analyses, while we used both left-(longer) and right-censored (shorter than the event time) intervals (see Section 2 for more details).
Conflicts of Interest: The authors declare no conflict of interest. ( v0. 3. 0b )[42] for basic sequence alignment and tree manipulations, as well as ETE3 framework [43] for basic tree manipulations (format conversion, pruning, etc. ). Survival analysis was performed with Python3 package SurPyval (github. com/derrynknife/SurPyval). Sierra web service [23] for ARV and DRM detection was used via Python3 package sierrapy (github. com/hivdb/sierra-clientaccessed on 1 April 2023).
Quantitative image analysis Granule total number and movement are quantified using the TrackMate plugin for ImageJ software  (Tinevez et al. , 2017) , as shown in figure  3  and in reference  (Juzans et al. , 2020). This open-source plugin allows automated particle tracking. Granules are detected in individual cells based on their estimated diameter (0. 5 µm approximately). Individual granule position is detected at each time point, allowing the automated linking between positions and missing detection filling. We recommend to set the maximum gap-closing frame on 1 to avoid linking detection between two different granules. Generated Excel files provide the number of granules present in the TIRF zone and their LysoTracker mean intensity at each time point. We recommend to consider only strongly fluorescent granules, more likely closer to the plasma membrane.
Image analysis Images captured by the InCell Analyzer 6000 ® were analyzed using Columbus software. Both cell-and parasite-based analyses were performed, and were supplemented by vacuole structure and organization assessments in certain cases. First, a local maxima detection algorithm was applied to the DAPI channel to detect the nuclei, which were then thresholded by size, roundness and DAPI mean intensity to avoid any false detection (as parasite vacuoles could be misidentified as nuclei). Next, cytosols were detected based on the DAPI background signal in this region, after which a spot detection on the DAPI channel was processed to detect parasites in the entire cell layer. Based on this parameter, parasite count per field or per HFF cell monolayer was obtained. In parallel, as parasites were gathered in vacuoles, the vacuole population was split into positive and negative groups.