dataesr/scientific-paragraphs-categorization huggingface dataset:
Phylogeny based on AlbFlavi2 and AlbFlavi36 DNA sequences from AlbFlavi2 and AlbFlavi36 were aligned with MUSCLE algorithm  [48]. Exogenous virus sequences were used as outgroup to determine the direction of character transformations. Phylogenetic trees were obtained by parsimony analysis implemented in the PAUP* software package (version 4. 0), by using gap as 5th state or not (data not shown) and nearest neighbour interchange (NNI) for tree rearrangements  [49]. Pilot experiment assessing the association between AlbFlavi2/Albflavi36 and vector competence of Ae. albopictus at the individual level Logistic linear regression models were used to test association between the presence of AlbFlavi2 and Alb-Flavi36 and viral dissemination in single mosquitoes from the Foshan colony and a field-collected African population (Tibati collected in Tibati, Cameroon in 2018).
Foshan and Tibati mosquitoes were experimentally infected with DENV (DENV-2, Accession number: MK268692  [50] ) and CHIKV (06. 21, Accession number: AM258992  [51] ) provided in blood meals as described in Amraoui et al. 2019  [52]. All surviving mosquitoes were examined individually at 14 days post-infection to define (i) the infectious status; presence of infectious particles in heads was estimated by focus fluorescent assay on C6/36 Ae. albopictus cells  [52] , and (ii) the presence of AlbFlavi2 and/or Alb-Flavi36 by PCR on DNA extracted from bodies and thorax. Statistical analyses were conducted using the Stata software (StataCorp LP, Texas, and USA). Pvalues <0. 05 were considered significant.
NIRVS analysis We further chose to study NIRVS distribution in 13 Ae. albopictus populations selected previously for assessing genetic diversity based on epidemiological data of arboviral outbreaks and Ae. albopictus widespread out if its native range  [17, 34, 35] ; we analyzed 10 females and 10 males per population (Table  1 ). NIRVS distribution was not homogenous across populations (Figure  2 , Supplementary Figure  2 ). AlbFlavi1 (Figure  2 (A)) and CSA-JJL (Figure  2 (H)) were the most widespread NIRVS, being detected in 85-100% of tested mosquitoes, respectively. On the opposite, AlbFlavi10 was the rarest NIRVS, being found in 17% of the tested mosquitoes. Despite displaying the highest presence of AlbFlavi2 and AlbFlavi4 (present in most populations), Brazilian populations showed the lowest number of NIRVS, with the absence of Alb-Flavi36 and AlbFlavi10 and the lower frequency of Alb-Flavi41 compared to the other Ae.
Sequence polymorphism and evolution of AlbFlavi2 and AlbFlavi36 AlbFlavi2 and AlbFlavi36 of individual mosquitoes were successfully amplified in 10 out of the 13 Ae. albopictus populations studied (Figure  2 ). Sequences of AlbFlavi2 and AlbFlavi36 amplified fragments were aligned with the closely related exogenous viral sequences: the Kamiti River virus (KRV) for AlbFlavi2 and the Aedes Flavivirus (AeFV) for AlFlavi36. The p-distance calculated between AlbFlavi2 sequences revealed a moderate diversity, with values ranging from 0% to 3. 9% (data not shown). All Alb-Flavi2 sequences underwent deletion events that were sometimes shared by several individuals from different populations or observed in a specific Ae. albopictus population (Figure  4 ; Supplementary AlbFlavi36). The p-distance values calculated between Alb-Flavi36 sequences were low and varied only from 0% to 0.
9% (data not shown). Therefore, AlbFlavi36 sequences mapping to an intergenic region  [24]  appeared more similar between individuals than Alb-Flavi2 sequences (which are integrated in piRNA cluster 27  [53] ). 2 ). Populations were characterized as high or low frequencies for each NIRVS by using the median and were analyzed with DE data by logistic regression models to find any association. Odds ratio (OR) >1 and <1 indicated positive and negative association respectively. 95% CI: 95% confidence intervals; N: sample size. We further performed phylogenetic analysis to describe the evolutionary history of AlbFlavi2 and Alb-Flavi36 (Figure  4 ; Supplementary Figure  4 ). As expected, AlbFlavi2 displayed higher divergence than AlbFlavi36. The resulted trees based on AlbFlavi2 supported four major clusters sharing sequences with the same deletion events (Figure  4 ).
Funding This study was supported by the  Institut Pasteur , the  French Government 's  Investissement d'Avenir program ,  Laboratoire d'Excellence "Integrative Biology of Emerging Infectious Diseases " (grant number  ANR-10-LABX-62-IBEID  to A-BF), the  European Union 's  Horizon 2020 research and innovation program ,  ZIKAlliance program  under ZIKAlliance grant agreement no. 734548  to  A-BF , and the  European Research Council, NIRV_HOST_INT Council  (grant agreement number  682394 -NIRV_HOST_INT  to VH and MB).
2E, H, K, N ) and newborn pig trachea (Fig. 2F, I, L, O ). Hybrid cells were also detected by qRT-PCR in a fully independent HAEC culture, after isolation of the cells using C1 technology (Fluidigm) and quantification of gene expression with a Biomark (Fluidigm). Cells isolated with the C1 were visually inspected, and these experimental settings ensured the absence of cell doublets. Four cells out of 74 expressed GC-specific genes (namely MUC5AC, MUC5B and TFF3), together with MCC-specific genes (FOXJ1), and more specifically, immature MCC genes (PLK4, MYB and CDC20B)  (Revinski et al. , 2018)  (Fig. S7A, B ). This result was confirmed after re-analyzing a recently published dataset  (Plasschaert et al. , 2018)  (Fig. S7C, D ). A further confirmation came from the detection at the protein level of cells that were simultaneously labeled for MUC5AC and acetylated tubulin, a specific protein marker of the cilia (Fig.
2Q ). We then explored the same dataset with Palantir, another algorithm that models cell trajectory, with which we confirmed the presence of GCs on the MCC branch (Fig. S7E ). The score for differentiation potential was highest for cycling basal cells. A high score was also found in the MCC branch in a region containing both CCs and GCs, before the gap separating them from MCCs (Fig. S7F ), further suggesting a high probability to differentiate into at least two distinct trajectories. Estimation of gene expression trends showed an upregulation and then a downregulation of both MUC5AC and MUC5B along the pseudotime in cells committed to the MCC lineage (Fig. S7G ). Finally, computing branch probabilities of randomly selected GCs on the MCC branch showed that some of them have between 24. 7% and 49.
, 2018). Our repertoire of KRT expression during airway regeneration was based on pseudotime ordering in our Pneumacult ALI 28 dataset. Our analysis confirmed the presence of KRT5 and KRT14 in BCs, of KRT4 and KRT13 in supraBCs, and the expression of KRT8 in luminal cell types (CCs, GCs and MCCs) (Fig. 4A, E ). Unlike recent data obtained by Plasschaert et al. under similar conditions  (Plasschaert et al. , 2018) , who showed parallel RNA expression of KRT13 and KRT4, we consistently noticed that expression profiles of KRT13 and KRT4 were slightly de-correlated, with KRT13 detected at earlier pseudotimes than KRT4. This was confirmed at the protein level by a quantification of immunostainings of the proportion of KRT5+/KRT13+ and KRT5+/KRT4+ double-positive cells (Fig. 4B ). Fig. 4C  shows that there were more KRT5+/KRT13+ (7.
Conclusions We provide several novel insights in the dynamics of airway differentiation by positioning goblet cells as possible precursors of multiciliated cells: this illustrates how cells carrying specialized function, i. e. club and goblet cells, can still constitute differentiation intermediates for other specialized cells, i. e. multiciliated cells. We also identify subpopulations of basal, suprabasal, club and multiciliated cells. Our dataset also provides extensive characterization of the deuterosomal cell population, an intermediate state before the formation of multiciliated cells. After establishing a comprehensive repertoire of keratin expression, we show that monitoring 'keratin switch' during differentiation could be self-sufficient to establish the different cell identities. Our improved characterization of the different signaling pathway components detects putative Notch repressors that probably contribute to Notch signal shutdown at the deuterosomal stage, and details Wnt pathway activity within the basal cell compartment.
Single-cell RNA-seq We followed the manufacturer's protocol (Chromium Single Cell 3′ Reagent Kit, v2 Chemistry) to obtain single cell 3′ libraries for Illumina sequencing. Libraries were sequenced with a NextSeq 500/550 High Output v2 kit (75 cycles) that allows up to 91 cycles of paired-end sequencing: read 1 had a length of 26 bases that included the cell barcode and the UMI; read 2 had a length of 57 bases that contained the cDNA insert; index reads for sample index of eight bases. Cell Ranger Single-Cell Software Suite v1. 3 was used to perform sample demultiplexing, barcode processing and single-cell 3′ gene counting using standard default parameters and human build hg19, pig build sus scrofa 11. 1 and mouse build mm10. All single-cell datasets that we generated, and the corresponding quality metrics are displayed in Table  S6  and were deposited on the Gene Expression Omnibus portal under the series number GSE121600.
Clustering robustness In order to best determine the key steps in the differentiation process, a customized method was implemented to analyze clustering robustness to dataset perturbation. For all possible numbers of clusters (from 2 to 9), multiple subsets of the studied datasets were created (10 subsets with 10% of the cells randomly removed each time) and clustering was performed multiple times on each subset with changing settings of the seed parameter. The result of those clusterings were stored in a (n cells)² stability matrix, containing for each pair of cells 1 or 0 depending on whether the cells are clustered together (1) or not (0). This stability matrix was then transformed in a Euclidean distance matrix between cells and then divided into the used k number of clusters k using hierarchical clustering (hclust with 'average' method).
To identify the optimal number of clusters, a visual inspection of the elbow plot of the average intra-stability (mean stability within each cluster) and the average inter-stability (mean stability between each cluster) was carried out. Cells with a stability metric less than 70% were labeled as 'unassigned', owing to the high clustering variability between each round of clustering, then removed from further analysis of the time course data. Cell clustering was performed using SIMLR ( package version 1. 4. 1) (Wang et al. , 2017). Heatmaps for the clustering of each dataset are shown in Table  S8.
Differential analysis To further analyze the robustness of each step of the differentiation process, we tested the robustness of the cell type marker gene identification through differential gene expression analysis. Differential expression analysis was performed using edgeR (package version 3. 22)  (Robinson et al. , 2010). In a one versus all differential analysis, a pool of 100 cells from one cluster were analyzed against an equal mixture of cells from all other clusters. In a one versus one differential analysis, pools of cells of the same size were compared. Those differential analysis were performed multiple times (10 times) on different pool of cells and the DEG identified were compared between each pool of cells using the rank-rank hypergeometric overlap algorithm  (Plaisier et al.
Trajectory inference Trajectory inference was performed using monocle 2 ( package version 2. 8)  (Qiu et al. , 2017). Cell ordering was based on highly variable genes (∼200-500 genes) selected by their expression dispersion. Monocle analysis on the aggregated time points was carried out on raw counts after library size correction (downsampling). Branch building was performed using BEAM analysis from Monocle, and corresponding differential analysis was carried out after a cross comparison of a group of cells along the pseudotime (before branching, after branching and at the branch end) using Seurat 1 versus 1 differential analysis.
Data visualization All graphs were generated using R (ggplot2). Heatmaps were obtained using pheatmap (no clustering used, genes ordered by their expression in pseudotime or in cluster, cells ordered by pseudotime or cluster). Heatmaps show smoothed gene expression values: for each gene, normalized gene expression values were first transformed into z-scores, then averaged across 10 neighboring cells in the chosen ordering ( pseudotime only or pseudotime in clusters). Single gene representation: for the sake of clarity, only cells with expression levels above the top 50 percentiles for that gene are represented.
Individual sample analysis Each sample of our study was reanalyzed with less stringent parameters to identify rare or transitory cell types or gene expression events Preprocessing, normalization and clustering Individual dataset analysis was performed using Seurat standard analysis pipeline. Briefly, cells were first filtered based on number of expressed features, dropout percentage, library size and mitochondrial gene percentage. Thresholds were selected by visually inspecting violin plots in order to remove the most extreme outliers. Genes expressing fewer than five UMI across all cells were removed from further analysis. Cell-level normalization was performed using the median UMI counts as a scaling factor. Highly variable genes were selected for following analyses based on their expression level and variance.
PCA analysis was performed on those genes, the number of PCs to use was chosen upon visual inspection of the PC variance elbowplot (∼10 to 20 PCs depending on the dataset). Clustering was first performed with default parameters and then by increasing the resolution parameter above 0. 5 to identify small clusters (but with the knowledgeable risk of splitting big clusters due to high gene expression variability). Differential analysis was again performed using Seurat FindAllMarkers and FindMarkers functions based on non-parametric Wilcoxon rank sum test. Gene Set Enrichment analysis was performed using fgsea R package with the following gene sets reactome. db (R package) and GO cellular component (Broad Institute GSEA MSigDB) genesets. Molecular function enrichment analysis was performed using Ingenuity Pathway Analysis (Qiagen).
Velocity RNA velocity was calculate using latest release of velocyto pipeline (velocyto. org/) using standard parameters: GTF file used for Cell Ranger analysis and the possorted_genome_bam. bam, Cell Ranger output alignment file. From the loom file that contains a count table of spliced and unspliced transcripts, the gene. relative. velocity. estimates function was used on cell type marker genes. The resulting expression pattern of unspliced-spliced phase portraits shows the induction or repression of those marker genes from one cell type to the next. We used velocyto package version 0. 5  (La Manno et al. , 2018).
Trajectory inference using Palantir algorithm Palantir analysis was used as an integrated function of the Scanpy workflow  (Wolf et al. , 2018). The filtered raw count matrix was loaded into Scanpy, along with the cell type annotation (Scanpy v1. 4, Python 3. 7); each cell was normalized to the total count over all genes (without log transform) before running Palantir  (Setty et al. , 2019). The first 14 principal components were used to compute the diffusion map. The corresponding t-SNE embedding was obtained using the first two diffusion components. A start cell was randomly selected among the cycling basal cell cluster to infer trajectories and the associated terminal states. In the process, each cell of the dataset was associated with a probability to differentiate into each of the terminal states identified.
Funding This project was funded by grants to P. B. from the  Fondation pour la Recherche Medicale  ( DEQ20180339158 ), the  Labex Signalife  ( ANR-11-LABX-0028-01 ), the  Association Vaincre la Mucoviscidose  ( RF20180502280 ) and the  Chan Zuckerberg Initiative  ( Silicon Valley Community Foundation ,  2017-175159-5022 ). The UCAGenomiX platform, a partner of the  National Infrastructure France Geńomique , is supported by the  Commissariat aux Grands Investissements  ( ANR-10-INBS-09-03  and  ANR-10-INBS-09-02 ) and  Canceropôle PACA. Deposited in PMC  for immediate release.
Acknowledgments Authors' roles: RC has designed and collected the data of the QUALYOR cohort, participated in data analysis, drafted the manuscript, analysis design. ES-R has designed and collected the data of the OFELY cohort, participated in data analysis, reviewed and corrected manuscript. RZ has developed the strax software, reviewed and corrected manuscript. PDD has designed the OFELY study. DL has conducted the cost effectiveness analysis. MB has conducted the statistical analysis. EL has designed and collected the data of the QUALYOR cohort, reviewed and corrected manuscript. ES has designed the analysis drafted manuscript, participated in data analysis.
In fact, precise gene identifications and detailed virulence portraits can be determined in a single session of sequencing, while bioinformatics analysis for foodborne pathogens provides valuable resources for databases and collections of genomic features  [22]. One of the most common foodborne related diseases is the staphylococcal food poisoning (SFP), resulting from the consumption of foods containing sufficient amounts of staphylococcal enterotoxins (SEs)  [23]. SFP is associated with enterotoxigenic S. aureus (S. aureus) strains that may produce one or more heat-stable enterotoxins in food, causing vomiting and abdominal pain in humans  [24]. In 2015, 9. 9% of all European outbreaks were caused by SEs  [25] , indicating the importance of surveillance plans and publication of updated reports for risk assessment evaluations.
The genes encoding for S. aureus VFS are carried by mobile genetic elements (MGEs) and consist of prophages, plasmids, transposons, and S. aureus pathogenicity islands (SaPIs)  [42, 43] , which encode, among the others, enterotoxins and adhesins  [23, 44]. MGEs mediate their own transfer and integration into new genomic sites and, with a phenomenon called horizontal gene transfer (HGT), also among other bacteria, cause adaptive consequences, such as the transfer and acquisition of antibiotic resistance genes  [45]. In this study, six enterotoxigenic S. aureus strains were genome-sequenced. A report of their virulence portraits is presented, and the databases used for the detection of the genomic features are described. The strains were also used for the production of naturally contaminated cheeses, suggesting new opportunities for the production of reference materials for inter-laboratory tests, required by regulations  [46].
Identification of Virulence Factors and Genomic Analysis The genomes were interrogated for a pool of 1300 genes, including VFs reported for staphylococci and determined using a combination of a database built for this study (Supplementary Material Database S1) and the PATRIC tool (http://patricbrc. org) (v3. 5. 41)  [59]. The VFs were classified into six functional categories: genes involved in adherence, exoenzymes, host immune evasion factors, iron uptake and metabolism, toxins and transcriptional regulatory elements. Antibiotic resistance genes were detected from the assemblies using ABRicate v0. 7 (https://github. com/tseemann/abricate) interrogating the comprehensive antibiotic resistance database (CARD) (https://card. mcmaster. ca). ABRicate was used for the creation of three novel databases built for the detection of SaPIs (Supplementary Material Database S5), plasmids (S7), and genes coding for enterotoxins (S3).
In order to avoid ambiguities in taxonomic inferences, the average nucleotide identity (ANI) comparison of the six genomes was carried out using two independent ANI determination methods: ANI Calculator (http://enve-omics. ce. gatech. edu/ani/index) and EzGenome (http://www. ezbiocloud. net/ezgenome/ani). The genomes of 29 enterotoxigenic S. aureus strains were used for the genomic comparison of the six strains presented in this study and included in the constructing of the pan-genome, which was annotated with the use of the Clusters of Orthologous Groups (COGs) database  [66]. The genomes were obtained from the assembly database (https://www. ncbi. nlm. nih. gov/assembly/), and the details of the strains are reported in the Supplemental Table  S2. All the tools and software were used with default settings if not specified, and data in figures or intext are presented as mean +/one standard error of the mean unless otherwise indicated.
Enterotoxin Insight Further analyses were performed for the gene coding staphylococcal enterotoxins. A novel database was built for ABRicate starting from amino acid sequences extracted from GenBank nucleotide and uniprot databases  [72]. The sequences manually annotated (Swiss-prot) and unreviewed enterotoxin variants (TrEMBL) were translated with the EMBOSS tool [73] and combined to build the final database (Supplementary Material Database S3), which included, in total, 99 sequences of staphylococcal enterotoxins and their variants. The strain ProNaCC1 was positive for the genes of the enterotoxins SED, SElJ, SER, SelX, and SElW; the strains ProNaCC2, ProNaCC4, and In this category, a total of 39 genes were identified and, selected for their involvement as transcriptional regulators of staphylococcal virulence factors the genes agrA, garL, graS_1, graS_2, lytR, pyrR, saeR, saeS, sarA, sarR, sarS, sarV, sarX, sarZ, srrA, srrB tarL, vraR, vraS walk, walR, ypdA, yycHm and yycI were detected in all the genomes.
Enterotoxin Insight Further analyses were performed for the gene coding staphylococcal enterotoxins. A novel database was built for ABRicate starting from amino acid sequences extracted from GenBank nucleotide and uniprot databases  [72]. The sequences manually annotated (Swiss-prot) and un-reviewed enterotoxin variants (TrEMBL) were translated with the EMBOSS tool [73] and combined to build the final database (Supplementary Material Database S3), which included, in total, 99 sequences of staphylococcal enterotoxins and their variants. The strain ProNaCC1 was positive for the genes of the enterotoxins SED, SElJ, SER, SelX, and SElW; the strains ProNaCC2, ProNaCC4, and ProNaCC6 were positive for both the gene of the enterotoxins SElX and SElW and for the genes of SEE, SEA, and SEH respectively.
The strain ProNaCC5 was positive for the genes of SEY, SElX, SEG, SEN, SElU, SEI, SEM, and SEO, and the strain ProNaCC7 was positive for the gene coding for SED, SElJ, SER, SEA, SelW, and SElX (Table  6 ). The genomes were submitted to the PHASTER web server for the identification of prophage sequences and the output was tested for the presence of SEs with the database developed for ABRicate. The intact prophages regions and the hypothetical proteins sequences associated with prophages were considered in the study for the analysis of SEs similarities. According to the score assigned by the software, 20 of the prophages most common regions within the six enterotoxigenic genomes were identified: three regions resulted intact, 13 regions were incomplete, and four were questionable (Table  S4 ).