dataesr/scientific-paragraphs-categorization huggingface dataset:
Results 2. 1. Prevalence and Description of ARGs Found in the New Bacterial Species from the Human Microbiota When using the reciprocal BLAST program, we found a total of 883 ARGs encoding for enzymes that may be involved in the development of antimicrobial resistance to the nine studied antibiotic classes (Figure  1 ). Details of all ARGs by classes found per bacteria is provided in (Supplementary Table  S1 ). From the 335 studied bacteria, 278 have at least one ARG and 57 bacteria did not harbor ARGs (29 Firmicutes, 15 Actinobacteria, 7 Bacteroidetes/Chlorobi group, 3 Proteobacteria and 2 Fusobacteria) (Supplementary Table  S1 ). The Firmicutes contained 654 ARGs that were found in 173 of the 203 analyzed genomes; they encode for ARGs that belong to the nine antibiotic classes. The most frequent ARGs were for Macrolides Lincosamides Streptogramins (MLS) (n = 212;32%), aminoglycosides (n = 161;25%) and tetracyclines (n = 100;15%) resistance genes (Figure  1 ).
Prevalence and Description of ARGs found in the New Bacterial Species from the Human Microbiota When using the reciprocal BLAST program, we found a total of 883 ARGs encoding for enzymes that may be involved in the development of antimicrobial resistance to the nine studied antibiotic classes (Figure  1 ). Details of all ARGs by classes found per bacteria is provided in (Supplementary Table  S1 ). From the 335 studied bacteria, 278 have at least one ARG and 57 bacteria did not harbor ARGs (29 Firmicutes, 15 Actinobacteria, 7 Bacteroidetes/Chlorobi group, 3 Proteobacteria and 2 Fusobacteria) (Supplementary Table  S1 ). The Firmicutes contained 654 ARGs that were found in 173 of the 203 analyzed genomes; they encode for ARGs that belong to the nine antibiotic classes. The most frequent ARGs were for Macrolides Lincosamides Streptogramins (MLS) (n = 212;32%), aminoglycosides (n = 161;25%) and tetracyclines (n = 100;15%) resistance genes (Figure  1 ).
The protein sequences of ARGs potentially conferring β-lactam resistance were aligned with the already known β-lactamase domains in the Conserved Domain Database (CDD)  [29]. We searched for motifs, i. e. , sequence patterns that occur repeatedly, in these proteins by MAST scanning within the MEME suite  [30]. This structural and functional analysis showed that 44 sequences conserved at least two amino-acid motifs, allowing its classification as the serine β-lactamase, including 41 in class A, two in class C and one in class D (Table  2 ). We identified fourteen ARGs encoding for potential amino acid bonds to zinc molecules that can be considered as metallo-β-lactamase, including six in the subclass B1, seven in the subclass B3 and one with the HARLDQ motif (Figure  2 , Table  2 ). The phylogenetic tree based on the analysis of β-lactams-resistance proteins showed that fifteen proteins from twelve species represented a single line descent together with already known ARGs for β-lactams, including the Cfx, CepA, CGA, CIA, Bla and BCL in class A, to IND in class B, to Amp in class C and OXA in class D.
Transposases were found in the genomes of thirteen species that formed a robust phylogenetic group with pathogenic bacteria in the ARG protein-based trees (Table  3 ). These mobile genetic elements were located in the region coding for the ARGs (not more than ten genes up and down the gene) which suggests possible gene acquisition by transposition. Thus, a gene cluster surrounded by an integrase and a transposase containing a tet(M) gene (Figure  6 ) found on Staphylococcus aureus ST398 (Accession number: AM990992. 1) chromosome was found in eight Gram+ bacteria analyzed, five of which were isolated from the gut, two from the vagina and one from urine (Anaerococcus jeddahensis, Anaerococcus mediterraneensis, Ndiopella massiliensis, Massiliomicrobiota massiliensis, Polynesia massiliensis, Peptoniphilus colimassiliensis, Khoudiadiopa massiliensis and Urinacoccus massiliensis) (Figure  6  and Figure  S8 , Table  3 ).
Genome were isolated from the gut, two from the vagina and one from urine (Anaerococcus jeddahensis, Anaerococcus mediterraneensis, Ndiopella massiliensis, Massiliomicrobiota massiliensis, Polynesia massiliensis, Peptoniphilus colimassiliensis, Khoudiadiopa massiliensis and Urinacoccus massiliensis) (Figure  6  and S8, Table  3 ). A plasmid search in these bacteria showed that seven bacteria had a repUS43 plasmid on the same contig as the tet(M) gene. The S. aureus genome showed an origin of replication CAQ49392. 1 described as transcriptional regulator, Cro/CI family with a replication initiation domain (Figure  6 ). [33]  between the genes in the genetic environments studied. Another example includes Bacteroides cutis, Butyricimonas phoceensis and Prevotella lascolaii that harbored two genes, ermF and tetX, next to each other that were found to be similar to genes on a plasmid from Bacteroides fragilis strain FDAARGOS (Accession number CP054002) (Figure  7 ).
Likewise, Psychrobacter timonensis has two aminoglycoside resistance genes, aadA1 and aac(2 ), next to transposases (Figure  8 ). The aadA1 gene has already been characterized in an Int1 integron of Pseudomonas aeruginosa (Genbank: AJ584652. 2) isolated from a lower respiratory tract infection and the following aac gene in an Escherichia coli MS 175-1 from the human gut. The two genes aadA1 and aac2 are surrounded by an insertion sequence and two tyrosine recombinases in Psychrobacter massiliensis and E. coli genomes.
Urinicoccus massiliensis transposase IS4351 followed by tet(X) and erm(F) in Butyricimonas phoceensis genome, whereas Prevotella lascolaii has in 5′ an insertion sequence (IS1380) followed by an erm(F) and tet(X). These data suggest that these ARGs can be mobilizable. Likewise, Psychrobacter timonensis has two aminoglycoside resistance genes, aadA1 and aac(2'), next to transposases (Figure  8 ). The aadA1 gene has already been characterized in an Int1 integron of Pseudomonas aeruginosa (Genbank: AJ584652. 2) isolated from a lower respiratory tract infection and the following aac gene in an Escherichia coli MS 175-1 from the human gut. The two genes aadA1 and aac2 are surrounded by an insertion sequence and two tyrosine recombinases in Psychrobacter massiliensis and E. coli genomes.
Likewise, Sommer et al. , (2009) characterized antibiotic resistance genes in human microflora and obtained an average of 69. 5% nucleotide coverage with GenBank sequences  [38]. Altogether, these results demonstrate the great diversity of resistance mechanisms and that the genes known so far in the bibliography and referenced are only a tiny part of the repertoire of ARGs existing in nature  [39]. We found evidence of putative transfers of ARGs between the studied bacteria and pathogenic bacteria. Phylogenetic trees with known ARGs and the hits resulting from ARG blasts highlighted some potential HGTs with pathogenic bacteria. The difference between the GC content of ARGs and their host genome indicate possible lateral acquisition. We estimated the transferability of these ARGs and their mechanism of integration into these new bacterial species genomes through the presence of transposon and plasmids.
Materials In this study, we analyzed the genome sequences of 335 new bacterial species that have been isolated from human microbiota for the first time in the Institut Hospitalo-Universitaire Méditerranée Infection using culturomics  [28]  (Supplementary Table  S1 ). The studied bacteria were identified in different biological specimens, specifically in stool (76%), vaginal swabs (9%), urine (4%), skin swabs (3%) and sputum and nasal swab (3%) samples. They were distributed as follows: 203 Firmicutes, 73 Actinobacteria, 43 Bacteroidetes, 14 Proteobacteria and 2 Fusobacteria. The genomes of these bacteria were previously sequenced using different whole genome sequencing (WGS) strategies and reads were assembled by the IHU bioinformatic team. Genome sequences have been deposited in the database of the European Bioinformatics Institute (EMBL-EBI) https://www.
Computational Method for Predicting ARGs Genome-based ARG tracking was performed by using a BLAST-based approach to query input amino acid sequence data from the 335 bacteria for the presence of a pre-determined set of ARG determinants contained in AMR, CARD and ARG-ANNOT reference databases  [51] [52] [53]. The output was sorted to keep the hits corresponding to enzymatic mechanisms of resistance, i. e. , inactivating genes or target-site alteration. Hits obtained with a minimum threshold of 50% identity and 70% coverage were considered as significant. Protein sequences corresponding to the obtained hits from this first BlastP were individually compared with all proteins in each bacterial genome by reciprocal BlastP. The best BLAST for each protein was extracted using a Perl script. The best hits for each resistance gene-bacteria and bacteria-resistance gene pairs were compared in order to determine the number of reciprocal best hits for each pairwise comparison.
Phylogenetic Analysis The potential ARGs were used to construct individual protein trees. The amino acid sequences were retrieved from the different databases, aligned using Clustal  [54]  and further optimized visually. Maximum likelihood (ML) analyses were conducted in the MEGA v7 program with default settings  [55]. Clade robustness was assessed using a bootstrap analysis with 1000 replicates  [56]. We used Bioedit  [57]  to calculate the identity matrix of the alignments of fosfomycin and β-lactamase amino acid sequences. This allowed us to determine the rate of similarity between well-described reference sequences in order to set the parameters for defining new ARGs or new variants. A clustering of two or more new ARGs can be considered as putative new resistance gene families.
Bioinformatic Characterization of Conserved Protein Domains and Motifs The conserved domain database, CDD  [29]  (https://www. ncbi. nlm. nih. gov/Structure/ cdd/cdd. shtml accessed on 15 March 2021), was used to find the protein domains in order to characterize the enzyme functional class. The MEME/MAST Suite was used to identify the patterns/motifs which are specific to a particular family or subfamily of β-lactamases  [30]  (http://meme. nbcr. net accessed on 15 May 2021). The conserved patterns were derived using the pattern-search tool with default parameters of fingerprint width of 10 residues and minimum occurrence of two patterns per sequence. Motif scanning was performed by the Motif Alignment and Search Tool (MAST)  [30]. The determination of family specific patterns/motifs helps to assign the newly identified β-lactamases to one or the other family or subfamily; it can be very useful to reveal gene novelty.
Detection of Mobile Genetic Elements Associated with ARGs and Determination of GC Content Further in silico analysis of the genetic context of the predicted ARGs was performed to search for the presence of mobile elements. Transposases were searched up to ten genes upstream and downstream for the predicted ARGs using a locally developed program that considered the annotation of the genome to locate the positions of potential mobile genetic elements. The plasmids were searched within the PlasmidFinder database  [58]  with a minimum threshold of 90% identity and 90% coverage. The sequence of the ARGs and its corresponding genomes in FASTA format were used to calculate the GC%. The start and end codons were not deleted from the sequence for all the tested genes. Comparative analysis between the studied genomes containing ARGs and genomes of pathogenic bacteria with homologs to these ARGs was generated and drawn using Easyfig software  [33].
1C-D ). Nonmalignant B cells clustered depending on their GC, Mem or PC type and expressed polyclonal BCR (Fig. 1C-D, S2 ). We also performed high-throughput droplet-based 10x 3' (Fig. 1E-G ) and 10x 5' scRNA-seq (Fig. 1H-J ) on all live cells from FL cell suspensions. We first analyzed the 10x 3' and 10x 5' datasets separately to avoid technology-related batch effects. We defined TME non-B cells, malignant and non-malignant B cells and visualized them in UMAP. Non-malignant B cells and TME cells did not cluster on the basis of donor origin (Fig. 1F-G, I-J ). Conversely, malignant FL B cells displayed strong inter-patient heterogeneity and only paired diagnosis/relapse samples (FL068 and FL338) mapped together (Fig. 1I-J ). Thus, malignant FL B cells exhibit patient-specific gene expression profiles distinct from non-malignant B cells.
Non-malignant B cell signatures Beyond inter-patient divergence, we sought to identify patterns of intra-patient gene expression heterogeneity in malignant B cells. To characterize potential similarities between malignant and non-malignant B cells, we first computed four lists of differentially expressed genes (DEGs) associated with PC, non-PC, GC and Mem gene expression patterns (detailed in Computational Methods and Fig. S3A-F , Supplementary Data Table  1 ). Then for each cell in a scRNA-seq dataset, we reported the expression of genes from a given signature with a single value defined as the signature score. That scoring strategy accurately classified normal GC, Mem and PC subsets analyzed by FB5P-seq (Fig. S3G-H ) or by 10x 3' and 10x 5' scRNA-seq (Fig. S3I-J ). We then used those signatures to detect malignant B cells expressing gene expression programs close to those of non-malignant B cell types.
PC-like FL B cells FL B cells are not expected to resemble antibody-producing PC in most cases, but some cases of FL with plasmacytic differentiation have been described  41. We first intended to rule out plasmacytic differentiation as a source of transcriptional heterogeneity in malignant B cells. We computed PC and non-PC signature scores for malignant FL B cells in the FB5P-seq, 10x 3' and 10x 5' datasets. Almost all FL B cells were transcriptionally distinct from PC cells (Fig. 2A-B ). However, we detected rare FL B cells carrying a PC-like gene expression profile in all datasets. PC-like cell frequency reached 2-3% of malignant FL B cells in 2/15 patients (FL2 and FL539), but those cells were absent or present in frequencies <0. 3% in other samples (Fig. 2C , Table  S7 ). PC-like cells had low expression of MS4A1 and high expression of PC markers CD27, XBP1, MZB1 and PRDM1 (Fig.
Main source of intra-tumoral transcriptional heterogeneity We then investigated intra-sample heterogeneity in each sample individually. We identified the GC-like to Mem-like gene expression continuum in every FL patient. In UMAP embeddings of sample-specific malignant B cells gene expression datasets, ScoreGCtoMem values varied continuously from one end of the cloud of cells to another (Fig. 4A ). That observation suggested that the GC-like to Mem-like continuum was the main axis of transcriptional heterogeneity within FL B cells. In principal component analysis (PCA), the first principal component (PC1) cumulates the largest proportion of gene expression variance and represents the main axis of gene expression heterogeneity in a scRNA-seq dataset. In 13/15 samples, PC1 scores were tightly correlated with ScoreGCtoMem values (Fig.
4B-C ). Unsupervised pseudo-time (PT) analysis, such as the one performed in the Monocle2 algorithm  43  , may identify cellular trajectories representing the most likely axes of variation (interpreted as differentiation axes) in a scRNA-seq dataset. When analyzing malignant B cells of each sample individually, PT ordering was also tightly correlated to the ScoreGCtoMem continuum for 12/15 FL samples (Fig. 4D-E ). Together, PCA and unsupervised PT analyses clearly demonstrated that the GC-like to Mem-like continuum was the main source of intra-tumoral transcriptional heterogeneity in malignant B cells.
Impact of TFH-derived activation signals Interactions with TME may drive FL B cell plasticity. Non-supervised clustering of TME cells from the integrated 10x 3' and 10x 5' datasets defined 15 distinct clusters (Fig. 6A, S6A ). Parallel analysis of TME cell types by flow cytometry (Fig. S1D ) demonstrated tight correlations between cell subset frequencies in scRNAseq and FACS (Fig. S6B ). The TME composition was variable and characterized by the predominance of CD4 + T cells, with high frequencies of T central memory (TCM), TFH and T regulatory (TREG) subsets (Fig. 6B, S6C ). TFH cells were the main producers of CD40LG, IL4 and IL21 mRNA (Fig. 6C, S6D ), encoding critical positive regulators of B cell activation, proliferation and differentiation. TFH cells control normal B cell differentiation at multiple stages of GC responses  45.
GC-like cells cultured with TFH-mimicking signals, but not autologous TFH cells, also showed a lower ScoreGCtoMem. Thus, TFH-derived signals impacted B cell transcriptional states in vitro, "pushing" malignant B cells towards a GC-like state. To track whether such dynamics also occurred in vivo, we defined a set of 25 genes whose expression was robustly increased in malignant B cells upon co-culture with autologous TFH in vitro (Fig. S8B-C ). We computed the corresponding signature score (ScoreTFH-activated) in malignant B cells from 10x 3' and 10x 5' datasets. In 10/15 samples, a minor subset of malignant B cells with high ScoreTFH-activated values clustered together in sample-specific UMAP (Fig. 6F ), suggesting those cells had been recently activated by TFH cells in vivo. Unsupervised marker genes-based metaclustering (detailed in Computational Methods) of the 74 sample-specific malignant B cell clusters from the 15 scRNA-seq datasets identified 6 distinct metaclusters (Fig.
7B, S9A ): centrofollicular (CF), perifollicular (PF), interfollicular (IF) and extrafollicular (EF). ST slides capture mRNA on 55-µm diameter spots, thereby averaging the transcriptome of several cells (5 to 20 cells in lymphoid tissues). We constructed a reference scRNA-seq dataset of FL cells by integrating our full 10x 3' and 10x 5' datasets (Fig. 7C ), and applied reference-based deconvolution to infer the presence of malignant and non-malignant cell subsets in each spot  49. The deconvolution algorithm cell2location outputs q05_spot_factors, probabilistic values reflecting the number of cells of a given cell type in ST spots  49. For most of the cell types/states assessed, there was a moderate to good correlation between q05_spot_factors and signature scores based on marker genes (Fig.
DISCUSSION Tumors are dynamic ecosystems where malignant and non-malignant cell types co-exist. Here we have combined supervised and non-supervised analyses of integrative single-cell transcriptomics datasets to uncover the conserved states of malignant B cells in FL. Beyond inter-patient divergence  37, 38  , our work reveals five recurrent malignant B cell states: GC-like, Mid, Mem-like, TFH-activated and PC-like states. Our analyses support a model where FL B cell states span a continuum from quiescent inter-and peri-follicular Mem-like to proliferative centro-follicular GC-like state, with most malignant B cells being in a quiescent nonactivated Mid state within follicles. In that model, activation of Mem-like and Mid cells through CD40dependent interactions with TFH cells, which may occur in peri-follicular FDC-rich areas, induces a specific TFH-activated state which likely serves as an intermediate towards the GC-like state.
Cells were then stained with intracellular marker fluorophore-conjugated antibodies (BCL6 and Ki67, Table  S3 ) for 30 min on ice. After a final wash in FACS Buffer, cells were resuspended in FACS Buffer and samples were acquired on a BD LSR Fortessa or a BD FACS Influx. Gating strategies for analyses and cell sorting are described in Fig. S1 , S4, S7. Flow cytometry data were analyzed with FlowJo software.
FB5P-seq We used a custom bioinformatics pipeline to process fastq files and generate single-cell gene expression matrices and BCR sequence files as previously described  39. Detailed instructions for running the FB5P-seq bioinformatics pipeline can be found at https://github. com/MilpiedLab/FB5P-seq. Quality control was performed on each dataset independently to remove poor quality cells. Cells with less than 250 genes detected were removed. We further excluded cells with values below 3 median absolute deviations (MADs) from the median for UMI counts, for the number of genes detected, or for ERCC accuracy, and cells with values above 3 MADs from the median for ERCC transcript percentage. For each cell, gene expression UMI count values were log-normalized with Seurat NormalizeData with a scale factor of 10,000  73.
We further excluded bad quality cells with values below 2 median absolute deviations (MADs) from the median for UMI counts or for the number of genes detected, and cells with values above 1. 5 MADs for mitochondrial transcript percentage for most datasets. The threshold for percentage of mitochondrial transcripts was adjusted to 2 MADs for the dataset including hashtags. HTO barcodes for sample demultiplexing after hashing were counted using CITE-seq-count and were normalized for each cell using a centered log ratio (CLR) transformation across cells implemented in the Seurat function NormalizeData. Cells were demultiplexed using MULTIseqDemux function  75  and barcodes assigned as doublets or negative were excluded from further analysis. The resulting filtered UMI count matrices were log-normalized with Seurat NormalizeData with a scale factor of 10,000.
10x 5' Raw fastq files from gene expression libraries were processed using Cell Ranger software (v3. 0. 1), with alignment on the GRCh38 reference genome. For each experiment, we retained cell barcodes if they had at least 500 genes with non-zero expression and genes if they were detected in at least 3 cells. We further excluded bad quality cells with values below 2 median absolute deviations (MADs) from the median for UMI counts or for the number of genes detected, and cells with values above 1. 5-2 MADs for mitochondrial transcript percentage. HTO barcodes for sample demultiplexing after hashing were counted using CITE-seqcount and were normalized for each cell using a centered log ratio (CLR) transformation across cells implemented in the Seurat function NormalizeData. For experiments where samples were tagged with a single hashtag barcode, cells were demultiplexed using MULTIseqDemux function and barcodes assigned as doublets or negative were excluded from further analysis.
For the experiment where samples were tagged with combinations of 3 hashtags (Table  S5  and Fig. 6D-E ), we manually set thresholds for hashtag positivity based on visual inspection of RidgePlots for all normalized hashtag signals, and selected cells which were positive for triplet combinations of hashtags used to label the different culture conditions. All other cells were excluded from further analysis. The resulting filtered UMI count matrices were log-normalized with Seurat NormalizeData with a scale factor of 10,000. BCR-seq raw fastq files were processed with the FB5P-seq pipeline  39  as described above for FB5P-seq datasets, omitting the part of the pipeline related to gene expression analysis, and using the list of cellassociated 10x barcodes from CellRanger analysis as inputs for splitting bam files upstream Trinity assembly of BCR contigs.
Segregation of TME, malignant and non-malignant B cells For each sample in the 10x 3' and 10x 5' datasets, B cells were distinguished from TME non-B cells by applying low resolution graph-based cell clustering (Seurat FindNeighbors and FindClusters, resolution 0. 1) and scoring B cell marker genes (CD19, CD79A, CD79B, MS4A1) with Seurat AddModuleScore. Clusters with a high B cell expression score were defined as B cells, other clusters were defined as TME non-B cells. In B cell clusters, non-malignant and malignant B cells were then distinguished by light chain isotype restriction. Briefly, cells were assigned as IGKC if the normalized expression of IGKC was superior to the sum of the normalized expressions of IGLC genes (IGLC [2-7]), or IGLC if the summed expression of IGLC genes was superior to the expression of IGKC.