dataesr/scientific-paragraphs-categorization huggingface dataset:
Samples Selection and Data Quality All mammalian species with transcriptome annotations in the 99th release of the Ensembl Compara Database  [61] , for which RNA-Seq samples of healthy liver, kidney, lung, brain, or heart organs were available in the NCBI Sequence Read Archive, were selected in this study. The full list of species and sequencing run IDs can be found in Supplementary Table  S1. For each species, samples from juvenile, diseased, and very old animals were discarded in order to control for signals pertaining to the developmental or potentially pathological state of the included organisms. Only species that had at least two RNA-Seq samples were kept. In addition, to account for the heterogeneity in expression values, each dataset has been checked for anomalous distributions by exploratory analysis; extreme outliers were manually analyzed and removed.
Orthology Orthology relationships and transcriptome annotations have been obtained from Ensembl, release 99th, preprocessed, and imported to GraphDB, which was further used for the analysis of orthology and for expression value extraction. The preprocessing scripts have been developed in-house and are available as notebooks in our group's repository at https://github. com/antonkulaga/species-notebooks. For the linear models, expression values for genes in species with more than one paralog were considered as missing values to reduce ambiguity. The final dataset for linear models included 11,831 genes. For the LightGBM models, 12,323 orthologous coding genes were selected by combining one-to-one and high-confidence one-to-many orthologs. For this analysis, genes for which orthologs were not present in more than 90% of the species were excluded.
RNA-Seq Pipeline For quality control, adapter cutting and trimming, Fastp (version 0. 20. 1) [63]  was used. Transcript quantification was done with salmon (version 1. 4. 0) [64]. Transcript expressions were aggregated at the gene level with tximport (version 3. 12). Raw read counts were normalized on a per-sample basis, using the transcript per million (TPM) normalization. For performing the comparative analysis, the gene expression levels were normalized with TPM, which accounts for potential gene length differences across species.
Linear Models The linear models' analysis was performed in Python (version 3. 8), using several libraries, including statsmodels (version 0. 11. 1) [65] , pandas (version 1. 1. 1), seaborn (version 0. 10. 1), pyUpsetPlots (version 0. 4. 0), Venn (version 0. 1. 3), etc. The code necessary to reproduce the analysis can be found at https://github. com/ursueugen/cross-specieslinear-models. To investigate the association between gene expression and MLS, organ-specific linear models were built, allowing for the selection of genes highly associated with MLS, in each organ. Single-variable linear models of the form GeneExpression ~β0 + β 1 x SpeciesTrait were fitted independently for every gene and every life-history trait. For this analysis, the variables were log 2 -transformed and normalized to z-scores, prior to fitting.
For all results, please see Supplementary Table  S2. The expression activity across pathways was evaluated by using the single-sample gene set enrichment analysis (ssGSEA), also known as the signature projection method, as was first described by Barbie et al. [66]. For pathways data, we used the KEGG pathway database  [67]. The signature projection method was applied using the ssgsea function from the gseapy Python package (version 0. 9. 18)  [68] [69] [70]. After obtaining the enrichment score for each pathway, linear models of the form EnrichmentScore ~β0 + β 1 x SpeciesTrait were fitted and interpreted following the same approach as described for single genes. For the full pathway results, please see Supplementary Table  S3.
The hyper-parametric optimization process was performed independently from the gene selection process. The cross-validation configuration in the processes was similar: 5-fold stratified cross-validation was used in both, except that one fold was excluded in the hyper-parametric optimization (kept as a hold-out and used for evaluating the entire hyperparameter optimization process). Optuna sqlite databases with trials are provided with a source-code repository. To investigate the effects of gene expressions on MLS and other life-history traits, we used 6 separate regression LightGBM models for predicting from gene expression each of the species' features: maximum lifespan, body mass, gestation days, temperature, metabolic rate, and mtDNA GC%. For each of the models, the other five species variables were not used in the analysis, such that only genes would provide the prediction power for the target.
Bayesian Networks A Bayesian network encodes the conditional independence structure among a set of variables, in our case among the input genes and MLS. A potential causal relationship between a gene and MLS is identified when the gene belongs to the Market blanket of MLS in the Bayesian network. This is used as a basis for a variable selection approach that can detect causal relationships under strict assumptions. Implementation-wise, as part of the variable selection methodology, SES, a constraint-based variable selection algorithm  [43] , was used. This algorithm has proved to be appropriate in the context of high dimensional datasets  [43] , since the selection is not based on the optimization of an objective function (i. e. , loss function), so it is not prone to overfitting. In this study, the MXM R package was used.
e. , frequency of appearance of a gene across all identified gene signatures) was considered. Constructing the composite ranking included joining all the above-mentioned metrics (number of repeats, mean Kendall's tau, mean absolute SHAP value, maximum linear R 2 , and relative frequency), and additionally, GenAge mentions were also accounted for for each selected gene (boolean metric: 0 or 1, depending on whether the gene was reported in the GenAge database)  [3]. To form a list of genes considered most predictive by all models along with GenAge mentions, the composite rank was computed, as follows: each selected gene was assigned with 6 different ranks, whereby each rank indicated the rank of a gene in a value space of each of the 6 metrics mentioned above. Ranks were calculated with Pandas rank function with the method's parameter set to dense  [76].
An Explanatory Multilevel Linear Model for Composite Integration The linear model used for the explanatory analysis of the final selection of genes was a Bayesian organ-wise multilevel model with random coefficients. A normal likelihood, uninformative normal prior for each coefficient, and a gamma uninformative prior for the precision parameter, were used. The regression parameters were fitted using the Gibbs sampling algorithm implemented in the R Package RJags  [77]. The ORQ transformation method  [78]  was used to do a "transform-both-sides" regression with the imputed dataset since the variables (i. e. , gene expressions) and the target variable (i. e. , MLS) were nonnormal. The penalized deviance (2e4 iterations) for the multilevel models with the top 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 and 13 genes was 770.
Data Availability Statement: The preprocessing scripts are available at: https://github. com/ antonkulaga/species-notebooks; The bioinformatic WDL pipeline that was used for downloading, quality control and qualification is available at: https://github. com/antonkulaga/rna-seq/tree/ master/pipelines/quantification; Code for linear models analysis for genes and pathways is available at: https://github. com/ursueugen/cross-species-linear-models; Code for LightGBM+Shap analysis, the intersection of the results from multiple models and ranking is available at: https://github. com/ antonkulaga/yspecies; Code for Bayesian networks analysis and multilevel Bayesian linear modeling available at: https://github. com/rodguinea/bayesian_networks_and_bayesian_linear_modeling.
Materials and Methods Human Languages. We used the WCS database (www1. icsi. berkeley. edu/ wcs/). Two languages with extremely sparse information (judgments from 1 speaker only for at least some chips) were removed, resulting in 108 analyzed languages. English, which is not in WCS, was approximated based on the relevant figures from the study by Zaslavsky et al. (7). Agent Architecture and Training. Both agents are feed-forward NNs. Speaker contains 3 hidden layers, each of size 1,000 and with leaky-ReLU (rectified linear unit) activations. For each color, the Speaker's output layer defines a Categorical distribution over its vocabulary V. Listener is modeled as a linear layer of hidden size 5. The impact of agents' capacity on results is discussed in SI Appendix, Supplementary Information Text, 8.
The latter is an iterative method that alternates between evaluating P(w|c) and m(c) (Speaker's meaning for each c) until convergence. We refer readers to Zaslavsky et al. (SI Appendix, section 1. 4 in ref. 7)  for more details about this two-step process. When annealing β according to Agglomerative IB, the IB solution is initialized with the one found with the previous value of β. Optimization ends when β = 1. Data Availability. The models reported in this paper have been deposited in GitHub (https://github. com/rahmacha/EGG).
ACKNOWLEDGMENTS. We thank  Emmanuel Chemla ,  Thomas Brochhagen ,  Roger Levy ,  Louise McNally ,  Rachid Riad ,  Noga Zaslavsky , the PNAS reviewers, and, especially,  Gemma Boleda  and  Diane Bouchacourt  for feedback. Ted Gibson  and  Bevil Conway  generously shared their data with us. This research was funding by  Agence Nationale pour la Recherche Grants   ANR-17-EURE-0017   Frontcog ,  ANR-10-IDEX-0001-02 PSL *  , and  ANR-19-P3IA-0001 PRAIRIE 3IA.
Additionally, we analyzed 1472 unique plasma proteins measured by Proximity Extension Assay using the Olink platform (a dataset that we previously published)  26  , quantified cell-free DNA, and performed high-throughput antibody profiling (partial cohort recently published)  [27] [28] [29]  (Figure  1A , Method Details). After quality control analysis of the RNA-sequencing data  30  , we retained 698 samples from 370 unique patients for analysis (Figure  S1B-H ).
Since 100% neutrophil purity could not be guaranteed following enrichment for bulk RNA-seq, we next sought to determine the cell type composition of our samples using CIBERSORTx  38. We used the COVID-19 Bonn Cohort 2 fresh whole blood single-cell dataset as a single-cell reference for the deconvolution of our bulk data  11  (Supplementary Table  S1 ). We generated cell-type-specific signatures and used these to estimate fractions of mature neutrophils, immature neutrophils, monocytes, T and NK cells (together denoted T/NK), B cells, and plasmablasts for each sample (Figure  S2 ). We defined total neutrophils as the sum of the mature and immature neutrophil fractions. We found that lower estimated total neutrophil content was associated with lower ANC (Fisher's test, p = 1. 2x10  -17  ) obtained clinically, and lower estimated T/NK fraction was associated with lower absolute lymphocyte count (ALC, Fisher's test, p = 0.
Our score correlated strongly with a previously-defined NETosis gene signature (Figure  S14B )  51. The analysis revealed a significant enrichment of the NETosis metagene score on Days 3 and 7 in severe patients compared to non-severe patients (Figure  4A ). Additionally, across NMF clusters, we found higher NETosis scores in the immature clusters NMF1 (Pro-Neu) and NMF4 (Immature Activated) (Figure  4A ). Though the transcriptional analysis of NETosis revealed an association between an expression signature of NETosis and disease severity, many factors promoting NETosis, such as histone modification, are post-transcriptional  52  and would not be captured by RNA-seq. Therefore, we next searched for protein markers of NETosis in the matched plasma sample proteomic dataset  26. Several known protein markers of NETosis were found to be significantly associated with severe disease across all time points and were measured at different levels across neutrophil subtypes, including MPO, CXCL8, TNF, PADI4, HGF, and CD177 (Figure  4B , Figure  S14F ).
Plasma proteomics identifies neutrophil-driven secreted proteins and potential ligand-receptor interactions driving phenotypes To further understand the role of neutrophils in COVID-19 in relation to other blood and immune cells, we analyzed the plasma proteome using our existing Olink Proximity Extension Assay dataset for this cohort. We began by searching for protein markers which were significantly and uniquely associated with neutrophil NMF clusters (Figure  6A , Figure  S19 , Supplementary Table  S5 ). NMF5 (G-MDSC) in particular had strong upregulation of markers of disease severity and neutrophil activation such as S100A12, HGF, IL1RL1, IL1R2, DEFA1/1B, PADI4, and TGFB1 (Figure  6B ). Of note, TGF-β has been shown to influence B cells to class-switch to IgA when stimulated with LPS in vitro  65.
NMF4 (Immature) had the highest levels of ACE2 potentially indicating tissue damage, while NMF3 (PD-L1 + ISG + ) unsurprisingly showed enrichment for IFNL1, CXCL10, and IFNG, likely explaining the phenotype. Next, we sought to determine which severity-associated proteins in the plasma likely originated from neutrophils. To this end, we compared the log 2 (fold change) values between severe and non-severe patients on each separate time point on both the RNA and protein levels (Figure  6C -D, Figure  S20A-B ). For both protein and RNA analyses, we included age, sex, ethnicity, heart disease, diabetes, hypertension, hyperlipidemia, pulmonary disease, kidney disease, immuno-compromised status as covariates to control for confounding (and CIBERSORTx cell type fractions for the RNA). We identified several components of neutrophil granules (CD177, MMP8, MMP9, ARG1, S100A12, TGFA), factors involved in clotting (F3, SERPINE1), chemoattraction (CXCL8, IL4R), and inflammation/activation (FKBP5, FCAR, IL18R1, CLEC4D) upregulated in severe disease in both data types, suggesting that neutrophils are key contributors to the severity-associated plasma proteome.
Many studies have shown that NETosis is a defining feature of severe disease  49, 72, 73  , and here, we demonstrate that NETosis can be strongly induced by IgA antibodies, which may be the dominant effector function occurring in patients with high ratios of IgA1-to-IgG1 antibodies in plasma. Several potential therapeutics have been suggested for use in autoimmune disease contexts aimed at inhibiting NETosis such as PAD4 inhibition or administration of recombinant human thrombomodulin  74  , and similar strategies could be applied to NETosis in the context of severe COVID-19. Additionally, clinical trials targeting IL1B aimed at decreasing NETosis are underway (ClinicalTrials. gov identifier: NCT04594356). We also hypothesize that patients with high abundance of G-MDSCs would benefit the most from these therapeutics, as the NMF subtype ligand-receptor interaction analysis suggests that G-MDSCs may be induced by a positive feedback loop of IL1B signaling in circulating neutrophils.
First, we performed bulk transcriptomics rather than single-cell RNA-sequencing, so the neutrophil state gene signatures reflect a mixture of true neutrophil subtypes. Second, our samples were enriched for neutrophils via magnetic bead negative selection, and high purity of all samples across the cohort could not be guaranteed; thus, we used estimated cell type proportions as covariates in all analyses, but the expression of contaminating cell type-specific genes and cytokines cannot be excluded. Third, our time course data was collected on Days 0, 3, and 7 of hospitalization, though patients may have been infected for varying amounts of time prior to enrollment, so this cannot be considered an exact time course. Fourth, we only collected longitudinal samples from hospitalized patients, so we were unable to study patients pre-hospitalization or patients who were infected but never hospitalized.
Fifth, sample collection at later time points was biased towards sicker patients as they needed to stay in the hospital for a longer period of time. Sixth, we did not collect a validation cohort, so our findings will need to be validated in external cohorts with similar multimodal data structures. Seventh, our study relies on blood draws and thus only provides insights into circulating factors that play a role in COVID-19 disease pathology, yet the findings regarding type III interferons and the IgA-dominant humoral response strongly implicate mucosal immunity; future studies should focus on longitudinal SARS-CoV-2 immunity occurring along mucosal barriers. Finally, samples were collected at the onset of the pandemic before there were recommended treatment options, and thus it is not known how different drug treatments such as dexamethasone, tocilizumab, or other monoclonal antibody therapeutics would affect neutrophil phenotypes; furthermore, with the advent of COVID-19 vaccines, it is not known how the adenovirus or mRNA vaccines modulate neutrophil function by the antibody profiles that they produce.
Data and code availability Due to IRB consent limitations, raw sequencing data is not publicly available. However, the read count matrix and TPM matrix used in this study will be available in GEO and Supplementary Table  S1. All code used for the analysis is deposited in GitHub at https://github. com/lasalletj/COVID_Neutrophils and all additional files required to run the code will be deposited in Zenodo upon acceptance. Additionally we also collected blood samples from 8 healthy donors. For each patient, medical history and clinical data were collected and are presented in Table  S1  and as previously described  26. Samples were collected at three different time points: Day 0 upon admission to the ED (n=374 samples); Day 3 (n=212 samples) and Day 7 (n=143 samples) for COVID-19 + hospitalized patients.
MPO ELISA Supernatants from the NETosis assay were collected gently without disturbing cell adhesion and were stored at -80°C until the assay was performed. Samples were thawed at RT and spun down at 1500rpm for 5 minutes, and supernatants were transferred to a new plate, leaving a small amount of liquid at the bottom of each well. The ELISA for quantification of MPO was performed using the LEGEND MAX TM Human Myeloperoxidase ELISA Kit (BioLegend, Cat#440007) according to manufacturer specifications. The standard curve was fitted with a 4-parameter logistic curve-fitting algorithm using the dr4pl package in R.
Institute GTEx-TOPMed RNA-seq pipeline (https://github. com/broadinstitute/gtex-pipeline/), with an appended SARS-CoV2 genome. GENCODE v35 with the appended SARS-CoV2 GTF was used for annotation. Raw FASTQ files were aligned to the custom genome FASTA in the Terra platform with the Broad Institute GTEx pipeline using STAR v2. 5. 3a, and expression quantification based on a collapsed annotation was performed using RSEM v1. 3. 0.
Quality control RNA-SeQC 2  30  (https://github. com/getzlab/rnaseqc) was used to calculate quality control metrics for each sample. Samples were excluded if they did not meet the following criteria: 1) percentage of mitochondrial reads less than 20%, 2) greater than 10,000 genes detected with at least 5 unambiguous reads, 3) median exon CV less than 1, 4) exon CV MAD less than 0. 75, 5) exonic rate greater than 25%, 6) median 3' bias less than 90%. This filtration kept 698 out of 781 samples (89. 4%) (Figure  S1 ). Genes were included in the analysis if they were expressed at a level of 0. 1 TPM in at least 20% of samples and if there were at least 6 counts in 20% of samples. In total, 20283 genes passed the filtration criteria.
Neutrophil fraction estimation and contamination control CIBERSORTx  38  was used to estimate the proportions of mature neutrophils, immature neutrophils, T/NK cells, B cells, plasmablasts, and monocytes in each sample. To generate the signature matrix for deconvolution, we utilized the single-cell RNA-seq data of PBMCs and neutrophils from whole blood from Cohort 2 of the Schulte-Schrepping et al dataset  11. Using the designations provided in the public data, we created pseudobulks for each cell type per patient by summing the counts of a given cell type, and we excluded pseudobulked cell type samples from individual patients if the cell type had less than 5000 counts. To generate the CIBERSORTx signature matrix, we set limits of 50 to 100 marker genes per cell type, and filtered for only hematopoietic genes.
Gene set enrichment analysis We performed gene set enrichment analysis using the fgsea package in R using the following pathway sets from MSigDB Release v7. 2: H, C5 GO BP. We also performed a search of MSigDB using the keyword "neutrophil" and added the following pathways: In addition to these pathways, we added gene sets corresponding to various neutrophil states and signatures: genes up-or down-regulated more than threefold in blood neutrophils from ARDS patients  44  , single-cell neutrophil clusters in blood or lung tissue of patients with lung cancer  42  , single-cell neutrophil clusters from blood of patients with sepsis  43  , and single-cell neutrophil clusters from COVID-19 patients and healthy controls  11. For single-cell cluster markers, if there were more than 100 marker genes per cluster, gene sets were selected as the top 100 genes ranked by p-value for enrichment in a given cluster.
Neutrophil state network analysis Neutrophil state gene signatures were taken from the same GMT file used for GSEA analysis in Figure  2F. The network was built using the igraph package in R. Edges were drawn between nodes if the Jaccard index between the two gene signature lists was greater than 0. 05. Edge width was scaled according to the overlap coefficient between the gene sets, and nodes were scaled according to gene set size and colored according to the number of neighbors in the graph.
Logistic regression models to predict severe COVID-19 on Day 0 Logistic regression models were built using the glm package in R. In order to ensure the stability and interpretability of the coefficients in the model, we included only COVID-19-positive patients on Day 0 who were not immediately discharged from the ED (Acuity Max 1-4) and who had complete data for ANC, ALC, D-dimer, CRP, LDH, and BMI measured at Day 0. For patients with Acuity Max 1-4, 7 patients had missing clinical data, and these 7 missing patients were not biased towards a particular severity according to Fisher's exact test. All parameters used were broken into discrete quintiles unless insufficient samples belonged to one category, in which case factor levels were combined in order to minimize the standard error of the coefficient estimation.
We combined factor levels for age, LDH, and BMI, leaving 4 factor levels for age and LDH, and 5 factor levels for BMI (BMI was the only category scored from 0 to 5). Models were built according to three tiers of parameters. Model 1: clinical characteristics (age, sex, ethnicity, heart disease, diabetes, hypertension, hyperlipidemia, lung disease, kidney disease, immunocompromised status, BMI), Model 2: clinical characteristics plus clinical laboratory values (ANC, ALC, Creatinine, CRP, D-dimer, LDH), and Model 3: clinical characteristics plus clinical laboratory values plus neutrophil gene signature scores (NMF1, NMF2, NMF3, NMF4, NMF5, NMF6, ARDS Up -Juss, ARDS Down -Juss). ROC curves and AUC values were calculated using the pROC package in R. Significance of model improvement was determined using the likelihood ratio test using the lrtest package in R.
Feature selection for the best predictors on Day 0 of severity within 28 days among the variables used in Model 3 was performed using LASSO with the glmnet package in R with 100 repeats of 5-fold cross validation  26. Model tuning was performed using the caret package in R. We ranked features according to the number of cross-validation folds in which they were selected for the LASSO model (Figure  3D ).