dataesr/scientific-paragraphs-categorization huggingface dataset:
The targeted cell input was 10,000 cells per condition. The Chromium Next GEM (Gel Bead-In Emulsion) Single Cell 3' Gel beads v3. 1 kit (10X Genomics, Pleasanton, CA, USA) was used to create GEMs following manufacturer's instruction. All samples and reagents were prepared and loaded into the chip and ran in the Chromium Controller for GEM generation and barcoding. GEMs generated were used for cDNA synthesis and library preparation using the Chromium Single Cell 3' Library Kit v3. 1 (10X Genomics) following the manufacturer's instruction. Generated libraries were sequenced on NovaSeq 6000 system using HiSeq 100 base pair reads and dual indexing. Cells were sequenced to an average depth of 31,383 reads per cell. The human genome, Ensembl GRCh38. 98. gtf, and the SARS-CoV-2 genome, NCBI Genome database accession MT020880.
Cell type annotation We used the standard scRNA-seq analysis pipeline for clustering  [41]. Briefly, to account for transcript dropout inherent to scRNA-seq, we removed genes that were expressed in fewer than 3 cells and removed cells that expressed fewer than 200 genes. Next, we filter out cells with more than 10 percent of mitochondrial genes. We did not find a correlation between viral copy number and mitochondrial expression. The resulting raw unique molecular identifier (UMI) counts in each cell were normalized to their library size. Then, normalized counts were square-root transformed, which is similar to a log transform but does not require addition of a pseudo count. Data pre-processing was performed in Python (version 3. 7. 4) using Scanpy (version 1. 4. 6) [42]. We visually observed batch effects between conditions in 2-dimensional cellular embeddings.
This value was determined by a threshold of viral counts in the mock condition. While the mock condition is not expected to have viral counts, we did observe a small number that we attribute to misalignment. We observed only 5 mock cells with full SARS-CoV-2 viral genome transcript counts ≥ 10 transcripts. These criteria allowed us to find 144 infected cells at 1 dpi, 1428 cells at 2 dpi and 3173 cells at 3 dpi. To quantify the extent to which an individual cell is transcriptionally similar to an infected cell, we used a previously developed graph signal processing approach called Manifold Enhancement of Latent Dimensions (MELD)  [51]. We encoded a raw experimental score for each cell in the dataset such that -1 represents a bystander or uninfected cell and +1 represents an infected cell. Using the kernel from the BB-kNN graph (described above), these raw scores were smoothed in the graph domain, yielding an "infection score" per cell that represents the extent to which an individual cell is transcriptionally similar to infected cells.
Viral genome read coverage analysis To visualize the viral read coverage along the viral genome we used the 10X Genomics cellranger barcoded binary alignment map (BAM) files for every sample. We filtered the BAM files to only retain reads mapping to the viral genome using the bedtools intersect tool  [52]. We converted the BAM files into sequence alignment map (SAM) files in order to filter out cells that were removed in our single cell data preprocessing pipeline. The sequencing depth for each base position was calculated using samtools count. To characterize read distribution along the viral genome we counted transcripts of 10 different ORFs: ORF1ab, Surface glycoprotein (S), ORF3a, Envelope protein (E), Membrane glycoprotein (M), ORF6, ORF7a, ORF8, Nucleocapsid phosphoprotein (N) and ORF10.
CONDENSE overview and notation The CONDENSE algorithm for analyzing gene dynamics in single-cell data is as follows: 1. Obtain the condition of interest using pseudotime or other trajectory inference method 2. Obtain a joint density estimate over the cells for each gene as a function of time using kernel density estimation 3. Normalize the per gene density estimate by its temporal condition to obtain a conditional density embedding of genes 4. Use the gene dynamics embeddings for downstream analyses such as visualization, clustering, archetypal analysis, and gene-set enrichment While the method was designed with single-cell expression data in mind, CONDENSE can be applied to any dataset of m examples X = {x (1) ,. , x (m) } with a condition c. Here, we treat each gene as a sample, with expression of that gene across the scRNA-seq dataset represented by a vector x (i) for a particular gene i.
Archetypal analysis of CONDENSE matrix To find extremal dynamics on an individual gene basis between strata, we performed archetypal analysis on the conditional density estimate matrix, P S k. Archetypal analysis estimates the principal convex hull of a d-simplex dataset  [59]. The convex hull of d dimensions can be thought of as an envelope around the data. We implement archetypal analysis by defining the principal convex hull as the (d -1)-dimensional convex hull that best accounts for the data according to, arg min Q,Y D(P S k | P S k QY) ( 7 ) where D is a measure of distortion and Q and Y are the optimal matrices that minimize the distortion D subject to the constraints, | q d | 1 = 1 and Q ≥ 0, and | y m | 1 = 1 and Y ≥ 0, respectively. These constraints yield archetypes, A = P S k Q (8) which are convex combinations that can be approximated by a weighted average of the feature vectors P S k ∈ R cxd  [60].
Biological insights from CONDENSE After identifying distinct gene dynamics by archetypal analysis, we further characterize the set of genes expressed in our dataset by visualizing the multidimensional, conditional density estimates in a single dimensionally-reduced embedding. Any manifold learning or dimesionality-reduction method is possible but we use the popular UMAP algorithm to learn a mapping u : R c×k -→ R 2  [44]. We fit this function on the conditional density embedding, P S k and transform P S k to get a conditional density embedding, P S k in UMAP space. In order to characterize the dynamics according to archetypal or extremal dynamics within the HBEC dataset, we assign each gene to an archetype by minimizing the Euclidean distance of each gene's conditional density embedded point to the nearest archetypal point in UMAP space, a i = arg min (A T ) P S k -(A T ) 2 2 (9) where a i represents the archetype that satisfies the minimization problem, (A T ) represents the archetypes in UMAP coordinates and P S k denotes the conditional density embedding of genes per strata k in UMAP coordinates.
Code and Data Availability All differential gene expression analyses and their associated metrics are publicly available at the Van Dijk Lab GitHub: https://github. com/vandijklab. Software for the CONDENSE method and scripts used to analyze the scRNA-seq data are also available at the Van Dijk Lab GitHub: https://github. com/vandijklab. The annotated scRNA-seq data can be browsed with an interactive web-tool, courtesy of the Chan-Zuckerberg Initiative. The raw scRNA-seq data will also be made publicly available and deposited in NCBI's GEO database.
Acknowledgments We thank the  Yale Center for Genome Analysis  for help with sequencing and alignment. This work was supported by  NIH  grants  K08 AI128043  (CBW), R21AI133440 (AW),  R01AI141609  (AW), and  K08AI119139  (EEF). CBW was also supported by a  Burroughs Wellcome Fund Career Award for Medical Scientists and Robert Leet and Clara Guthrie Patterson Trust  Award. CBW and DvD were supported by a  Fast Grant (Emergent Ventures ). AI and AMP are Investigators of the  Howard Hughes Medical Institute. Author contributions: conceptualization, MMA, JW, NGR, VG, EFF, DvD, CBW; methodology, MMA, NGR, EFF, DvD, CBW; computational analysis NGR, VG, DvD; validation, MMA, CBW; formal analysis NGR, VG, DvD; investigation, MMA, RBF, JW, CBW; resources, BW, EFF; writing-original Draft, MMA, NGR, VG, DvD, CW; writing-review & editing, all authors; visualization, NGR, VG, DvD; supervision, DvD, CBW; funding acquisition, DvD, CBW.
Methods Data merging and SNP filtering. We combined Illumina BovineSNP50 BeadChip genotypes of 57 Brune de l' Atlas individuals (TUNIND) sampled from our previously published data  4, 52  with data already available for 221 animals belonging to seven North African populations (BAL, BIS, CHE, CHF, GUE, TID and OUL) obtained from Flori et al. 16  and Gautier et al. 53. We also included genotyping data belonging to 9 other populations, representatives of European taurines (EUT) (four breeds: ANG, HOL, JER and MON), African taurines (AFT) (three N'Dama populations: ND1, ND2 and NDA) and indicine (two populations: GIR and NEL) from Matukumalli et al. 54. All genotypes were recovered from the web-interfaced genetic Diversity Exploration (WIDDE) database  55. We performed a relatedness test between individuals within each population using PLINK  56.
First, we excluded rare SNPs with low minor allele frequencies (MAF) < 0. 05. Then, the whole genotype dataset was subjected to linkage disequilibrium (LD) pruning using the default parameters of PLINK (SNP window size:50, step 5 SNPs, r 2 : 0. 5). In total, 38,464 SNPs spread over all autosomal chromosomes were finally considered for population structure analyses. Population structure and genetic relationship analyses. Population structure was inferred by PCA for African, European, indicine and North African populations using the adegenet R package  57. Unsupervised hierarchical clustering was carried out for all populations using ADMIXTURE 1. 23 software  58. We ran ADMIX-TURE with cross-validation for values of K from 2 through 17 (the number of populations) to identify the best value of K clusters.
DISTRUCT software  59  was then used to graphically display ancestry within each individual. The pairwise fixation index (F ST ) between populations was estimated using Genepop 4. 6 software 60. The patterns of population splits and mixtures were inferred using TreeMix  61. First, we built a maximum likelihood tree of the 17 populations of the study with no migration events allowed and setting GIR as outgroup. Then, we built a phylogenetic tree of these populations and started adding migration events (modeled as edges) sequentially to the phylogenetic model. The migration edges were added until 99. 93% of the variance in ancestry between populations was explained by the model. The residuals from the fit of the model to the data were visualized using the R script implemented in TreeMix.
Outliers are those loci that require the locus-specific component to explain observed genetic diversity. For the Markov chain Monte Carlo (MCMC) algorithm we used 20 pilot runs of 5,000 iterations, a burn-in of 50,000 iterations, a thinning interval of 10 (5,000 iterations were used for the estimation of posterior odds) with a resulting total number of 100,000 iterations. To control the number of false positives, significant SNPs were defined by applying a q-value threshold of 0. 05. Haplotype extended patterns were investigated using three metrics implemented in rehh package  63  : the iHS within-population statistic  64  and two between-population methods: Rsb  65  and XP-EHH  66. In iHS computation, the information on the ancestral and derived allele state is needed for each SNP because this statistic is based on the ratio of the EHH associated to each allele.
We used 1-Mb sliding windows that partially overlapped 10 kb with adjacent windows to perform selection signature detection. A window is classified as putatively under selection when it contains at least 3, 4 and 4 markers exceeding the significance threshold of -log10 (p value) = 3 for iHS, Rsb and XP-EHH tests, respectively. Finally, we checked the overlap of the candidate genomic regions detected with at least two EHH-based approaches with the previously identified bovine Quantitative Trait Loci (QTL) available in the cattle QTL database (https ://www. anima lgeno me. org/cgi-bin/QTLdb /BT/index ). The overlaps were checked using QTL coordinates according to the Bos taurus genome assembly ARS-UCD1. 2.
Gene identification and functional enrichment analysis. Candidate genome region intervals detected by at least two EHH-based methods (iHS, Rsb, XP-EHH) were interrogated for genes annotated to the Bos taurus genome assembly ARS-UCD1. 2 using BioMart tool of Ensembl (https ://www. ensem bl. org/bioma rt/ martv iew/c8fe3 a6996 1a408 8a55b 7a249 db7e2 fa). Cattle structural variants which overlapped the genomic coordinates (in bp) of these relevant candidate selective sweep regions were retrieved using the same database. We have only considered structural variants of less than 8 Mb which corresponds to the maximum size that can be identified, from whole genome sequence data, by the pindel software (https ://gmt. genom e. wustl. edu/packa ges/ pinde l/user-manua l. html). We used the online tool, Database for Annotation, Visualization and Integrated Discovery (DAVID) software version 6.
Outlierwindows from iHS, Rsb and XP-EHH tests include 57, 581 and 305 known genes, respectively (Table1). Genes identified with Rsb and XP-EHH are distributed as follows: 151 and121, 264 and 127, 166 and 57 for AFT/ North African, EUT/North African, IND/North African comparisons, respectively (Table1). Thirty genes were common to both iHS and EUT/North African comparison (either Rsb or XP-EHH). Similarly, 109, 143 and 65 genes were jointly identified by Rsb and XP-EHH for each of the AFT/North African, EUT/North African, IND/North African comparisons, respectively of which 74, 97 and 50, respectively, could be mapped by DAVID Bioinformatics resources (https ://david. ncifcrf. gov/). Gene Ontology (GO) analysis showed that AIG1 (IPR006703, n = 6, Benjamini-corrected p value = 4. 45 × 10 -7 ) and P-loop containing nucleoside triphosphate hydrolase (IPR027417, n = 14, Benjamini-corrected p value = 0.
Acknowledgements This study was supported by grants from  International Foundation for Science (IFS  grant  B/5478 ), and by the  Basic Science Research Program  through the  National Research Foundation of Korea (NRF) , funded by the  Ministry of Science, ICT & Future Planning  ( 2013K2A4A1044519 ). The authors would like to thank  Dominique Rocha  and  Emmanuelle Rebours  ( GABI unit ) for their help in data acquisition and genotyping.
19 47,120,000 49,140,000 2. 02 MARCHF10, TANC2, CYB561, ACE, KCNH6, DCAF7, TACO1, MAP3K3, LIMD2, STRADA, CCDC47, DDX42, FTSJ3, PSMC5, SMARCD2, bta-mir-10173, TCAM1, GH1, CD79B, SCN4A, ERN1, SNORD104, SNORA50C, TEX2, PECAM1, MILR1, POLG2, DDX5, bta-mir-3064, CEP95, SMURF2, KPNA2, C17orf58, BPTF Test BTA Start (pb) End (pb) length (Mb) Genes 12 28,400,000 30,490,000 2. 09 PDS5B, N4BP2L2, N4BP2L1, BRCA2, ZAR1L, FRY, RXFP2, bta-mir-2299, B3GLCT, HSPH1, TEX26, MEDAG, ALOX5AP, USPL1, HMGB1, KATNAL1 ATP2C2, MEAK7, COTL1, KLHL36, USP10, CRISPLD2, ZDHHC7, KIAA0513, FAM92B, GSE1, GINS2, EMC8, COX4I1, IRF8, FOXF1, MTHFSD, FOXC2, FOXL1, FBXO31, MAP1LC3B, ZCCHC14, JPH3, KLHDC4, SLC7A5, CA5A, BANP, ZNF469, 18 10,590,000 15,060,000 4. 47 ZFPM1, ZC3H18, CYBA, MVD, SNAI3, CTU2, RNF166, PIEZO1, bta-mir-2327, CDT1, APRT, GALNS, TRAPPC2L, CBFA2T3, ACSF3, CDH15, SLC22A31, ANKRD11, SPG7, RPL13, SNORD68, CPNE7, DPEP1, CHMP1A, CDK10, SPATA2L, VPS9D1, ZNF276, FANCA, SPIRE2, TCF25, MC1R, TUBB3, DEF8, DBNDD1, GAS8, Rsb IND versus North AFT U1, SHCBP1, VPS35, ORC6, MYLK3 MEIOC, CCDC43, DBF4B, ADAM11, GJC1, HIGD1B, EFTUD2, bta-mir-2343, 19 44,340,000 45,970,000 1.
56 TPGS1, CDC34, GZMM, BSG, HCN2, POLRMT, FGF22, RNF126, FSTL3, PRSS57, PALM, MISP, PTBP1, PLPPR3, AZU1, PRTN3, ELANE, CFD, MED16, R3HDM4, KISS1R, ARID3A, WDR18, GRIN3B, TMEM259, CNN2, ABCA7, ARHGAP45, POLR2E, GPX4, SBNO2 XP-EHH EUT versus North AFT 15 33,830,000 35,760,000 1. 93 GRAMD1B, SCN3B, ZNF202, SAAL1, TPH1, SERGEF, KCNC1, MYOD1, OTOG, USH1C, ABCC8, KCNJ11, NCR3LG1, NUCB2, PIK3C2A, RPS13, SNORD14, PLE- KHA7, U6, C15H11orf58 20 57,910,000 59,770,000 1. 86 U6, ANKH, OTULIN, OTULINL, TRIO, DNAH5 21 14,800,000 16,630,000 1. 83 SLCO3A1, SV2B, U6, AKAP13, KLHL25 12 28,400,000 30,490,000 2. 09 PDS5B, N4BP2L2, N4BP2L1, BRCA2, ZAR1L, FRY, RXFP2, bta-mir-2299, B3GLCT, HSPH1, TEX26, MEDAG, ALOX5AP, USPL1, HMGB1, KATNAL1 GSE1, GINS2, EMC8, COX4I1, IRF8, FOXF1, MTHFSD, FOXC2, FOXL1, FBXO31, XP-EHH IND versus North AFT 18 11,580,000 14,350,000 2.
and functional annotation of the genes within the candidate regions. Outlierwindows from iHS, Rsb and XP-EHH tests include 57, 581 and 305 known genes, respectively (Table1). Genes identified with Rsb and XP-EHH are distributed as follows: 151 and121, 264 and 127, 166 and 57 for AFT/ North African, EUT/North African, IND/North African comparisons, respectively (Table1). Thirty genes were common to both iHS and EUT/North African comparison (either Rsb or XP-EHH). Similarly, 109, 143 and 65 genes were jointly identified by Rsb and XP-EHH for each of the AFT/North African, EUT/North African, IND/North African comparisons, respectively of which 74, 97 and 50, respectively, could be mapped by DAVID Bioinformatics resources (https ://david. ncifcrf. gov/). Gene Ontology (GO) analysis showed that AIG1 (IPR006703, n = 6, Benjamini-corrected p value = 4.
45 × 10 -7 ) and P-loop containing nucleoside triphosphate hydrolase (IPR027417, n = 14, Benjamini-corrected p value = 0. 0031) InterPro protein functional groups were the two significantly enriched functional classes identified in the AFT/North AFT comparison (Supplementary TableS7). Sensory perception of smell (GO:0,007,608, n = 18, Benjamini-corrected p value = 3. 23 × 10 -14 ) and G-protein coupled receptor signaling pathway (GO:0,007,186, n = 19, Benjamini-corrected p value = 3. 48 × 10 -6 ) were the most enriched biological process (BP) terms identified in the EUT/North AFT comparison. Olfactory receptor activity (GO:0,004,984, n = 21, Benjamini-corrected p value = 6. 66× 10 -7.
Results Homologous sequences of mobile colistin resistance genes in available database. As a starting point, the first MCR-1. 1 protein described (NG_050417) was used as query in a BlastP analysis against the NCBI database to fish out a total of 13'658 protein hits with aa identity ranged from 30% to 100% and alignment ≥30% (Fig. 1A ). Results include all MCR-variants (from MCR-2 to MCR-9) with aa identity ranged from 30. 82% with MCR-4. 1 to 82. 66% with MCR-6. 1 (Fig. 1B ). Interestingly, reference MCR variants exhibit almost the same sequence size (1'635-bp in average) but show a significant difference in their GC content ranged from 40. 1% for MCR-4. 1 to 56. 17% for MCR-7. 1 (Fig. 1B ), likely suggesting different origins. As shown in Fig. 1A , inferred phylogenetic tree with all obtained sequences reveals a total of 494 bacterial genera (only representative genera are shown on this figure) highlighting a high diversity of PET proteins among bacteria.
Rhizome analysis of MCR-1 sequence. To investigate the origin of MCR-1 protein, this latter was split into 11 fragments (50 aa in size) and each fragment was blasted against the NCBI database to retrieve the 100 best BlastP hits. Then, phylogenetic tree analysis was assessed for each fragment to determine its origin. As shown in Fig. 1C , 8 out of the 11 MCR-1 fragments appear to be from Moraxella genus, especially Moraxella pluranimalium. Subtree analysis of each MCR variants. As seen on Fig. 1A , the phylogenetic tree analysis shows that the 9 MCR variants appear in different branches of the tree. So, subtree containing each MCR variant was retrieved for further detail analysis: MCR-1, 2, and 6. These three MCR sequences share more than 80% of aa identity (Fig. 1B ) and therefore appear clustered together on the inferred phylogenetic tree (Fig.
MCR variants in sequenced bacterial genomes. A total of 64'628 Gram-negative bacterial genomes were downloaded from the NCBI RefSeq database and include 60 bacterial genera constituted by 1'047 bacterial species (Fig. 3 ). As presented in this figure, a RpoB-based phylogenetic tree was inferred to establish their phylogenetic relationship. Only 4'214 out of the 64'628 genomes (6. 52%) were complete genomes and the remaining genomes were a Whole-genome sequences (WGS). BlastP analysis of the nine MCR sequence variants against all downloaded data set provided significant positive results for MCR-1. 1, MCR-2. 1, MCR-3. 1, MCR-4. 1, MCR-5. 1, MCR-8. 1 and MCR-9. 1 with aa identity values ranging from 91 to 100% and alignment ≥98%. A total of 1'386 MCR hits were identified in 15 out of the 60 genera analysed.
Methods Looking for homologous MCR-1 sequences from the NCBI database. Reference sequence of MCR-1 (AKF16168. 1) was used as query for BlastP analysis against the NCBI database using as threshold e-value 10e-5. All homologous sequence with identity ≥30% and alignment ≥30% were kept for further analysis. Protein alignments were performed using Mafft tool (https://mafft. cbrc. jp/alignment/software/) and phylogenetic trees were inferred using neighbor-joining method in FastTree program (http://www. microbesonline. org/fasttree/) and visualized with FigTree v1. 4. 2 (http://tree. bio. ed. ac. uk/software/figtree/), Archaeopteryx 0. 9920 (http:// en. bio-soft. net/tree/ATV. html), and iTOL, an interactive online program (https://itol. embl. de). A rhizome analysis  49  was carried out for MCR-1 to investigate putative ancestral recombination events of bacterial sequences leading to creation of chimerical gene.
For this purpose, protein sequence of MCR-1 (AKF16168. 1) was split into 10 fragments of 50 amino acids. Each fragment was blasted against the NCBI database using as threshold e-value 10e-5 to select the 100 best homologous sequences according to the aa identity ≥50% and coverage ≥80%. From these sequences, phylogenetic tree with each fragment was constructed to infer the closest ancestor. Looking for all MCR variants in bacterial genomes. Whole genome sequences of the most common Gram-negative bacteria were retrieved from RefSeq data: NCBI Reference Sequence Database (ftp://ftp. ncbi. nlm. nih. gov/genomes). This access to the repository is made through version 0. 2. 7 of the NCBI genome download scripts, which allowed the download of the entire genome assembly, including all required data for each genome, genomic sequences, annotated proteins and assembly structure reports.
Reference protein sequence of each MCR variant was also retrieved from the NCBI database. To the best of our knowledges, today (July 10, 2019), nine MCR variants have been described including MCR-1, MCR -2, MCR -3, MCR -4, MCR -5, MCR -6, MCR -7, MCR-8 and MCR -9 with respectively 18, 2, 29, 6, 3, 1, 1,2 and 1 subvariants. BlastP analysis was conducted to investigate the presence of the nine MCR variants (i. e. MCR-1. 1, MCR-2. 1, MCR-3. 1, MCR-4. 1, MCR-5. 1, MCR6. 1, MCR-7. 1, MCR-8. 1 and MCR-9. 1) in the large downloaded data. MCR hits were considered significant on the basis of 90% aa identity and 98% of alignment. Reference sequences of MCR variants were obtained from the Pathogen Detection database (https://www. ncbi. nlm. nih. gov/pathogens/). Analysis of genetic environment of mcr gene variants.
Genetic environments of mcr genes were investigated in the selected genomes using the Easyfig version x2. 2. 3 program (http://mjsull. github. io/Easyfig/) and CGView program (http://wishart. biology. ualberta. ca/cgview/). An annotation using Prokka (http://www. vicbioinformatics. com/software. prokka. shtml) and BlastP were made for genetic environment verification. Metadata and network analysis. Metadata of each bacterial genome containing significant hits of mcr gene variants (≥90% identity, ≥98% alignment) including bacterial strain name, bioproject number, collection date, geographical location, isolation source and host organisms, were collected from PATRIC database (https://www. patricbrc. org/) and detailed in Supplementary Table  2. To investigate putative links between bacterial strain harboring mcr gene variants, metadata informations are subjected to networking analysis using the cytoscape software v.