dataesr/scientific-paragraphs-categorization huggingface dataset:
For sequence analysis and annotation, the BLAST algorithm (www. ncbi. nlm. nih. gov/BLAST), the ISfinder database (www-is. biotoul. fr/), and the ORF finder tool (www. bioinformatics. org/sms/) were utilized. Comparative genome alignments were performed using Mauve, v2. 3. 1  (60). Plasmid typing. The PCR-based replicon typing method for A. baumannii developed by Bertini et al. (22)  was used to determine the plasmid content of the 41 isolates. With this method, the 27 known A. baumannii replicase (rep) genes are grouped into 19 homology groups according to their nucleotide sequence similarities. Six multiplex PCRs were done using primers designed by Bertini et al. (22). Accession number(s). The draft genomes were deposited in the NCBI databases under the accession numbers listed in Table  S3  in the supplemental material.
ACKNOWLEDGMENTS This work was partially financed by the  School of Arts and Sciences Research and Development Council at the Lebanese American University  and by the  National Sustainability Program I (NPU I ; grant number  LO1503 ), provided by the  Ministry of Education Youth and Sports of the Czech Republic , the  Charles University Research Fund-PROGRES  (grant number  Q39 ), and project no. CZ. 02. 1. 01/0. 0/0. 0/16_019/0000787 , Fighting Infectious Diseases. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
COVID-19 is caused by the coronavirus SARS-CoV-2, which jumped into the human population in late 2019 from a currently uncharacterised animal reservoir. Due to this recent association with humans, SARS-CoV-2 may not yet be fully adapted to its human host. This has led to speculations that SARS-CoV-2 may be evolving towards higher transmissibility. The most plausible mutations under putative natural selection are those which have emerged repeatedly and independently (homoplasies). Here, we formally test whether any homoplasies observed in SARS-CoV-2 to date are significantly associated with increased viral transmission. To do so, we develop a phylogenetic index to quantify the relative number of descendants in sister clades with and without a specific allele. We apply this index to a curated set of recurrent mutations identified within a dataset of 46,723 SARS-CoV-2 genomes isolated from patients worldwide.
This is being supplemented by increasing raw sequencing data available through the European Bioinformatics Institute and NCBI Short Read Archive, together with data released by specific genome consortiums, including COVID-19 Genomics UK (https://www. cogconsortium. uk/data/). Research groups around the world are continuously monitoring the genomic diversity of SARS-CoV-2, with a focus on the distribution and characterisation of emerging mutations. Mutations within coronaviruses, and indeed all RNA viruses, can arrive as a result of three processes. First, mutations arise intrinsically as copying errors during viral replication, a process which may be reduced in SARS-CoV-2 relative to other RNA viruses, due to the fact that coronavirus polymerases include a proof-reading mechanism  8, 9. Second, genomic variability might arise as the result of recombination between two viral lineages coinfecting the same host  10.
gisaid. org) and provided in Supplementary Data 1. The map in Fig. 1b  was created using the R package rworldmap using the public domain Natural Earth data set. the entire set of 185 recurrent mutations seem to fit the expectation for neutral evolution with respect to transmissibility, with a mean and median overall log 10 RoHO score of -0. 001 and -0. 02, respectively. Moreover, the distribution of individual site-specific RoHO scores is symmetrically distributed around 0 with 97/185 mean positive values and 88/185 negative ones. To summarise, we would expect that recurrent mutations should be the best candidates for putative adaptation of SARS-CoV-2 to its novel human host. However, none of the recurrent mutations in circulation to date shows evidence of being associated with viral transmissibility.
Discussion In this work, we analysed a data set of >46,700 SARS-CoV-2 assemblies sampled across 99 different countries and all major continental regions. Current patterns of genomic diversity highlight   multiple introductions in all continents (Fig. 1  and Supplementary Figs. 1 2 3 ) since the host-switch to humans in late 2019  [1] [2] [3] [4]. Across our data set, we identified a total of 12,706 mutations, heavily enriched in C→U transitions, of which we identified 398 strongly supported recurrent mutations (Supplementary Data 3 and Supplementary Figs. 4  and 5 ). Employing a phylogenetic index (RoHO) to test whether these recurrent mutations contribute to a change in transmission, we found no candidate convincingly associated with a significant increase or decrease in transmissibility (Figs.
Methods Data acquisition. We downloaded 48,454 SARS-CoV-2 assemblies were downloaded from GISAID on 30/07/2020 selecting only those marked as 'complete', 'low coverage exclude' and 'high coverage only'. To this data set, all assemblies of total genome length <29,700 bp were removed, as were any with a fraction of 'N' nucleotides >5%. In addition, all animal isolate strains were removed, including those from bat, pangolin, mink, cat and tiger. All samples flagged by NextStrain as 'exclude' (https://github. com/nextstrain/ncov/blob/master/defaults/exclude. txt) as of 30/07/2020 were also removed. Twenty one further accessions were also filtered from our phylogenetic analyses as they appeared as major outliers following phylogenetic inference and application of TreeShrink 51 despite passing other filtering checks.
This left 46,723 assemblies for downstream analysis. A full metadata table, list of acknowledgements and exclusions is provided in Supplementary Data 1. Multiple sequence alignment and maximum likelihood tree. All 46,723 assemblies were aligned against the Wuhan-Hu-1 reference genome (GenBank NC_045512. 2, GISAID EPI_ISL_402125) using MAFFT v7. 471  52  implemented in the rapid phylodynamic alignment pipeline provided by Augur 6. 3. 0 (github. com/ nextstrain/augur). This resulted in a 29,903 nucleotide alignment. As certain sites in the alignment have been flagged as putative sequencing errors (http://virological. org/t/issues-with-sars-cov-2-sequencing-data/473), we followed two separate masking strategies. The first masking strategy is designed to test the impact of the inclusion of putative sequencing errors in phylogenetic inference, masking several sites within the genome (n = 68) together with the first 55 and last 100 sites of the alignment (the list of sites flagged as 'mask' is available at https://github.
com/W-L/ ProblematicSites_SARS-CoV2/blob/master/problematic_sites_sarsCov2. vcf, accessed 30 Jul 2020)  30. We also employed a less stringent approach, following the masking strategy employed by NextStrain, which masks only positions 18,529, 29,849, 29,851 and 29,853 as well as the first 130 and last 50 sites of the alignment. A complete list of masked positions is provided in Supplementary Data 5. This resulted in two masked alignments of 46,723 and 46,745 assemblies with 12,706 and 12,807 SNPs, respectively. Subsequently, for both alignments, a maximum likelihood phylogenetic tree was built using IQ-TREE 2. 1. 0 COVID release (https://github. com/iqtree/iqtree2/ releases/tag/v2. 1. 0) as the tree-building method  53. The resulting phylogenies were viewed and annotated using ggtree v1. 16.
6  54  (Fig. 1  and Supplementary Fig. 1 ). Site numbering and genome structure are provided for available annotations (nonoverlapping open reading frames (ORFs)) using Wuhan-Hu-1 (NC_045512. 2) as reference. Recombination analysis. In order to test for the presence of recombination, we performed a LD analysis considering whether the correlation between SNPs tends to disappear with an increase in the distance separating them on the genome. A classical approach to do so is the use of LD decay curves, which represent LD as a function of the distance separating each SNP pair. We calculated metrics of LD (r 2 and D') across all pair-wise combinations of variant sites using Tomahawk 0. 7. 0 (https://github. com/mklarqvist/tomahawk). The relationship between LD and distance yielded a regression coefficient of 3.
2 ). Homoplasy screen. The resulting maximum likelihood trees were used, together with the input alignments, to rapidly identify recurrent mutations (homoplasies) using HomoplasyFinder v0. 0. 0. 9  1, 56. HomoplasyFinder employs the method first described by Fitch  57  , providing, for each site, the site-specific consistency index and the minimum number of changes invoked on the phylogenetic tree. All ambiguous sites in the alignment were set to 'N'. HomoplasyFinder identified a total of 5710 homoplasies, which were distributed over the SARS-CoV-2 genome (Supplementary Fig. 4 ). For the less stringent masking of the alignment, HomoplasyFinder identified a total of 5793 homoplasies (Supplementary Fig. 5 ). As previously described, we filtered both sets of identified homoplasies using a set of thresholds attempting to circumvent potential assembly/sequencing errors (filtering scripts are available at https://github.
We defined homopolymer regions as positions on the Wuhan-Hu-1 reference with at least four repeated nucleotides. While homopolymer regions can arise through meaningful biological mechanisms, for example polymerase slippage, such regions have also been implicated in increased error rates for both nanopore  58  and Illumina sequencing  59. As such, homoplasies detected near these regions (±1 nt) could have arisen due to sequencing error rather than solely as a result of underlying biological mechanisms. If this were true, we would expect the proportion of homoplasic sites near these regions to be greater than that of homopolymeric positions across the entire genome. We tested this by identifying homopolymer regions using a custom python script (https://github. com/cednotsed/genome_homopolymer_counter) and performing a binomial test on the said proportions.
Visualisation of PCA projections (Supplementary Fig. 10a ) suggested that there was no hidden structure introduced by our homoplasy filtering steps. The first two PCs accounted for 56% of the variance and were mostly loaded by the variables encoding the reference and variant nucleotides (Supplementary Fig. 10b ). Annotation and characterisation of homoplasic sites. All variable sites across the coding regions of the genomes were identified as synonymous or nonsynonymous. This was done by retrieving the amino acid changes corresponding to all SNPs at these positions using a custom Biopython (v. 1. 76) script (https://github. com/cednotsed/nucleotide_to_AA_parser. git). The ORF coordinates used (including the ORF1ab ribosomal frameshift site) were obtained from the associated metadata according to Wuhan-Hu-1 (NC_045512.
2 ). We define this as the RoHO index (full associated code available at https:// github. com/DamienFr/RoHO)  38. HomoplasyFinder 56 flags all nodes of a phylogeny corresponding to an ancestor that acquired a homoplasy. We only considered nodes with at least two descending tips carrying either allele and with no children node embedded carrying a subsequent mutation at the same site (see Fig. 2 ). For each such node in the tree, we counted the number of isolates of each allele and computed the RoHO index. We finally restricted our analysis to homoplasies having at least n = 3 individual RoHO indices (i. e. for which three independent lineages acquired the mutation). The latter allows us to consider only nodes for which we have multiple supported observations within the phylogeny, conveniently accounting for phylogenetic uncertainty.
Paired t tests were computed for each homoplasy to test whether RoHO indices were significantly different from zero. To validate the methodology, this analysis was carried out on data analysed using two different masking strategies (Fig. 3 , Supplementary Fig. 13  and masked sites available in Supplementary Data 5). Full metadata associated with each tested site, including the number and associated countries of descendant offsprings, are provided in Supplementary Data 4. Assessing RoHO performance. To assess the performance of our RoHO index, we performed a set of simulations designed to test the distribution of RoHO values under a neutral model. We simulated a 10,000 nucleotide alignment comprising 1000 accessions using the rtree() simulator available in Ape v5. 3  60  and genSeq from the R package PhyTools v0.
2 ). We further enforced a minimum number of three replicates (Fig. 3 , Supplementary Fig. 13  and Supplementary Data 4). While we discarded homoplasies located on 'embedded nodes' to avoid pseudoreplication (see node 1 in Fig. 2 ), we note that including such sites has no impact on our results (Supplementary Fig. 11b ). In addition, we assessed the statistical power to detect significant deviations from neutrality of the RoHO index according to (i) the number of independent emergences of a homoplasy in the phylogeny and (ii) the imbalance between offspring number for each allele (i. e. fitness differential conferred by the carriage of the derived allele). To do so, we generated 1000 replicates for each combination of independent emergences (counts) of a homoplasy and corresponding fitness differential values using results from both masked alignments.
Acknowledgements L. v. D. and F. B. acknowledge financial support from the  Newton Fund UK-China NSFC initiative  (grant  MR/P007597/1 ) and the  BBSRC  (equipment grant  BB/R01356X/1 ). L. v. D. is supported by a  UCL Excellence Fellowship. D. R. is supported by an  NIHR Precision AMR  award. L. P. S. acknowledges funding from the  Antimicrobial Resistance Cross-council Initiative  supported by the seven UK research councils (grant  NE/ N019989/1 ). Computational analyses were performed on UCL Computer Science cluster and the  South Green  bioinformatics platform hosted on the CIRAD HPC cluster. We additionally wish to acknowledge the very large number of scientists in originating and submitting laboratories who have readily made available SARS-CoV-2 assemblies to the research community. We additionally wish to thank helpful comments on Twitter and bioRxiv, in particular from  Harald Ringbauer ,  Palle Villesen ,  Sally Otto  and  Daniel Falush.
Code availability All codes used to generate the RoHO method presented in this manuscript are available at https://github. com/DamienFr/RoHO  38. In addition, we provide links to the code used to conduct homoplasy filtering (https://github. com/liampshaw/CoV-homoplasyfiltering), per site annotations (https://github. com/cednotsed/nucleotide_to_AA_parser. git) and to assess the number of homopolymer regions in the SARS-CoV-2 genome (https://github. com/cednotsed/genome_homopolymer_counter).
The acceleration of DNA sequencing in samples from patients and population studies has resulted in extensive catalogues of human genetic variation, but the interpretation of rare genetic variants remains problematic. A notable example of this challenge is the existence of disruptive variants in dosage-sensitive disease genes, even in apparently healthy individuals. Here, by manual curation of putative loss-of-function (pLoF) variants in haploinsufficient disease genes in the Genome Aggregation Database (gnomAD) 1 , we show that one explanation for this paradox involves alternative splicing of mRNA, which allows exons of a gene to be expressed at varying levels across different cell types. Currently, no existing annotation tool systematically incorporates information about exon expression into the interpretation of variants.
Our annotation is fast, flexible and generalizable, making it possible for any variant file to be annotated with any isoform expression dataset, and will be valuable for the genetic diagnosis of rare diseases, the analysis of rare variant burden in complex disorders, and the curation and prioritization of variants in recall-by-genotype studies. A primary challenge in the use of genome and exome sequencing to predict human phenotypes is that our capacity to identify genetic variation exceeds our ability to interpret their functional impact  3, 4. One underappreciated source of variability for variant interpretation involves differences in alternative mRNA splicing, which enables exons to be expressed at different levels across tissues. These expression differences mean that variants in different regions of a gene can have different phenotypic outcomes depending on the isoforms they affect.
Isoform diversity and variant interpretation We find that isoform diversity is a contributor to the paradoxical finding of disruptive variants in dosage-sensitive disease genes in ostensibly healthy individuals. In the gnomAD database, we identify 401 high-quality pLoF variants that pass both sequencing and annotation quality filters in 61 haploinsufficient disease genes in which heterozygous pLoF variants are established to cause severe developmental delay phenotypes with high penetrance (Methods). Given the severity of these phenotypes and their extremely low prevalence worldwide, ranging from 1 in 10,000 to less than 1 in a million, very few, if any true pLoF variants would be expected to be found in the gnomAD population. As such, most or all of these observed pLoF variants are likely to be sequencing or annotation errors  12.
pext score summarizes isoform expression The advent of large-scale transcriptome sequencing datasets, such as GTEx  2  , provides an opportunity to incorporate cross-tissue exon expression into variant interpretation. However, the current formats of these databases do not readily allow for unbiased estimation of exon expression. The GTEx web browser offers information on exon-level read pileup across tissues, but this approach is confounded by technical artefacts such as 3′ bias  13  (preferential coverage of bases close to the 3′ end of a transcript) (Supplementary Fig. 2a ). Such systematic biases mean that simple exon-level coverage in a transcriptome dataset cannot be used as a reliable proxy for exon expression, especially in longer genes (Fig. 2a , Supplementary Fig. 2b ). Isoform quantification tools provide estimates of isoform expression levels that correct, albeit imperfectly  13, 14  , for confounding by 3′ bias as well as other technical artefacts such as isoform length, isoform GC content, and transcript sequence complexity  [15] [16] [17].
Here, we use isoform-level quantifications from 11,706 tissue samples from the GTEx v7 dataset to derive an annotation-specific expression metric. For each tissue, we annotate each variant with the expression of every possible consequence across all transcripts, which can be used to summarize expression in any combination of tissues of interest. We first compute the median expression of a transcript across tissue samples, and define the expression of a given variant as the sum of the expression of all transcripts for which the variant has the same annotation (Fig. 2a , Supplementary Fig. 3a ). By normalizing the expression of the annotation to the total gene expression, we define a metric (proportion expression across transcripts, or 'pext'), which can be interpreted as a measure of the proportion of the total transcriptional output from a gene that would be affected by the variant annotation in question (Supplementary Fig.
3b ). The pext metric allows for quick visualization of the expression of exons across a gene. In Fig. 2b , transcript-expression based annotation is shown for TCF4, a haploinsufficient gene in which heterozygous variants result in Pitt-Hopkins syndrome  18  , a highly penetrant disorder associated with severe developmental delay. This gene contains 20 unique high-quality pLoF mutations across 56 individuals in the gnomAD database. All 20 variants lie on exons with no evidence of expression across the GTEx dataset (Fig. 2b , Supplementary Fig. 4 ), which indicates that functional TCF4 protein can be made in the presence of these variants. This visualization is now available for all genes in the gnomAD browser (https://gnomad. broadinstitute. org), and can aid in the rapid identification of variants occurring on exons with little to no evidence of expression in GTEx.
Functional validation of pext To explore whether expression-based annotation marks functionally important regions, we compared the distribution of the pext metric in evolutionarily conserved and unconserved regions using phyloCSF  19. Exons with patterns of multi-species conservation consistent with coding regions have higher phyloCSF scores, and should exhibit detectable expression patterns, whereas regions with lower scores will be enriched for incorrect exon annotations, which are expected to have little evidence of expression in a population transcriptome dataset. As expected, we observe significantly lower expression for unconserved regions, and near-constitutive expression in highly conserved regions (Fig. 3a , Supplementary Fig. 5a ). This difference remains statistically significant after correcting for exon length (logistic regression P < 1.
Nonetheless, we find support for some highly conserved CDS regions, several of which show evidence of transcription in fetal tissues, underlining the importance of incorporating several isoform expression datasets for interpretation (Supplementary Fig. 6d ). Non-synonymous variants found on constitutively expressed regions would be expected to be more deleterious than those on regions with no evidence of expression. To test this, we defined expression bins based on the average pext value across GTEx tissues, in which an average pext value less than 0. 1 was defined as low (or unexpressed), above 0. 9 as high (or near-constitutive) and intermediate values as medium expression. We compared the mutability-adjusted proportion singleton (MAPS), a measure of negative selection on variant classes  21  , partitioned on the loss-of-function observed/expected upper-bound fraction (LOEUF) decile, a measure of constraint against pLoF variants in the gnomAD dataset 1 in each of these expression bins.
Number of reads Transcript-level expression Fig. 2  | Summary of transcript-expression based annotation method. a, Overview of transcript-aware annotation. Most genes have many annotated isoforms, which can have varying expression patterns across tissues. Using the number of reads aligning to exonic regions in transcriptome datasets as a proxy for exon expression (top, black) has confounding effects, due to 3′ bias. In this example, although exons 3 and 8 have markedly different expression levels in brain cortex, the number of reads aligning to the two exons is similar, and this masks the differences in exon usage. Transcript-aware annotation defines the expression of every variant as the sum of transcripts that have the same annotation. The resulting transcript-level expression plots do not exhibit 3′ bias, and reveal differences in exon usage, such as those in exons 3 and 8, across tissues.
b, Example of utility of transcript-expression based annotation. There are 20 high-quality pLoF variants in the haploinsufficient developmental delay gene TCF4 in gnomAD, annotated as dashed lines and arrows. All 20 variants have no evidence of expression in the GTEx dataset, which suggests that functional TCF4 protein can be made in the presence of these variants. Supplementary Fig. 5c , Supplementary Table  5a , b). This information is complementary to existing variant prioritization tools such as PolyPhen-2  22  (Supplementary Fig. 5d , Supplementary Table  5c ). This skew of non-synonymous variation in high-expressed regions suggests that variation arising in such exons tends be more deleterious, whereas non-synonymous variants on regions with low expression are similar to missense variants in their inferred deleteriousness.