dataesr/scientific-paragraphs-categorization huggingface dataset:
Nuclei filtering and quality control Nuclei filtering was performed manually on each of the four samples using Seurat v3. 1  (24). Ranked barcode plots of UMIs and gene counts were used to identify the lower 'elbow point' and cell barcodes with UMI count or gene counts below the elbow point removed as empty droplets. Further steps were used to identify additional empty droplets: the "SCTransform" function  (25)  was used to normalise the data, prior to centering, scaling, principal component analysis (PCA) and high-resolution clustering (using 30 PCs and resolution of 2). Wilcoxon rank sum differential gene expression tests were used to identify up-regulated genes in each cluster (see Methods Section 'Differential gene expression tests'). Cell clusters that both lacked distinguishing markers and had a low median UMI or gene counts (typically 2 median absolute deviations lower than the median across all nuclei) were removed as likely empty droplets or poor-quality nuclei.
Assignment of cellular identity Each cell cluster was assigned to one of the major liver cell lineages using a priori marker genes (Supplementary Table  2 ). Populations identified as erythrocytes and thrombocytes were not included in downstream analyses due to their likely origin from contaminating blood. Populations identified as hepatocytes, cholangiocytes, mesenchymal cells, endothelial cells, and immune cells were merged into five separate Seurat objects (using the Seurat "merge" function) for separate analyses of each cell lineage. Batch effects across samples were removed in each merged object with Harmony  (26). Gene annotations were taken from the Ensembl annotation for Atlantic salmon, and in cases where the gene was not assigned a name, the name of the nearest annotated putative orthologue in rainbow trout (Oncorhynchus mykiss), zebrafish or mouse (Mus musculus) was used, or in some cases informed by BLASTp searches against the NCBI non-redundant database.
Differential gene expression tests Throughout this study, differentially expressed genes were defined using the Seurat function "FindAllMarkers", applying the Wilcoxon rank sum test with default cut-offs (multi-test adjusted p-value < 0. 05, log2-fold change > 0. 25, with expression of the gene in at least 20% of nuclei in the cluster tested). For the global analyses of major cell types, differential gene expression tests were performed for each defined cluster in turn versus all other nuclei in the dataset as the background. For the hepatocyte-specific analysis we compared each hepatocyte cluster to all other hepatocyte nuclei. For the immune cell nuclei, the comparisons made were between each T cell subcluster and all other T cell nuclei, each myeloid sub-cluster and all other myeloid nuclei, and between each of the other immune cell clusters (i.
Hepatocytes comprised most of the profiled nuclei (88. 1%), followed by cholangiocytes (4. 3%), immune cells (3. 5%), mesenchymal cells (2. 6%) and endothelial cells (1. 5%) (Supplementary Figure  2 ; Supplementary Table  1 ). Clustering and identification of these major cell lineages was repeatable across individual samples, with each sample contributing a large proportion of the nuclei (Supplementary Figure  2 ). Differential gene expression analysis A B D C FIGURE 1 Major cell types in the Atlantic salmon liver defined by 47,982 nuclear transcriptomes. (A) UMAP highlighting five main liver cell-type clusters according to a priori defined marker genes (see Supplementary Table  revealed markers for each major liver cell type (Figures  1C, D ; Supplementary Table  3 ). Marker genes for candidate cell types and sub-populations are hereafter reported according primarily to annotations provided by Ensembl, or in cases where no Ensembl annotation was available, using supplementary BLAST homology searches against the NCBI database to support our inferences.
While Ensembl annotation utilizes phylogenetic information to inform homology relationships, it may nonetheless fail to correctly capture orthology of Atlantic salmon genes to mammalian species, particularly for fast evolving and complex gene families. Adding to this challenge, genetic orthology is not a prerequisite for conservation of gene function or expression, which limits our ability to transfer knowledge about cell marker genes from mammalian studies to Atlantic salmon. While the reader must be aware of these caveats, they represent a general issue in functional genomics studies using non-model taxa like salmonids. While it was possible to identify the five major liver cell lineages of Atlantic salmon using orthologues to marker genes defined in mammalian liver scRNA-Seq datasets (e.
g. 6), many markers were notably absent or expressed at very low levels in our snRNA-Seq dataset (Supplementary Figure  3 ). For example, Atlantic salmon orthologues of the widely used epithelial marker Epcam did not show expression in the cholangiocyte cluster (Supplementary Figure  3 ). Likewise, salmon orthologues of Pecam1 and Pdgfrb, which are excellent markers of mammalian endothelial and mesenchymal cells, were not detected at significant levels in these cell types in our dataset (Supplementary Figure  3 ). These results may be explained by differences in transcriptome composition between snRNA-Seq and scRNA-Seq datasets from liver, as identified in humans (30). In support of this idea, epcam and pecam1 were respective markers of cholangiocytes and endothelial cells in a recently published zebrafish scRNA-Seq dataset  (31) , making it less likely that the lack of expression of these markers in our Atlantic salmon dataset represents a true evolutionary difference between mammals and teleosts.
Significant markers for the latter were dominated by genes encoding APPs, many of which were observed to be expressed across all cell types (Figure  2B ). This is likely a consequence of the numerically dominant hepatocytes (Figure  1A ) leaking mRNA from highly expressed genes into the ambient RNA. Similar examples of this type of leakage can be observed in many liver scRNA-Seq datasets with, for example, hepatocyte expressed genes Alb and Hp being widely 'expressed' across non-hepatocyte cells in past mammalian studies (e. g. 3, 32). The immune compartment contained identifiable T and B cells, along with candidate populations of neutrophils, myeloid cells, and NK-like cells (Figure  2B ). The transcriptome of cholangiocyte nuclei was homogeneous, while limited heterogeneity was identified in the endothelial and mesenchymal cells (Figure  1A ; Supplementary Figure  5  and 6 ).
In our Atlantic salmon dataset, one predicted orthologue of Clec4f is expressed across all four myeloid populations (ENSSSAG00000076658), showing h i g h e s t e x p r e s s i o n i n M y 2 , w h i l s t t h e o t h e r (ENSSSAG00000040735) was expressed across all immune cell types (Supplementary Table  11 ). Furthermore, both genes were among the top up-regulated genes in hepatocytes following Aeromonas infection (Supplementary Table  8 ). These findings highlight the potential for major differences in marker gene cellspecific expression and presumably function for homologous genes shared by mammals and teleosts. Having robust marker genes for different cell-types is clearly essential for accurate studies of cell biology in any species. However, single cell transcriptomics studies performed to date, including this study, demonstrate that marker genes can vary markedly across species, especially for immune cells  (96, 104).
Funding This study was funded by grants from the  Scottish Universities Life Sciences Alliance (Technology Seed Funding Call) , the  University of Edinburgh's Data Driven Innovation Initiative (Scottish Funding Council Beacon 'Building Back Better' Call) , and the  Biotechnology and Biological Sciences Research Council , including the institutional strategic programme grants  BBS/E/D/10002071  and  BBS/E/D/20002174  and the responsive mode grant  BB/W005859/1. NH is supported by a  Wellcome Trust   Senior Research Fellowship in Clinical Science  (ref. 219542/Z/19/Z ).
| DuraClone B27™ T cells were selected on a CD3 PE-Dyomics649 (B695/30 nm)/SSC dot plot (Figure  1k ). Beads were used to standardize the detection settings at target values of 19 ± 0. 5 for B575 and 15 ± 0. 5 for B675 given by the manufacturer (Figure  1l ). A minimum of 3000 beads were acquired. HLA labeling was analyzed on a HLA-B27-PE (B575/30 nm) histogram (Figure  1m -n). Values in the "indeterminate zone", defined by the manufacturer between 15 and 23, were not conclusive. Flow cytometric data were analyzed using the Kaluza™ software (Beckman-Coulter) and statistics were performed with Microsoft™ Excel software. Results were expressed as median fluorescence intensity (MdFI) + / À 1 standard deviation and compared using Student's t test.
| Routine immunophenotyping of HLA-B27 using the IOTest™ kit Over a period of 2 years, 848 samples were routinely tested using the IOTest™ kit. One hundred ninety seven (23. 2%) samples were immunotyped as HLA-B27 and 592 (69. 8%) samples came out as HLA-B27 negative. The HLA-B27 status of 66 (7. 8%) samples could not be determined because of double labeling with anti-HLA-B27 and anti-HLA-B7. Comparatively from HLA typing National Institute data base, 7. 1% had HLA-B27 and 17. 1% had HLA-B7 genotypes in the population from the same geographical area. Other most frequent genotypes from the CREG groups were the HLA-B13 (4. 1%), HLA-B41 (2. 2%), HLA-B55 (3. 5%) while the remaining genotypes represented less than 1% of the population. In our routine analysis, 504 patients were females (59. 4%) and 344 were males (40.
COVID-19 presents a complex disease that needs to be addressed using systems medicine approaches that include genome-scale metabolic models (GEMs). Previous studies have used a single model extraction method (MEM) and/or a single transcriptomic dataset to reconstruct context-specific models, which proved to be insufficient for the broader biological contexts. We have applied four MEMs in combination with five COVID-19 datasets. Models produced by GIMME were separated by infection, while tINIT preserved the biological variability in the data and enabled the best prediction of the enrichment of metabolic subsystems. Vitamin D3 metabolism was predicted to be down-regulated in one dataset by GIMME, and in all by tINIT. Models generated by tINIT and GIMME predicted downregulation of retinol metabolism in different datasets, while downregulated cholesterol metabolism was predicted only by tINIT-generated models.
Predictions are in line with the observations in COVID-19 patients. Our data indicated that GIMME and tINIT models provided the most biologically relevant results and should have a larger emphasis in further analyses. Particularly tINIT models identified the metabolic pathways that are a part of the host response and are potential antiviral targets. The code and the results of the analyses are available to download from https://github. com/CompBioLj/COVID_GEMs_and_MEMs.
2 model was reported by Yaneske et al. [8]. GEMs can be used in a combination with different model extractions methods (MEMs), which allow us to automatically adapt a model and its constraints to a specific context using omics data  [9]. One of the main problems is that different MEMs produce variable, sometimes contradictory results  [10]. Different datasets gathered either from COVID-19 patients or cell lines have been applied to the analysis of COVID-19 metabolic signatures by genome-scale metabolic modelling and automatically generated context-specific GEMs  [11, 12]. However, only selected MEM approaches have been applied in the present studies. Namely, Nanda and Ghosh  [11]  applied iMAT and Cheng et al. [12]  applied tINIT extraction method to variable scope of COVID-19 datasets. In this study, we describe the integration of various publicly available COVID-19 datasets using different MEMs in combination with the Human-GEM metabolic model  [13].
We reconstruct an individual metabolic model for each of the datasets and healthy groups and each of the MEMs, namely iMAT  [14] , INIT  [15] , tINIT  [16]  and GIMME  [17]. We perform a thorough comparison of obtained models and discuss their biological relevance. Moreover, we generate a set of flux samples  [18]  using each of these models, which are used to identify enriched metabolic subsystems in each of the models between the infected and healthy groups. These are then used to perform a comparative analysis of the obtained results, which we discuss in a broader biological context.
However, present studies have only applied a single MEM to the analysed datasets. Since each MEM is based on certain assumptions different MEMs might produce significantly different or even contradictory results  [10]. The predictive power of extracted models can be increased if different MEMs are analysed in a given context to guide the selection of a method that yields models with the largest significance (e. g. , models preserving the separation between the observed experimental groups)  [39]. Furthermore, the selection of pre-processing and omics data integration steps also has a significant effect on obtained results and the use of various tools, approaches and datasets significantly increase the accuracy, precision and robustness of analyses  [40] [41] [42] [43]. Setting up workflows combining the integration of different datasets and tools enabling systematic assessment of omics data is one of the main goals of systems medicine  [44, 45].
In this work, we describe an approach to perform a selection of the most suitable MEMs in a combination with different COVID-19-specific datasets. The proposed selection is guided by different analytical approaches including principal component analysis (PCA), evaluation of model distances based on Jaccard index metrics, and analysis of the reactions and their flux values observed in different settings. Moreover, using the selected MEMs we perform detailed analyses of obtained context-specific GEMs to compare the metabolic states between healthy and infected cells. There are two main contributions of this study. Firstly, we propose a set of analyses that can be used to identify the most suitable MEM for an exact GEM-based analysis of observed datasets. Secondly, we analyse the models produced by selected MEMs in a wider biological context and identify the metabolic pathways that are a part of the host response and thus present potential antiviral targets.
COVID-19 datasets We collected the data needed to reconstruct the models from publicly available repositories. We identified the repositories with transcriptomes of healthy and COVID-19-specific human bronchial epithelial (HBE), lung biopsy cells, human embryonic kidney (293T), Calu-3, and adenocarcinoma human alveolar basal epithelial (A549) cell lines as described by Blanco-Melo et al. [46]  and Weingarten-Gabbay et al. [47]. Since each dataset belongs to a different cell type, it also defines the biological context in which the measurement of the transcriptome was performed. Next, we used the SRA toolkit  [48, 49]  to download the raw transcriptomic data, from each of the repositories. We used the Kallisto tool  [50]  to obtain the transcripts-per-million (TPM) values for the observed transcripts.
Extraction of COVID-19 specific models We used the Human-GEM model version 1. 6. 0 [13]  as a scaffold model for the extraction process. 1  Before reconstructing infected models, the model was augmented with the virus biomass objective function (VBOF) as described by Nanda and Ghosh  [11]. Different model extraction methods (MEMs) were used in the reconstruction, namely iMAT  [14] , INIT, tINIT  [16] , GIMME  [17]. Models were reconstructed for both the infected and the healthy group. 75th percentile of the average TPM values within a dataset and within a condition (infected or healthy) was used as a threshold value for each of the algorithms. During the reconstruction, the default Human-GEM model biomass function was used for the healthy state and VBOF for the infected state. After the reconstruction, we verified that a generated model could perform the essential metabolic tasks.
The reference model, which was used as a scaffold model in the extraction process, has more than 13,000 reactions. Extracted models eliminate inactive reactions for a given context but can still preserve up to more than 3,000 active reactions (see 2). Dimensionality reduction techniques such as principal component analysis (PCA) allow us to project high-dimensional data to a space with a manageable level of dimensions. This projection can be used to guide further analyses. For example, PCA plots in two-dimensional space can be used to assess how well obtained models cluster together in dependence on different factors. In our case, we aimed to find a configuration of the extraction process (namely, selection of the most suitable MEM) that would be able to separate the healthy models from the infected models.
Moreover, PCA can be used to assess the relative impact of each of the observed factors (namely MEM, dataset, and infection) on a model content  [10]. We thus opted to analyse the obtained models using PCA, the results of which were also confirmed by t-SNE (t-distributed stochastic neighbour embedding) plots  [52, 53]. The analysis was conducted in a similar manner as described by Opdam et al. [9]  and later applied by Walakira et al. [39]. We constructed a matrix describing the presence or absence of each metabolic reaction in each of the reconstructed models. We evaluated the amount of the variability between the models that is explained by each of the factors of extraction, namely MEM, dataset, and presence or absence of the infection. We observed the clustering of obtained models and their compliance with predefined groups using PCA and t-SNE plots.
82 (10-fold activation), and for the latter the cutoff value of at most -0. 82 (10-fold inhibition) between the infected and the healthy groups  [11]. GEMs are composed of different metabolic subsystems containing metabolic reactions with specific functionalities. For example, Human-GEM, which was used as a reference model in our analysis, is composed of 142 subsystems ranging from Glycolysis/Gluconeogenesis to Cholesterol Metabolism. We used the generated lists of significantly up-and down-regulated reactions for each of the MEMs and each of the datasets to perform the enrichment analysis of the metabolic subsystems with the hypergeometric test as described by Walakira et al. [39]. We identified the metabolic subsystems that were consistently enriched in the same dataset using different MEMs and/or consistently enriched in different datasets using the same MEMs.
These results were then used to perform the biological interpretation of the analysed datasets and compared with related COVID-19 studies. Reconstruction of context-specific models, their curation, and sampling of model fluxes were performed with the COBRA  [55]  and RAVEN Toolboxes  [56]  in Matlab R2019b (MathWorks Inc. , Natick, Massachusetts, USA) using the Gurobi solver (Gurobi Optimisation, LLC, Beaverton, Oregon, USA). Further analyses of reconstructed models, visualisation, and enrichment analysis were performed in Python 3. 2.
Models vary by size in dependence on different factors We analysed how the size of the extracted models depends on different factors, namely on a MEM, a dataset, and infection. The distributions of model sizes in dependence on different factors are presented in Fig. 2(A ). Differences in model sizes were significant only between different MEMs (p < 10 -5 , Kruskal-Wallis H-test). Models produced with GIMME pertained to the largest number of reactions (mean model size was 7134. 6 reactions) and models produced with iMAT the smallest number of reactions (mean model size was 3035. 4 reactions). To identify the number of reactions that carry nonzero flux in simulations, we generated a set of flux samples for each of the models. We evaluated the size for each of the models by counting the reactions that had nonzero flux in at least one of the samples.
Metabolic reactions should be analysed on the quantitative level We proceeded with the analysis of the obtained models based on the generated flux samples. We analysed the specificity of active reactions in dependence on each factor using the Jaccard index metric  [51]. We analysed the distribution of this metric within each factor and between the factors (see Fig. 3 ). We further analysed the significance of specificity of active reactions of each of the factors using the two-sample Kolmogorov-Smirnov test. We adjusted the calculated p-values for multiple testing using the Benjamini-Hochberg procedure. Active reactions were significantly specific for each of the MEMs and all the datasets except for the HBE data. Reactions were not significantly specific for neither healthy nor infected models.
2 ). We further analysed the models produced with each MEM to evaluate the amount of variability explained by the remaining factors. In all cases, a dataset explained a significantly larger amount of variability than infection (see Supplementary Table  3 ). Within the PC1 the former equalled almost 94% for iMAT (see Supplementary Figure  3 ), 82. 3% for GIMME (see Supplementary Figure  4 ), 54. 6% for INIT (see Supplementary Figure  5 ), and 71. 5% for tINIT (see Supplementary Fig. 6 ). However, even though a dataset explained a larger amount of variability, models extracted using GIMME could be separated by infection using the PC2 explaining 21. 4% of variability in the data, which was almost as much as PC1, which explained 25. 3% of variability (see Supplementary Fig. 4 ). Models were also partially separated by infection using tINIT and PC1, which explained around 23.
3% of variability in the data. The only outliers were healthy models obtained using lung biopsy cell data (see Supplementary Fig. 6 ). However, models produced by tINIT were separated also by the datasets, which represent separate biological entities. These were cells or tissue infected by the SARS-CoV virus, which differed in cell-specific metabolic steady-states, levels of viral load and host response to the virus. GIMME models did not separate by the dataset, while PCA analysis of the original data showed strong separation by the dataset/cell type except for lung biopsy cell data. Separation by infection and cell type in GIMME and tINIT models indicated that these two MEMs can generate models which reflect large compliance with predefined groups of data and should have a larger emphasis in further analysis.
Results of the enrichment analysis vary with a MEM We performed the enrichment analysis using the two-sample Kolmogorov-Smirnov test to identify up-and down-regulated metabolic reactions. Full results of the enrichment analysis are available at the supplementary GitHub repository. Furthermore, we used the lists of enriched metabolic reactions in a combination with metabolic subsystems as defined in the reference model  [13]  to perform the hypergeometric test, and to identify up-and down-regulated metabolic subsystems in each of the models (see also  Section 3). Different MEMs and different datasets produce different, in some cases contradictory results  [9]. We analysed a specific MEM by extracting the metabolic subsystems that were consistently enriched (either up-or down-regulated) in a predefined number of models obtained with all of the observed datasets and with a given MEM.