key: cord-1011339-n5bnpzj6
authors: nan
title: A blood atlas of COVID-19 defines hallmarks of disease severity and specificity
date: 2022-01-21
journal: Cell
DOI: 10.1016/j.cell.2022.01.012
sha: 3c53e3017a78ef35020ec5fb98b553487ce357a4
doc_id: 1011339
cord_uid: n5bnpzj6

Treatment of severe COVID-19 is currently limited by clinical heterogeneity and incomplete description of specific immune biomarkers. We present here a comprehensive multi-omic blood atlas for patients with varying COVID-19 severity in an integrated comparison with influenza and sepsis patients versus healthy volunteers. We identify immune signatures and correlates of host response. Hallmarks of disease severity involved cells, their inflammatory mediators and networks, including progenitor cells and specific myeloid and lymphocyte subsets, features of the immune repertoire, acute phase response, metabolism and coagulation. Persisting immune activation involving AP-1/p38MAPK was a specific feature of COVID-19. The plasma proteome enabled sub-phenotyping into patient clusters, predictive of severity and outcome. Systems-based integrative analyses including tensor and matrix decomposition of all modalities revealed feature groupings linked with severity and specificity compared to influenza and sepsis. Our approach and blood atlas will support future drug development, clinical trial design and personalized medicine approaches for COVID-19.

Correlation network: GM-CSF, CXCL10, TREM-1, CCL2/19, TF, IL-6/15, MPO, S100A9 (calprotectin) Machine learning: SAA2, CRP, IGHG4, CCL20, CCL2, IL-6 and C5a

Clustering by plasma protein biomarkers

Further genes, processes and pathways COVID clonotype

Plasmablast proliferation (COVID-19 specific) 

Neutrophilia CD274

CD64:10 sepsis>COVID-19 CD49d:CD43 in COVID-19 cMono (cycling/prolif; S100A8/9 hi ) ncMono pDC and cDC intMono in COVID-19 activated CD4 + ,CD8 + in COVID-19

To further deconvolute biological pathways and cellular functions associated with COVID-19 we analyzed gene expression at single cell resolution (STAR Methods; Data S6; Table S3 ). Among hospitalized and community COVID-19 cases (prioritized sample set), we found the principal components of variance in gene expression involved cMono, naïve B cells, plasmablasts, CD4 + T cells (naïve, effector and effector memory) and cycling NK cells (Figures 2H, I and S4C) . We performed pathway analysis for differentially expressed genes in major cell types ( Figure S4D ; STAR Methods) (Liberzon et al., 2015) . This showed enrichment for type I and II interferon pathways in the less severe hospitalized COVID-19 patients across cell types. Redox state (reflected by MTORC1 signaling and oxidative phosphorylation) was enriched across mononuclear phagocytes (MNP), T cells, NK cells and plasmablasts in more severe COVID-19, as were cell cycle (MYC targets, E2F targets, G2M checkpoint) pathways (except for MNP) while IL2-STAT5 pathway enrichment in T cells was found in more severe disease. Interferon stimulated genes showed enrichment in a range of cell subsets in COVID-19 cases, notably cDCs ( Figure S4E -G).

We then determined network modules for major cell subsets by WGCNA ( Figure 3A ,B; Star Methods; Data S6). Analysis of module co-variation identified five distinct module sets. A set of type I IFN response modules was found across cell populations in which representative gene expression profiles (module eigengenes) correlated with milder disease, better oxygenation status and earlier sampling from symptom onset ( Figure 3A ,B). The second module set, discovered in all cell types except plasmablasts, showed strong enrichment for activator protein 1 (AP-1) (FOS, JUN, ATF family genes) and the p38MAPK cascade. The module eigengenes were highly expressed in all COVID-19 patient groups, including recovery phase community cases, were distinct from influenza and sepsis, and did not show a consistent relationship with severity or other clinical features ( Figure 3A -D). The third module set was enriched for classical (C2H2) zinc finger (ZNF) genes and contained IRF2 and IL16; expression of these eigengenes was lower in COVID-19 and influenza compared with healthy volunteers and sepsis cases ( Figure 3A ,B). The fourth set of modules involved ribosomal proteins and inflammasome function (top genes by membership included NLRP1, MAP3K14 and FOXP1) and was negatively correlated with COVID-19 severity in monocytes ( Figure 3A ,B,E). Finally, we found a set of "cycling" modules which in cMono correlated with severity, and included S100A8/9 encoding calprotectin, a known severity biomarker (Silvin et al., 2020) ( Figure 3A ,B,F). We also identified two cell-type-specific modules associated with severe disease, a JAK-STAT/interleukin signaling module in CD4 + T cells, and an EGFR pathway-enriched module in cMono including the stress response gene FKBP5 and scavenger receptor CD163 (Figure 3A ,B,G,H).

We further investigated signatures of severity for specific mononuclear phagocyte populations ( Figures 4A,B ; STAR Methods; Data S3-6). Hospitalized COVID-19 patients with more severe disease had a relatively higher frequency of cMono and fewer intermediate monocytes (CD16 + CD14 + ), ncMono (CD14 -CD16 + ) and DCs . With increasing disease severity, we found a shift in the phenotype of cMono to lower expression of HLA-DR, CD33 and CD11c, and evidence of proliferating monocytes based on expression of Ki-67 and DNA abundance, with comparable changes in sepsis patients (Figures 4C and S5A, B) . Lower levels of pDCs and CD33 low cDC2 were found in sepsis compared with severe/critical COVID-19 ( Figure 4A ).

Further analysis of cMono using CITE-seq showed that a cycling cluster, and a cluster with high expression of the anti-oxidant metallothionein genes (MT hi ), were significantly elevated in critical COVID-19 cases, influenza and sepsis compared to healthy volunteers ( Figures 4B and S5D ). We found that S100A8/9/12 hi HMGB2-expressing cMono correlated with COVID-19 severity, and were also increased in sepsis. cMono expressing VCAN, which is implicated in cytokine release, were specifically increased in COVID-19 and reduced in influenza, whilst complement component C1Q-expressing ncMono were increased in influenza and sepsis, but not in . pDCs showed reduced abundance in more severe COVID-19, influenza and sepsis, as did CD1c + cDCs. Consistent with the progressive changes in abundance according to COVID-19 severity, the frequencies of cycling cMono and S100A8/9/12 hi HMGB2 expressing cMono, and CD1c + cDCs, were associated with clinical variables relating to oxygenation and respiratory function in hospitalized cases ( Figure 4D ). Pathway enrichment analysis of differentially expressed genes identified inflammatory response/TNF signaling and interferon response in milder disease including community cases (cMono, ncMono and DC), hypoxia (cMono, DC) and IL2_STAT5 pathways (cMono) across severity groups, and complement, coagulation and cholesterol metabolism in more severe disease in cMono ( Figure S5E -G).

To investigate epigenetic correlates of the COVID-19 response we analyzed chromatin accessibility by single-cell ATAC-seq (Figures 1A and S5H, I; Star Methods) . Overall, 750 and 303 accessible sites were up-and down-regulated respectively in COVID-19 patients compared to healthy volunteers in myeloid cells ( Figure S5J ). Genes linked to top differentially open chromatin peaks included STK24 (MAPK promoting apoptosis) and FGFRL1 (cell adhesion promoting fibroblast growth factor receptor) ( Figure S5J ,K). We identified the most significant DNA binding motif enrichments in the differentially accessible sites involved AP-1, SW1/SNF and BACH transcription factor family members which regulate chromatin remodeling and immunity ( Figure 4E ). Moreover, motif footprint analysis revealed increased accessibility of genomic regions containing FOS and JUN motifs in COVID-19 patients relative to healthy volunteers in myeloid cells, a signal which was also seen in convalescence ( Figure 4F ).

Page 8 of 92 reductions in the percentage of repertoire occupied by TRAV10, specific to invariant NK T cells, and TRAV1-2 and TRAJ33 usage, in keeping with reductions in MAIT cells ( Figure S7H ).

To better understand the relationship between COVID-19 and T cell clonality, we calculated Shannon diversity indices across clones based on the TCRbeta chain, controlling for age. While CD4 + subsets showed higher diversity than CD8 + subsets, differences with disease severity were only seen in CD8 + T cells (Figures S7I, J; STAR Methods) . Across disease states, accounting for age, CD8 + T effector memory (CD8.TEM/TEMRA), CD8 + T central memory (CD8.TCM/CD8.TCM.CCL5) and MAIT cell diversity was reduced in COVID-19 severe and critical disease with comparable changes in sepsis ( Figure 5G ).

Recent evidence suggests that effective CD8 + T cell responses involve increased numbers of expanded clones (Fairfax et al., 2020) . Consistent with this, we found hospitalized COVID-19 patients with mild disease had higher numbers of expanded clones in both CD4 + and CD8 + subsets and the mean clone size was higher within the CD8 + subset (Figures S7K, L) . In keeping with the observation that expanded CD8 T cell clones show increased expression of cytotoxicity markers (Watson et al., 2020b) , using a composite gene score for cytotoxicity we found that the number of expanded clones correlated with the average cytotoxicity score across all cells for that individual in both CD4 + T effector (CD4.TEFF/TEFF.prolif) and CD8 + T effector (CD8.TEFF/TEFF.prolif) populations ( Figure S7M ), and was higher in mild and community COVID-19 cases with reduced cytotoxicity observed in critical and severe disease ( Figure 5H ).

To further explore whether COVID-19 leads to generalized signatures of antigen presentation with reciprocal effects on TCR sequence and corresponding CDR3 usage, we devised an approach to identify COVID-19 associated amino acid sequences (Kmers of 4 amino acids) within the beta chain CDR3 region (STAR Methods). These were compared with chains from healthy volunteers and sepsis patients to exclude sequences non-specifically associated with infection (Data S3). We identified 125, 4-amino acid Kmers (referred to as COVSeqs) enriched in COVID-19 (Pc <0.05 versus both groups), the vast majority in CD8 + T cells ( Figure S7N ), with the proportion of cells with TCRs containing at least one COVSeq in the beta chain specifically increased in all COVID-19 patients ( Figure 5I ). In hospitalized patients we found a lower proportion of CD8 + T effector memory cells with COVSeq containing TCRs with increasing disease severity ( Figure 5J ). Critical disease was associated with naïve CD8 + T cells containing COVSeqs, indicating failure of the SARS-CoV-2 reactive cells in critical patients to expand into the effector phenotype, or possibly a distinct redistribution of the expanded cells. Further supporting functionality of the COVID-19 Kmer containing cells, the proportion of COVSeq containing cells was correlated with the median cytotoxicity of cells per individual amongst the COVID-19 patients ( Figure 5K, S7N) . Notably, COVSeq positive CD8 + T effector cells from critical patients showed reduced cytotoxicity compared to mild disease ( Figure S7O ).

Finally, we addressed whether using previously published COVID-19 associated beta chain clonotypes could further resolve variation in the T cell response according to disease severity (STAR Methods). We observed many cells carrying such TCRs across the COVID-19 patients, often overlapping COVSeq-containing cells. Notably, the distribution of these cells across clusters varied markedly according to COVID-19 disease state ( Figure 5L ). Replicating the observations with COVSeq-positive cells, CD8 + T effector memory cells were relatively depleted for COVID-19 clonotypes in critical disease ( Figures 5M and S7P ).

Page 9 of 92

We aimed to complement our multimodal cellular profiling with analysis of the COVID-19 plasma proteome. To do this we performed high-throughput liquid chromatography with tandem mass spectrometry (LC-MS-MS), presenting data for 105 proteins on 257 individuals (340 samples) ( Figure 1A ; STAR Methods; Data S3). We found differences by severity and etiology by analyzing principal components of variance ( Figure 6A ) and on unsupervised hierarchical clustering and supervised correlation analysis (Data S3). Severe disease, reflected in variance component loadings, was associated with increased acute-phase proteins and complement system proteins, including recognized biomarkers of inflammation (SAA1, SAA2 and CRP), complement membrane attack complex components (C5, C6, C9 and CFB), and functionally related protein families such as protease inhibitors (SERPINA3, SERPINA1 and ITIH3) and serum amyloid P-component (APCS) ( Figure 6B , Data S3). We also found differential protein abundance involving markers of tissue injury and necrosis, notably reduced extracellular actin scavenger plasma gelsolin (GSN); increased fibrinogens (FGA, FGB and FGG) ; and an increase in proteins implicated in IL-6 mediated inflammation (LGALS3BP, LRG1, LBP, HP and ITIH4). We further identified protein clusters based on the protein-protein interaction network, including a large cluster enriched for biological processes involving cholesterol transport and fibrin blood clots within which individual proteins showed positive and negative correlations with disease severity (PC1). Two smaller clusters enriched for cytolysis and complement activation positively correlated with disease severity (both showing negative correlations for all constituent proteins with PC1) ( Figure 6C ).

We found the main processes associated with differences between samples were acute-phase response and inflammation, metabolic (retinoid and lipoprotein) and cholesterol transport ( Figure  S8A ). Reduced levels of proteins associated with lipoprotein and cholesterol metabolism included apolipoproteins A-I, A-II, C-I and C-II (APOA1/2 and APOC1/2) and transthyretin (TTR), consistent with their downregulation in systemic inflammation and differences in metabolic state specifically associated with disease severity. This was further evident on pairwise comparisons, with mild hospitalized COVID-19 patients differing from healthy volunteers in metabolic processes and vesicle transport of retinoid, cholesterol, lipoproteins, and fat-soluble vitamins; and from community cases by higher levels of complement activation and coagulation ( Figure  S8B ). Severe COVID-19 patients differed from mild and from critically ill patients in processes relating to platelet degranulation and neutrophil degranulation respectively ( Figure S8C ,D). When we compared severe and critical COVID-19 with sepsis, 19 out of 105 proteins showed changes specific to FC>1.5) , enriched in acute-phase response, complement activation, and receptor-mediated endocytosis ( Figure S8D ).

To characterize inflammatory mediators of the response to SARS-CoV-2, we analyzed 51 circulating cytokine and chemokine proteins using the Luminex assay for 171 individuals ( Figure  1A ; STAR Methods; Data S3). There was clear clustering of hospitalized COVID-19 cases by severity on analysis of principal components of variance, while community cases overlapped with heathy controls and sepsis cases clustered separately ( Figure S8E ). The major proteins contributing to these axes of variance between groups were CXCL10, CXCL5, EGF, CCL2, S100A9, IL6, LCN2, CCL20, LF and G-CSF ( Figure S8E ). Overall, we found 49% (25 of 51) analytes were differentially abundant in plasma from COVID-19 cases versus healthy volunteers (Figures 6D and S8F ,G, Data S3). Amongst these, CCL2, CCL19, CCL20, CXCL10, GM-CSF, IL-6, IL-8, IL-15, S100A9 and SCGF (all increased abundance) were strongly correlated with severity in hospitalized COVID-19 patients (r 2 >0.5, P <0.001).

We further compared with sepsis and influenza to investigate disease specificity and found the plasma levels of G-CSF, IL-8, LF, CD163, LCN2, CCL20, IL-6, IL-10, CCL4, CCL19, TNF and C5a were lower in critical and severe COVID-19 than sepsis (Figures 6E and S8G) . Compared with influenza, serum EGF, LF and CD40L were higher in serum from patients with critical COVID-19 while G-CSF was lower ( Figure S8H ). We then investigated protein-protein correlation network relationships of assayed plasma cytokines and chemokines. This identified S100A9, M-CSF and CCL2/19 as nodal proteins. When we performed protein-clinical trait correlation network analysis for COVID-19 severity we found strong correlations (|r|> 0.5) between clinical features (CRP, SaO2/FiO2 and ventilation days) and specific nodal proteins (GM-CSF, CXCL10, TREM-1, CCL2/19, TF, IL-6/15, MPO and S100A9) at the center of the network ( Figure 6F ).

We next investigated the utility of plasma proteins for patient sub-phenotyping within hospitalized COVID-19 cases (n=122 samples) by integrating the LC-MS-MS and Luminex datasets using Similarity Network Fusion (SNF) (Wang et al., 2014) ( Figure 6G , STAR Methods). We first constructed a sample-by-sample similarity matrix from which we derived a network for each of the two data types. Analyzing these individually in an unsupervised manner with spectral clustering, we could only discriminate a minority of cases (the most mild from all others). However, when we fused these networks into a single similarity network that maximized shared and complementary information, we discovered two clusters that separated by clinical measures of disease severity including inspired oxygen concentration and SOFA oxygen score (t-test Pc <0.05) ( Figure 6G ). Notably, this molecular classification stratified severe cases assigned based on WHO categorical criteria into those that group with more mild cases and those clustering with critical cases. We identified 11 proteins as the main discriminatory features distinguishing the clusters (mutual information score >= 0.15) ( Figure 6G ). The predictive protein set spanned key inflammatory mediators, including the cytokines and chemokines IL-6, IL-8 (CXCL8), CCL2, CCL19, CCL20 and CXCL10 together with S100A9 (calprotectin), the acute phase proteins serum amyloid protein (SAA1) and protease inhibitor (SERPINA3), GM-CSF and the C-type lectin CLEC11A. When we compared the two clusters, we found that membership of cluster two was associated with higher 28-day mortality ( Figure 6H ).

We validated the clusters in an independent acute hospitalized COVID-19 cohort assayed using a different technology, targeted proteomics by Olink (Filbin et al., 2020 ) (STAR Methods). Clustering analysis, using 7 of the 11 predictive proteins for which data was available, identified two optimal clusters (Figures 6I). These showed a clear relationship with measures of disease severity, including WHO ordinal score (maximum) ( Figure 6I ), and patients in cluster 1 had lower mortality at 28 days (5/164=3.0%) compared with cluster 2 (33/105=31.4%) (Chi-squared test P <0.0001), validating the findings from our discovery cohort. We extended the approach to include a combination of hospitalized COVID-19 and sepsis patients from COMBAT. This revealed three clusters, two corresponding to the clusters seen with COVID-19 cases analyzed alone, indicating a high level of specificity ( Figure S8I ,J). Features that separated COVID-19 and sepsis included lipocalin2 (LCN2) and CCL20, which were elevated in sepsis, and CXCL10, APCS and fibronectin (FN1) which were higher in COVID-19.

We next used machine learning to combine the two proteomics data types with whole blood total RNA-seq to determine which features were predictive of disease severity and their relative informativeness ( Figure S9A ; STAR Methods). We first identified assay type specific informative components of variance to reduce dimensionality for a training sample set, and then determined which were most informative and the genes/proteins maximally contributing to each ( Figure 7A ).

After feature elimination based on performance, we found the minimal set of cross-modality features to predict severity were the acute phase proteins SAA2 and CRP, an immunoglobulin (IGHG4), chemokines (CCL20 and CCL2), IL-6 and complement component C5a; the combined performance of these features in the hold-out validation set showed a balanced accuracy of 75-80% to predict WHO category group (Figures 7B and S9B, C) . We also used machine learning to search for features that distinguish hospitalized COVID-19 patients from sepsis. A multi-omic set of 81 features was discovered using SIMON (Tomic et al., 2021) (STAR Methods) (AUC=0.85, 95% CI=0.59-1), identifying specific differentially abundant genes, proteins (including FCN1 and APCS as higher in COVID-19) and significant pathway enrichment for hematopoietic cell lineage and the renin-angiotensin system ( Figure S9D ,E).

We sought to dissect the immune response to COVID-19 across all assay types using a multiomics tensor approach Fanaee-T and Thoresen, 2018; Taguchi, 2017) , specifically the sparse decomposition of arrays (SDA) algorithm (Hore et al., 2016) . We analyzed 152 samples assayed for cellular composition, gene expression and plasma proteomics, and found 381 latent SDA components, each comprising vectors of scores (loadings) that indicate the contribution of individual cell types, genes or proteins linked by that component, and thereby offering insights into shared mechanism (Figure 7C,D; Table S4 ; STAR Methods). We identified components associated with specific clinical covariates, severity or patient group, noting that while in some instances e.g. gender there was a single associated component, typically several components were associated ( Figure S9F ). The strongest association with COVID-19 severity was for component 171 (Pc=5.9x10 -14 , rho=0.74 Spearman) ( Figure 7E ) which was unusual in having a high feature contribution from plasma proteins whereas gene expression contributed most to the majority of the other components ( Figure S9G ). Contributing features to component 171 included raised plasma chemokines involved in chemotaxis and activation (CXCL8, CXCL10 and CCL20) and GM-CSF together with acute phase activating proteins (SAA1/2 and SERPINA3), LRG1 and LBP; reduced abundance of intermediate monocytes; high expression of cell stress chaperone CLU and methyltransferase METTL7B, and downregulation of IgE receptor and multiple HLA class II genes; and pathway enrichment for antigen presentation, TCR signaling and asthma ( Figure 7E ).

To further delineate COVID-19 associated SDA components, we performed pairwise contrasts and analysis of variance involving COVID-19 patient groups. Overall, 130 of 381 components were associated with COVID-19 versus healthy volunteers ( Figure 7D ; Figure S9H ). Components associated with mild and severe but not critical disease included component 256 (upregulation of interferon response genes and downregulation of genes such as catalase and cytochrome c oxidase) which was specific to COVID-19 cases ( Figure 7F ); and component 42 (features of monocyte/granulocyte proliferation and function, elevated plasma proteins G-CSF, IL-2, IL-8 and IL-15, and enrichment of cell division related pathways) ( Figure S9I ). A further component strongly associated with severe disease involved plasmablast proliferation, combined with increased MT hi cMono and a clear DNA replication signature (component 289) ( Figure 7F ). We found an innate response component specific to critical COVID-19 (component 247) with differential expression of granulocyte activation marker (CEACAM8), neutrophil elastase (ELANE) and defensins (DEFA1B/4), and increased soluble CD163 scavenger protein levels, reflected in pathway enrichment for neutrophil functions ( Figure 7F ). Neutrophil related features were also found in component 123 associated with COVID-19 severity, influenza and sepsis ( Figure S9I ).

This approach also identified a further COVID-19 specific component (187) with high loading scores in hospitalized and community COVID-19 patients across NK, B and T cells ( Figure 7G ). This was driven by upregulation of key stress and activation response genes including immediate early response protein (PMAIP1), AP-1 transcription factor genes FOS and JUN, the early activation marker, tissue residency and metabolic reprogramming gene CD69, and TNFAIP3, which limits NFkB mediated inflammation. The cytokine-induced STAT inhibitor (CISH) and immune checkpoint regulator of inflammation and metabolism TNFAIP8L2 were downregulated. Pathway enrichment was seen for type-2 inflammation (IL4 and IL13), TLR signaling and the ATF-2 network. We additionally identified COVID-19 and influenza associated components ( Figure 7H ,I) including widespread upregulation of immunoglobulin heavy/kappa/lambda genes, JCHAIN (regulating multimerization and mucosal secretion of IgM/IgA), and MZB1 (involved in antibody secretion and integrin-mediated cell adhesion) linking with possible antibody-dependent cellular toxicity (component 6); and significant upregulation of interferon pathway genes (component 235).

Overall, the latent component analysis identifies hallmarks of COVID-19 severity, and specificity with respect to sepsis and influenza. Our findings highlight key cellular populations such as proliferating monocytes and plasmablasts, and features of innate and adaptive mechanisms ranging from interferon signaling to myelopoiesis, immunoglobulin production and stress activation response signaling. The results prioritize and validate hallmarks seen on individual modality analysis such as AP-1, and generate hypotheses for how hallmarks may be related in terms of pathophysiology based on co-occurrence in a given component.

Our comprehensive multimodal integrated approach, applied to multiple well-defined cohorts, has identified blood hallmarks of COVID-19 severity and specificity involving particular immune cell populations and their development, components of innate and adaptive immunity, and connectivity with the inflammatory response (Graphical abstract).

Identified hallmarks of severity involving myeloid related features include emergency myelopoiesis, immature neutrophils, increased HSC and platelet/CD34megakaryocyte progenitors, with the latter associated with thromboembolism. These findings substantiate and add granularity to previous reports (Bernardes et al., 2020; Stephenson et al., 2021) . We find evidence ERG is central to a gene network linked to these populations, encoding a transcription factor important in determining lineage plasticity, modulating inflammation and maintaining an anti-thrombotic environment (Yuan et al., 2009) . We further identify hallmarks supporting the importance of mononuclear phagocyte dysfunction in severe disease (Bost et al., 2021; Schulte-Schrepping et al., 2020) , namely proliferating cMono, and specific monocyte populations showing reduced HLA-DR, CD33 and CD11c expression, high expression of antioxidant metallothionein and S100A8/9/12 (calprotectin), together with reduced pDCs.

Changes in the frequency of specific T cell subsets, and their activation and exhaustion, has been previously reported in severe COVID-19 (Chen and Wherry, 2020; Jouan et al., 2020; Parrot et al., 2020) . We found increased numbers of activated CD8 + T cells and NK cell populations in COVID-19, and, with increasing severity, failure of clonal expansion in CD8 + T effector and central memory cells and depletion of COVID-19 clonotypes. We further found association with severity for exhaustion markers and specific activated NK and CD69 + MAIT cell populations. In terms of adaptive immunity (Brouwer et al., 2020; Galson et al., 2020) , we find increased numbers of terminally differentiated plasmablasts, with expansion of unmutated B cells differing in selection and tolerance, and a higher proportion of clonally related B cells.

Persisting changes in cell composition were found in convalescent samples, including an increase in activated CD4 + and CD8 + T cells and unswitched memory B cells as well as reduced CD16 -CD56 dim NK cells. These differences were also seen in community cases when compared to healthy volunteers; upregulation of TNF signaling was a consistent feature across cell populations in community cases. Redox state and cell cycle associate with more severe disease across cell populations. Our data are also consistent with the importance of the hyperinflammatory state (Moore and June, 2020) and interferon response (Hadjadj et al., 2020; Lei et al., 2020) but as features of less critical disease and earlier phase of illness.

Our proteomic analysis identified plasma levels of specific cytokines and chemokines as biomarkers of severe disease with evidence for acute phase response, complement activation/attack, fibrin clots, proteases, serum amyloid, tissue necrosis, receptor mediated endocytosis and cholesterol transport as hallmarks. Moreover, we discovered plasma protein signatures that can be used to stratify acute hospitalized COVID-19 cases into disease subphenotypes not identified from WHO categorical criteria, with cluster membership informative for response state and associated with differential 28-day mortality. We validated this in an independent dataset using a predictive set of seven plasma proteins (cytokines IL-6, IL-8; chemokines CCL2, CCL19, CCL20, CXCL10; and C-type lectin CLEC11A, a key growth factor for primitive hematopoietic progenitor cells). Patient stratification is important given the observed clinical heterogeneity within severe COVID-19. Such variability has historically confounded clinical trials for targeted immune therapy in other severe infections (Davies et al., 2018; Marshall, 2014) .

This work has demonstrated the informativeness of a multi-linear tensor approach to stratify and interpret the complexity of multi-omic datasets. Application of the SDA algorithm identified latent components of variance that linked signals from across cellular, gene expression and plasma protein measurements. For example, the most significant tensor component associated with COVID-19 severity involves myeloid chemotaxis and activation, the acute phase response, HLA class II downregulation and TCR signaling. Our dataset provides an important resource to develop other approaches for identifying multimodal signals and associated mechanistic insights, such as leveraging algebraic systems biology (Gross et al., 2016) , multi-layer networks (Kivela et al., 2014) , topological data analysis (Camara, 2017) or tensor clustering (Seigal et al., 2019) .

Considering features specific to COVID-19, we found the CD49d:CD43 ratio is elevated across COVID-19 patients, including into recovery and in convalescence, suggesting this may be an informative index score for neutrophil activities specific to COVID-19 that could be prognostically useful acutely and in convalescence. A subset of tissue-infiltrating CD49d expressing neutrophils have been reported, activated by vascular endothelial growth factor A released from hypoxic sites (Massena et al., 2015) , while reduction in CD43 expression leads to neutrophil retention in the blood and increased adherence to vessel walls (Woodman et al., 1998) . This raises the possibility of a link to the enhanced thrombosis observed in COVID-19 patients, the basis of which remains unresolved (Ortega-Paz et al., 2021) .

Our epigenetic, gene expression and integrative SDA analyses all identified upregulation of the AP-1 p38 MAPK pathway as a specific feature of COVID-19 across different immune subsets.

While not associated with severity, this observation further supports systemic immune activation and proliferation as a hallmark of COVID-19. AP-1 is a pioneer transcription factor that reversibly imprints the senescence enhancer landscape following stress and can be modulated to reverse T cell exhaustion (Lynn et al., 2019; Martinez-Zamudio et al., 2020) . This raises the possibility of AP-1 acting as a mediator of inappropriate chromatin remodeling and cellular activation/senescence in multiple cell types, with evidence from our analysis of persisting differential chromatin accessibility, gene expression signatures and cell protein markers of activation that may contribute to both acute disease and post-COVID-19 syndrome. Our findings Page 14 of 92 are consistent with work involving convalescent COVID-19 patients showing evidence of trained immunity and persisting chromatin remodeling in myeloid and lymphoid cell populations involving high levels of activity of AP-1 transcription factors (You et al., 2021) . Persistent reduced chromatin accessibility in myeloid cells at AP-1 targeted loci has been reported following influenza vaccination associated with innate immune refractoriness (Wimmers et al., 2021) . A role for AP-1 in COVID-19 is also consistent with the reported efficacy of baricitinib improving recovery of hospitalized patients (Kalil et al., 2021) . Baricitinib acts upstream of AP-1 (Zarrin et al., 2021) and controls macrophage inflammation and neutrophil recruitment in COVID-19 (Hoang et al., 2021) .

Our findings involving immune activation and proliferation are of further relevance given the key nodal plasma cytokines, including GM-CSF and IL-6, we identify as hallmarks of disease severity in proteomic and integrative analysis, targeting of which could also benefit dysfunctional granulopoiesis and neutrophil subsets. This supports the efficacy of tocilizumab (IL-6 receptor inhibitor) in severe disease (Horby et al., 2021b) . Moreover, we find upregulation of IL6R in more severe disease. Our data supports inhibitors of GM-CSF in current clinical trials, use of anti-CCL2 (Gracia-Hernandez et al., 2020) , the repurposing of immune checkpoint inhibitors (Pezeshki and Rezaei, 2021) to restore/enhance CD8 + T cell cytotoxicity, and potential targets such as FGFRs, cofactors for early-stage viral infection upregulated by MERS-CoV2 leading to lung damage (Hondermarck et al., 2020) .

Establishing optimal biomarkers of response and relevant sample collection, as well as timely availability of results, are important considerations in current immunomodulatory trial design in COVID-19. We found the widely used clinical biomarker CRP has predictive utility but less than serum amyloid A and is one component of a discriminating biomarker marker set for severity. In terms of multi-omics, total RNA-seq of whole blood stabilized at the bedside, a tractable sample type for collection in a pandemic, was highly informative for understanding variance in disease severity while use of small volumes of rapidly fixed whole blood for later flow and mass cytometry gives complementary and highly granular resolution of the cellular immune response state and identifies specific cell populations important in severe disease.

An important question is how the hallmarks of severe disease in COVID-19 relate to non-SARS-CoV-2 sepsis (Beltran-Garcia et al., 2020; Olwal et al., 2021) . We compared patients with severe or critical illness, not only revealing shared features relating to emergency myelopoiesis and progenitors but also identifying discriminating neutrophil markers (CD49d:CD43 ratio) in COVID-19 (summarized in Graphical abstract). In terms of monocytes, upregulation of cMono and specifically higher abundance of CD14 + CD33 + HLADRcells and S100A8/9/12 expressing cMono in severe/critical COVID-19 and sepsis is consistent with similar changes in monocytic myeloidderived suppressor cells in both diseases, while similar reductions in ncMono and DCs were seen in both COVID-19 and sepsis. We found specific differences for other populations, for example increased complement expressing (C1Q) ncMono in sepsis and extracellular matrix protein versican expressing cMono (mito hi VCAN hi ) in COVID-19. Lymphocyte exhaustion is proposed as a common mechanism in COVID-19 and sepsis (Boomer et al., 2012; Diao et al., 2020) . We find overall higher levels of CD4 + and CD8 + T cell activation, and more marked changes in cell markers, in COVID-19, suggesting a greater degree of CD8 + T cell exhaustion. While we find no major changes in naïve and transitional B cells in either condition, in sepsis we observed an absence of terminally differentiated plasmablasts and fewer clonally related B cells, class switching and autoreactive BCRs. We further identified plasma proteins discriminating COVID-19 and sepsis, for example LBP, lactoferrin and lipocalin 2; and differential pathways including neutrophil degranulation, complement, AP-1/p38MAPK, TLR, renin-angiotensin and the HSC lineage. We note that COVID-19 induced downregulation of ACE2 may be involved in modulating both the renin-angiotensin system and activation of AP-1/p38MAPK signaling (Grimes and Grimes, 2020) .

Similarly, while differences in target cells and control of viral replication are recognized between COVID-19 and influenza, shared immune responses and mechanisms are reported (Flerlage et al., 2021) . Our data builds on previous blood immune signatures (Lee et al., 2020; Mudd et al., 2020; Zhu et al., 2020) to gain cellular and proteomic insights into more severe forms of these two infections. We find similar changes in frequency of major immune cell populations between COVID-19 and influenza for HSC, the majority of T and NK cell populations, classical/cycling monocytes and plasmablasts; but differences for many B cell populations (differential naïve, intermediate and memory B cell response) as well as specific cell subsets notably ncMono subpopulations (summarized in Graphical abstract). While the bulk of the circulating proteomic response is shared, specific differences involved EGF, G-CSF and IL-15. We further found extensive sharing of pathways and networks including modules involving JAK-STAT and zinc finger proteins, as well as differences, notably enrichment of AP-1/MAPK in COVID-19.

In conclusion, our multi-omic integrated blood atlas comprehensively delineates the host immune response in COVID-19 from the start of the UK pandemic, prior to clinical trial-led implementation of approved treatments or vaccination. This provides the community with a unique reference resource for replication and meta-analysis, to interpret datasets generated from interventional trials and includes a browser for direct visualization of data (https://mlv.combat.ox.ac.uk/). Integrative approaches such as we have applied here are essential to better differentiate patients according to disease severity, underlying pathophysiology and infectious etiology. This will be important as we seek novel therapeutic targets and the opportunity for a precision medicine approach to treatment that is appropriately timed and targeted to those patients most likely to benefit from a particular intervention.

This study is limited to analysis of peripheral blood. Future multi-tissue studies including lung samples using multi-omics approaches are needed to further delineate the host response in COVID-19. It will also be important to more clearly establish changes from baseline, such as will be afforded by challenge studies with SARS-CoV-2, and more detailed temporal profiling over the course of acute illness into convalescence. Further work is needed to delineate the role of AP-1/p38 MAPK pathway upregulation in COVID-19, for example the relationship with viral replication and specific disease mechanisms such as T cell exhaustion. AP-1 inhibitors may enable insights into reversibility, the extent and nature of persisting epigenetic differences in convalescence, and relationship with COVID-19 vaccination. Functional studies are needed to understand the relationship between COVID-19 associated Kmers, CD8+ T cell expansion and cytotoxicity in critical disease but TCR sequences, including an agnostic Kmer analysis, have been used to define T cell responses linked to SARS-CoV-2 (Shoukat et al., 2021 critical; CComm: COVID-19 community case in the recovery phase (never admitted to hospital); CConv: COVID-19 convalescence (survivors from 28 days after discharge); Flu: influenza inpatient critical; Sepsis: in-patient severe and critical sepsis; SeConv: sepsis convalescence. 

Further information and requests for resources and reagents should be directed to and will be fulfilled by Julian Knight (julian.knight@well.ox.ac.uk).

This study did not generate new unique reagents.

Derived and processed data for all the datasets generated during this study and reported in this paper are available including through this paper, the European Genome-phenome Archive (EGA), Zenodo and COVID-19 Cell Atlas (as detailed in Key Resources Table) . For sequence level RAW datasets deposited at EGA, access is managed by the COMBAT Consortium Data Access Committee. Web-based interfaces for visualizing COMBAT datasets and outputs reported here are available at https://mlv.combat.ox.ac.uk/ and https://shiny.combat.ox.ac.uk. All code used for every algorithm followed in data processing and analysis is fully referenced within the specific methods text sections and Key Resources Table. EXPERIMENTAL MODEL AND SUBJECT DETAILS

The study was designed to allow deep molecular, multi-omic and immunological profiling of COVID-19 in peripheral blood at the outset of the pandemic in the United Kingdom during a public health emergency. Samples used were derived from multiple sources to allow comparison between adult patients with varying severities of COVID-19 and comparator disease or health states. A summary of the cohorts included in this study is provided below followed by a description of the methods used to define clinical parameters and handling of the blood samples. Demographic information and clinical phenotyping for all study subjects is summarized in Table  S1 , together with numbers of patients/samples assayed by cohort.

Patients admitted to Oxford University Hospitals NHS Foundation Trust, UK were recruited into the Sepsis Immunomics study (a prospective observational cohort study applying an integrated immune -omic approach to understand why some patients have a severe response to infection) if they were found to have a syndrome consistent with COVID-19 and a positive test for SARS-CoV-2 using reverse transcriptase polymerase chain reaction (RT-PCR) from an upper respiratory tract (nose/throat) swab tested in accredited laboratories, with the majority of patients co-recruited using the ISARIC/WHO Clinical Characterisation Protocol for Severe Emerging Infections. Written informed consent was obtained from adults or personal/nominated consultees for patients lacking capacity, with retrospective consent obtained from the patient once capacity was regained. Ethical approval was given by the South March and 28 April 2020. A selection of survivors were approached and asked to provide samples with consent under both research protocols at day 28 or more following discharge from hospital.

In order to provide a time of symptom matched set of samples from individuals with mild COVID-19 disease in the community, healthcare workers based at Oxford University Hospitals NHS Foundation Trust, UK with symptoms consistent with mild COVID-19 and a positive test for SARS-CoV-2 using RT-PCR from an upper respiratory tract (nose/throat) swab tested in an accredited lab were recruited into the Gastro-intestinal illness in Oxford: COVID substudy [Sheffield REC, reference: 16/YH/0247]. Individuals were consented and sampled at or after 7 days from the start of symptoms when the participants were returning to work.

In order to include samples from patients with sources of severe sepsis other than COVID-19, patients older than 18 years of age admitted to Oxford University Hospitals NHS Foundation Trust, UK with symptoms and signs of established sepsis (suspected infection with an acute change in total Sequential Organ Failure Assessment (SOFA) score of ≥2 points or a change in quick SOFA score by ≥2 points) were consented into the Sepsis Immunomics project [Oxford REC C, reference:19/SC/0296]) between 12 October 2019 and 13 March 2020. Patients were sampled on days 1, 3 and 5 of either hospital or intensive care admission and a selection of survivors were sampled, with additional consent at up to 6 months post-discharge.

Following advertisement, interested individuals 55 years or over and self-reporting as healthy were consented and recruited into the Genetic diversity and gene expression in white blood cells study [South Central Oxford REC B, reference 06/Q1605/55]. Blood samples were collected on one occasion only.

In order to provide comparator samples from patients critically unwell with influenza, samples were included from patients ≥18 years old managed in an intensive care unit (ICU) for ≥24 hours requiring ventilator support with either a diagnosis of COVID-19 or severe influenza using the Aspergillosis in patients with severe influenza ( 

Clinical data capture Healthy volunteers and healthcare workers with COVID-19 had age, sex and, where possible, selfreported ethnic background information collected. Information on previous medical history was not collected or stored for the healthy volunteer cohort but participants were invited on the basis of their self-reporting being 'healthy' and were deemed capable of self-presenting for study assessment. Healthcare workers with COVID-19 were asked for the date of onset of symptoms ('when they started to feel unwell') and the number of days between this date and sampling was calculated for every participant. Any healthcare worker who was admitted to hospital or required J o u r n a l P r e -p r o o f

Page 27 of 92 oxygen had this information collected and retained and defined as maximal severity of illness using the World Health Organization Criteria (see below).

A range of clinical data was collected and stored for downstream analysis using structured methodology from the hospitalized patients included in the study (COVID-19, all-cause sepsis and influenza).

Sex and ethnicity were captured using electronic healthcare records (EHR). Age was calculated at the time of sampling or maximal severity of illness using dates of birth registered in EHR.

Smoking status was derived from clinical clerking or direct patient or next-of-kin questioning where possible. Estimated or calculated weight was available for the vast majority of patients. In London COVID-19 and influenza patients, body-mass index (BMI) was calculated through availability of estimated or measured height. In Oxford patients, height was rarely available and so BMI was estimated using sex and weight and correlation with measured BMI was acceptable in cases where BMI was available. BMI of 25 or over was defined as overweight. Index of multiple deprivation (IMD) quintile and Charlson comorbidity indices (2012 definitions) were derived from the EHR in Oxford and were not available for London patients. Previous medical problems in patients admitted to Oxford hospitals were defined using NHS UK Read Codes Clinical Terms Version 3 and converted into ICD-10 codes using 'TRUD' (https://isd.digital.nhs.uk/trud3/user/guest/group/0/pack/9). Pre-existing medical conditions defined through ICD-10 were defined through the consensus of at least one clinical opinion and in cases where there was uncertainty if a code was assigned to an acute presentation, pre-existing conditions were only assigned if they were present in the first episode of an admission, and not if they were only present in later admission episodes. For patients recruited in London, pre-existing diseases were defined by the study teams using clinical records and patient or relative questioning.

After defining pre-existing ICD-10 codes or reported diagnoses, patients were classified as having the following conditions: hypertension, chronic respiratory disease (not explicitly available for London patients), asthma (not explicitly available for London patients but chronic obstructive pulmonary disease was), chronic cardiovascular disease, diabetes, hematological malignancy, other malignancy, liver disease, neurological disease (not available for London patients), chronic kidney disease, solid organ transplant, rheumatological condition (not available for London patients), significant immunosuppression (not available for London patients) and stroke / dementia (not available for London patients). An OpenSafely (OS) score was calculated for all patients using these classifications combined with BMI and IMD (Williamson et al., 2020 ). An adjusted OS score was calculated for the London cases that did not include some missing pre-morbid conditions and IMD quintile and the correlation was very high with the standard OS score (r=0.98).

Length of hospital stay was defined using hospital records for all hospitalized patients. Length of intensive care admission and duration of vasopressor, ventilation or, where appropriate extracorporeal membrane oxygenation (ECMO) and/or renal replacement therapy were all defined for hospitalized patients using EHR. All intervention and maximal severity timepoints were defined according to the date of onset of symptoms for each patient. This was defined by at least 2 independent clinicians through review of the clinical notes or direct questioning of the patient according to any unusual symptoms related to the current clinical condition. COVID-19 was defined by presence of at least one symptom consistent with COVID-19 and a positive microbiological test. Patients without symptoms were not approached for recruitment. Time

Page 28 of 92 between symptom onset and sampling was measured in days. All London patients were classified as critical requiring mechanical ventilation. No specific physiological observations were available for the London patients at the time of sampling and were only available at the time of admission onto ICU for calculation of the Acute Physiology and Chronic Health Evaluation (APACHE) and SOFA scores. Physiological observations were available from EHR for all hospitalized Oxford patients and were defined on the day of sampling at midday or the closest time before or after midday for each patient and timepoint including oxygen saturation, delivered fraction of inspired oxygen (either exact or estimated depending on delivery method: 0.24 for nasal cannula, 0.28 for simple face mask, 0.8 for non-rebreathe mask), pulse rate and blood pressure. Oxygen saturation: fraction of inspired oxygen ratio (SaO2:FiO2) was calculated as a proxy for partial pressure of oxygen:FiO2 ratio in the absence of complete arterial blood sampling for all patients. Days of oxygen therapy were defined in relation to sampling and maximal illness and total hospital stay for hospitalized COVID-19 patients where 2 or more timepoints of oxygen therapy were recorded. Fever was defined as 38.0ºC and, if present at midday, was defined at present at time of sampling. Highest measured temperature 24 hours prior to midday was also captured using EHR for every timepoint and the presence of 'persistent fever' more measures of temperature ≥38ºC in 24 hours before midday than temperatures recorded <38ºC was also recorded for all sampling timepoints in hospitalized Oxford patients. In those patients where persistent fever was observed the time in days between onset of symptoms and end of persistent fever was calculated.

Maximal severity was defined for Oxford COVID-19 patients using SaO2:FiO2 ratio and the date of lowest value was defined through manual inspection of longitudinal physiological parameters available through EHR. The individual levels of SaO2, method and exact or estimated concentration of oxygen delivered and additional recruitment strategies such as paralysis or proning were captured. In cases where death occurred, this was defined as maximum severity of illness and time between symptom onset and maximum severity was defined thereafter.

All together this information was used to understand the illness response and severity including (1) to generate a correlation matrix of clinical features, severity scores and markers illness response in the hospitalized COVID-19 cohort (Data S2); (2) to perform unsupervised machine learning (see section Statistical and unsupervised analysis of clinical phenotyping using clustering; also Data S2); (3) to define WHO categorical (mild as no oxygen, severe as ≤93% oxygen saturation, and critical as requiring mechanical ventilation) and ordinal scales (https://www.who.int/blueprint/priority-diseases/key-action/COVID-19_Treatment_Trial_Design_Master_Protocol_synopsis_Final_ 18022020 .pdf) for each timepoint of sampling and maximal severity of illness and these times were all aligned in days from onset of symptoms (summarized in Data S2). Severity was also classified according to the WHO Ordinal scale (https://www.who.int/blueprint/priority-diseases/key-action/COVID-19_Treatment_Trial_Design_Master_Protocol_synopsis_Final_18022020.pdf) and Sequential Organ Failure Assessment score (Ferreira et al., 2001) together with the extent of lung disease and compromise (SOFA oxygenation score).

Full blood count differentials (hemoglobin concentration, platelet count, total white cell, neutrophil, lymphocyte, monocyte and eosinophil counts), highest lactate, C-reactive protein (CRP), alkaline phosphatase (ALP), alanine transaminase (ALT) and lowest bicarbonate were captured at the time of sampling for all hospitalized COVID, and sepsis patients using EHR (Data S2). Equivalent results at the time of sampling were not available for London patients. D-dimer, lactate dehydrogenase (LDH) creatine kinase (CK) levels at the time of admission, sampling and maximal severity of illness were available for hospitalized Oxford patients (Data S2) but not London patients. Decisions on maximal levels of care based on senior lead clinician decision (frequently J o u r n a l P r e -p r o o f Page 29 of 92 with consensus from senior intensive care clinicians) were captured and patients were defined as frail if the plan was not to offer multi-organ support.

Patients with computerized tomography (CT) images of the lungs and thorax had images assessed independently by three radiology or respiratory experts to define ground-glass or consolidation patterns of lung involvement, presence or absence of pulmonary embolus and extent of involvement.

The clinical suspicion and/or radiological confirmation of occurrence of any major arterial (e.g. stroke) or venous thromboembolic (deep venous thrombosis, pulmonary embolus) event was captured for each hospitalized Oxford patient and defined as a single binary outcome. Radiological diagnosis of pulmonary embolus was defined separately but compared to this definition to ensure consistency.

Co-prescription of agents hypothesized or proven to be effective in the management or treatment of COVID-19 were logged for all patients with COVID-19. These agents included dexamethasone, remdesivir, interferon 1beta, tocilizumab, anakinra, azithromycin, convalescent plasma or other experimental immunomodulatory agents.

Patients <18 years old or those with active malignancy or receiving significant immunosuppression (greater than an equivalent of 40mg once a day of prednisolone) prior to admission, or those with a clear alternative cause for symptoms and hospital presentation were excluded from analyses. For most modalities, samples were prioritized following a stepwise algorithm to match for age, sex and severity of illness. The steps for matching were as follows: Select severe/critical hospitalized COVID-19 patients first 1. As many critical samples were identified at the earliest point of their maximal severity and samples from patients sampled during severe disease (not defined as frail) were matched in terms of age, sex and time between onset of symptoms and sampling as closely as possible. 2. Rejection sampling was used to select individuals from the comparator groups (mild COVID-19 including both hospitalized and community healthcare workers), healthy volunteers, influenza and all-cause sepsis) to match confounder distributions.

Some downstream analyses required multiple samples from the same individuals transiting through different disease states. In cases where data was available from multiple samples for the same individual a cross-sectional analysis involved prioritizing a single sample based on closeness to maximal severity of illness.

Blood sample processing Whole blood from hospitalized Oxford patients, healthcare workers and healthy volunteers were sampled into Tempus tubes (Life Technologies Corporation) for extraction of whole blood total RNA for sequencing or DNA for genotyping or sequencing, or EDTA buffered vacutainers (Fisher Scientific) for processing within 4 hours of sampling. Processing consisted of mixing of whole blood in a 1:1 ratio with Cytodelics (Cytodelics), and then fractionated using Leucosep (Griener Bio-One) into peripheral blood mononuclear cells (PBMC) and EDTA-buffered plasma. PBMCs and plasma were frozen and thawed in batches for specific experiments to minimize the risk of batch effects. Tempus tubes were frozen at -80C until extraction of RNA performed in batches.

On the day of staining, Cytodelics stabilized samples were thawed, processed to remove red blood cells and fixed using Whole Blood Processing Kit (Cytodelics) as per manufacturer instruction. Two control samples from 2 healthy volunteers (Control A and Control B) were included to each batch to assess batch to batch variations. Fixed cells were distributed in 96 deep well V-bottomed plate, washed once with CSB containing Benzonase, and then barcoded using the Cell-ID 20-Plex Pd Barcoding Kit (Fluidigm).

Briefly, cells were washed once with Barcode Perm Buffer, resuspended in Barcode Perm Buffer supplemented with Heparin to reduce nonspecific eosinophil staining artifacts (Rahman et al., 2016) and loaded with half the amount of metal barcodes recommended by the manufacturer. After 30m incubation at room temperature, cells were washed twice in CSB, pooled and counted. 40-70% of cells in whole blood are CD66 + granulocytes, and their frequency can increase to up to >95% in septic patients. To maintain both the possibility to measure changes in the frequency of granulocytes and to analyze with a reasonable efficiency mononuclear cells in samples from COVID-19 and sepsis patients, each batch was split in two aliquots. The first aliquot was stained with no enrichment/depletion done (i.e. unaltered). The second aliquot was enriched for mononuclear cells using a granulocytes depletion kit based on the magnetic separation of CD66+ cells using beads coated with anti-CD66 antibodies (EasySep; 17882). Both granulocyte depleted and non-depleted samples were first stained with the surface antibody cocktail for 30m at room temperature), washed with CSB, then fixed with Nuclear Antigen Staining Buffer and permeabilized with Nuclear Antigen Staining Perm before staining with an intracellular antibody cocktail. After staining at room temperature for 45m cells were washed in CSB buffer, and fixed in 1.6% FA solution before overnight incubation with Cell-ID™ Intercalator-Ir. On the day of acquisition, cells were washed once in CSB buffer supplemented with Benzonase, resuspended in water and run at an acquisition rate of <300 events/second and acquired on a Helios CyTOF machine. An additional aliquot of whole blood was thawed and stained for those samples where yield was too low after staining and acquisition. In total 31 samples were acquired twice (in different batches).

For flow cytometry we analyzed one aliquot of the same PBMC samples at the same time they were being processed for 10x Genomics CITE-seq. Each aliquot was split into three (to have 5-10x10 5 cells/staining) and distributed in 96 well V-bottom plates. Cells were spun at 500g for 10 min, then pellets stained form 10min at RT in PBS with 15l live dead Aqua, diluted as per manufacturer's instructions) supplemented with 165g/well human Ig (Gammanorm, Octopharma), to block non-specific Fc receptors binding. Cells were washed with 200l FACS buffer (PBS 10% FCS), 10min 300g and pellets were stained with each antibody mix in Brilliant Stain Buffer (BD Biosciences), 20min at RT. Cells were with 200l FACS buffer, 10min 300g and fixed in IC Fixation Buffer (Thermo Scientific) diluted 1:2 in PBS, total 200l/well. Cells were fixed a minimum of 30min and a maximum of overnight before being acquired at the BD Symphony X50.

Total RNA-seq was performed with libraries prepared by Oxford Genomics Centre with the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina after rRNA and globin depletion. Libraries sequenced as a single pool of 144 samples (124 patients) on one NovaSeq S4 flow cell (4 lanes) with a target of 50M 100bp read pairs per sample. Due to limited material availability influenza samples were sequenced used a small-bulk RNA-sequencing method in order to facilitate genetic demultiplexing of the CITE-seq data.

Bulk isotype-resolved B cell receptor sequencing was performed on RNA from whole blood using a protocol adapted from (Bashford-Rogers et al., 2019) , with moderations including the use of a more sensitive reverse transcriptase (Superscript IV) and primer concentrations (primer sequences provided in Key Resources Table) . Bulk TCR sequencing was based on the same protocol adapted from (De Mattos-Arruda et al., 2019) , in which the TRB primers were redesigned to capture the full productive repertoire (primer sequences provided in Key Resources Table) . Sequencing libraries were prepared using Illumina protocols and sequenced using 300bp pairedended sequencing on a MiSeq (Illumina).

To minimize batch effects, cryopreserved PBMCs were processed in ten batches of 14 samples per batch, with each batch containing at least one sample from each comparator group and similar distribution of patients age and sex across the batches. After thawing, all 14 samples were mixed together to form a single pool (at a ratio that yielded the same number of live cells pooled per sample based on live/dead counting) followed by viability staining (7-AAD dye, BioLegend 420404) and live/dead sorting on a BD FACSAria Fusion sorter.

Post sort, each pool was incubated with FcX block (BioLegend) for ten minutes on ice, washed, and stained with a 192 TotalSeq-C antibody panel (BioLegend 99814) for 30 minutes on ice. Cells were washed three times in PBS + 1% BSA, counted on the BioRad TC20 Automated Cell Counter and loaded onto the 10X Genomics Chip G at 50,000 cells per channel. Each pool was loaded across seven channels.

scVDJ-seq data generation scVDJ data was generated on the same cells as the scCITE-seq and GEX datasets. Cell Ranger software (v3.1.0) was used to process the Chromium scRNA-seq output data. The FASTQ files were filtered per sequence library plex (i.e. per-pool) using the 10x Genomics index-hopping-filter (https://support.10xgenomics.com/docs/index-hopping-filter) that implements a strategy to mitigate the known index hopping issue with the Illumina machines that use patterned flow cells. The IMGT's reference genome was used as a reference for the BCR and TCR VDJ libraries. Cell Ranger was also used to filter read alignments, cell barcodes and UMIs, excluding those droplets with low numbers of reads (e.g. erythrocytes or droplets with encapsulating ambient RNA from dying cells).

10x Genomics scATAC-seq Similar to the CITE-seq approach, cryopreserved PBMCs were processed in two batches of five samples per batch with each batch containing at least one sample from each comparator group. Each sample underwent viability staining (7-AAD dye, BioLegend 420404) followed by live/dead sorting on a BD FACSAria Fusion sorter. Post sort, five samples of an equal number of cells were mixed together to form a single pool, and each pool underwent cell lysis and nuclear extraction according to the 10X demonstrated protocol Nuclei Isolation for Single Cell ATAC Sequencing CG000169 Rev D. Briefly, 200,000 cells from each sample were added to form the pool and cell lysis was performed with 100 μl chilled Lysis Buffer (10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl2, 0.1% Tween-20, 0.1% Nonidet P40 Substitute, 0.01% digitonin and 1% BSA) for three minutes on ice. The lysis reaction was quenched with 1 ml chilled Wash Buffer (10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl2, 0.1% Tween-20 and 1% BSA) and the nuclei were centrifuged (500g for 5 min at 4 °C). After removal of supernatant, the nuclei were resuspended in chilled diluted nuclei buffer (10X Genomics; 2000153) and concentration determined with the Bio-Rad J o u r n a l P r e -p r o o f Page 32 of 92 TC20 Automated Cell Counter. Each pool of nuclei was loaded across four channels of the 10X Genomics Chip E (15,000 nuclei per channel).

Sample preparation of non-depleted plasma/serum for LC-MS/MS analysis Samples Plasma and serum samples were stored at -80C then thawed overnight at 4C before use.

For the bottom-up proteomics approach, proteins were precipitated using Isopropanol (IPA) and loaded into the S-trap 96-well plates (Profiti, Huntington, NY, USA) where proteins were retained for subsequent trypsin digestion following S-trap manufacturer instructions (Zougman et al., 2014) . Briefly, five microliters of plasma or serum were pipetted under a category 2 fume hoods into standard 96-well plates pre-filled with 200 µl of isopropanol following the 96-well plate layouts. Precipitation plates were kept at room temperature for 1h before transferring the isopropanolplasma (serum) mixture into the 96-well S-trap plates. 200 µl of S-trap binding buffer (90% methanol, 100mM Triethylammonium bicarbonate (TEAB)) were added into each well of the precipitation plate to recover any protein precipitates and transferred into the corresponding 96well S-trap plates. Plates were then spun at 1500g for 2 min. Protein colloidal precipitates were retained into the S-trap mesh. S-trap plates were then subjected to four consecutive washes with S-trap binding buffer (200 µl for the first wash and 350 µl the three remaining), followed by a 2 min spin at 1500g. 125 µl 50mM TEAB containing trypsin (Worthington) to a 1:25 wt:wt ratio was added into each well and incubated overnight at 37C. Tryptic peptides were sequentially eluted from the S-trap into a clean 96-well plate with 80 µl of 50mM TEAB, 80 µl of 0.2% formic acid (FA) in LC-MS water and 80 µl of 0.2% FA in 50% acetonitrile with 2min spins at 1500g after each solvent addition. Finally, plates containing tryptic digests were dried in the speed vac and reconstituted with 70 µl of 0.1% FA for subsequent LC-MS/MS analysis after two serial dilutions (1/100 and 1/10) in 0.1% FA for a final 1/1000 dilution plate to achieve suitable peptide concentration for LC-MS/MS loading.

For quality control (QC) purposes and for repeat injections, a pool for each clinical sample group was created and processed as described above. In addition to the sample group pools, an overarching master pool for each of the plasma and serum cohorts was also built from the individual non-diluted tryptic peptide sample group pool with master pools built to reflect the contribution of each individual group to the overall study.

To increase the depth of the non-depleted plasma/serum proteome, master pools were subjected to a high pH fractionation using the SOLA HRP cartridges (Thermo Fisher) with the aim to create a protein library. In brief, an equivalent of 150µg of peptides from each master pool was fractionated into eight different fractions going from 0% to 70% acetonitrile in 10mM ammonium formate at pH 10. Prior sample loading, SOLA HRP columns were conditioned with 0.1% trifluoroacetic acid (TFA) in 70% acetonitrile and washed with 0.1% TFA. Following the sample loading, columns were washed with 0.1% TFA before starting the stepwise high pH fractionation elution with 200µl of 0%, 10%, 15%, 20%, 25%, 30%, 40% and 70% ACN in 10mM ammonium formate pH 10. Eluted fractions were then dried in a speed vac and reconstituted with 40µl of 0.1% FA prior to LC-MS analysis.

C18 Evotips (Evosep one, Odense Denmark) were prepared following manufacturer's instructions. Briefly, Evotips were conditioned with 1-propanol, washed with solvent B (0.1% FA in 100%

Page 33 of 92 acetonitrile) and equilibrated in solvent A (0.1% FA in LC-MS water) before loading them with a total volume of twenty microliters containing 10 µl of sample (peptides) from dilution 1/1000 plate (equivalent to an estimated 30 ng of peptides). Evotips were washed once with solvent A and stored in 100 µl solvent A till the sample injection.

Samples were analyzed using an Evosep One LC system connected to the TimsTOF Pro mass spectrometer (Bruker Daltonics). Peptides were analyzed using the pre-built 100 samples / per day method (Evosep) with an 11.5 min gradient (total cycle time of 14.4min) at a 1.2 µl/min flow rate (Bache et al., 2018) . In brief, tryptic peptides were transferred from the pre-loaded C18 Evotips with a pre-built gradient to a sample loop and separated on an 8cm C18 analytical column (Evosep Pepsep, 3um beads, 100 um ID) with an overall gradient from 3 to 40% acetonitrile.

Mass spectrometry data were acquired in PASEF mode (oTOF control v6.0.0.12). The ion mobility window was set to 1/k0 start = 0.85 Vs/cm2 to 1/k0 end = 1.3 Vs/cm2, ramp time 100 ms with locked duty cycle, mass range 100 -1700 m/z. MS/MS were acquired in 4 PASEF frames (3 cycles overlap). Target intensity was set to 6000 and threshold intensity 200. The data acquisition method has been deposited within the raw data in the ProteomeXchange data repository (PXD023175) (Perez-Riverol et al., 2019) .

To assess the technical reproducibility, master pools were run before each sample group and sample group pools were run after the master pools, every twenty sample runs and after the last sample of the group. Blanks were injected in between groups to control for carry over.

The concentrations of selected proteins in the plasma and serum were measured with Human Magnetic Luminex Kits (Bio-techne) with 3 panels containing total 51 analytes: C-C motif ligand (CCL)2/3/4/11/17/18/19/20, CD40 Ligand (CD40L), CD163, complement component 5a (C5a), C-X-C motif chemokine ligand (CXCL)1/5/10, epidermal growth factor (EGF), basic fibroblast growth factor (FGF2), granulocyte colony-stimulating factor (G-CSF), granulocyte-macrophage colonystimulating factor (GM-CSF), granzyme B (GrB), interferon (IFN)a/b/g, interleukin (IL)-1a/1b/2/3/5/6/8/10/12/13/15/17A/23/33, lactoferrin (LF), Lipocalin-2 (LCN2), Lymphotoxin-alpha (LT-a), macrophage colony-stimulating factor (M-CSF), Myeloperoxidase (MPO), beta-nerve growth factor (b-NGF), Oncostatin M (OSM), S100 calcium-binding protein A9 (S100A9), stem cell growth factor (SCGF), tissue factor (TF), tissue factor pathway inhibitor (TFPI), transforming growth factor alpha (TGF-a), Thrombopoietin (THPO), tumor necrosis factor (TNF) and triggering receptor expressed on myeloid cells 1 (TREM-1). The assays were conducted according to the manufacturer's instruction. Results were obtained with a Bio-Rad Bio-Plex® 200 Systems. Fluorescence intensity (FI) data from the assays were used for further analysis.

To analyze demographic and clinical features cohorts the following statistical tests were applied. For approximately gaussian distributed continuous data we used the anova test (using the function f_oneway from python's scipy.stats package). For other distributions of continuous variables, we used the Kruskal-Wallis-test (using the function kruskal from python's scipy.stats package). For categorical data we used the chi-squared test (using the function chi2_contingency from python's scipy.stats package).

Unsupervised and consensus k-means clustering were used to analyze the clinical data (Moniti et al., 2003; Seymour et al., 2019) . We prepared the data by applying quantile-normalization to counts and shift+rescaling to zero mean and unit variance to all features (applying quantile normalization to all features, not just counts, further improves the correlation of the ordered consensus matrix with the WHO severity_at_sample classes). We then ran consensus k-means clustering to identify the best number of clusters. This means we sample from our clinical dataset 50% of the entries, run k-means clustering (where we vary k, i.e. the number of clusters), and record how often two elements were clustered together out of the time where they were in the 50% of samples data points in a so-called consensus matrix. We repeated this process 600 times for k=2,3,4,5,6 clusters. For each k, we computed the empirical cumulative distribution (CDF) and the change in the area under the CDF curves for different k. The optimal cluster number is obtained for the k after which the CDF curves (and hence the areas) do not change significantly anymore. We then performed an unsupervised consensus k-means clustering analysis of (subsets) of the clinical variables (using python's sklearn library). In order to mitigate the question of which features are independent, we performed this analysis on the principal components (using sklearn). After consensus clustering, we sorted the consensus matrix using hierarchical clustering. The methods were implemented in python and used the packages numpy, scipy and pandas. We also performed hierarchical clustering directly on the PCs, rather than on the consensus matrix (using AgglomerativeClustering and dendrogram from sklearn).

CyTOF data pre-processing After acquisition, data was normalized and concatenated on CyTOF Software v7.0, compensated on CATALYST (version 1.5.3.23). Then, data was further processed for removal of beads and of DNA negative events and for Gaussian parameters using a manual gating strategy using FlowJo v10.6 as described in Data S3. Then, samples were manually debarcoded (FlowJo v10.6) and CD45 + events (non-depleted samples), or CD66-Siglec8-CD45 + (granulocytes-depleted samples) events were selected for further processing. For each batch, two control samples were run together with the patients' samples to allow for correction for batch effects. Batch correction was performed with CytoNorm software (version 0.0.5) separately for granulocyte depleted and nondepleted samples. For both analyses, the training of the algorithm was done with 20,000 cells per sample and 5 clusters (parameter nClus) using one of the control samples and the results were assessed using the other control sample. Doublet cells were removed by manual gating on FlowJo v10.6 (Data S3) and were deemed as Iridium high Ki67cells to avoid removing proliferating cells. Finally, to correct for potential biases due to highly varying cell numbers per sample, downsampling was performed to a maximum 75,000 cells per sample for the granulocyte depleted samples and 40,000 for the non-depleted samples.

Clustering analysis and the identification of the different immune cell population was done using the analytical pipeline described in Data S3. Clustering analysis was performed separately for granulocyte depleted and non-depleted samples using a self-organizing map algorithm through an implementation of FlowSOM via the CATALYST R package (version 1.10.3). The initial clustering was performed using a resolution of 144 clusters (dimensions 12x12) and metacluster merging by consensus clustering was performed at a resolution of 25 clusters for both datasets. The clusters were then manually annotated to define the populations. Three major populations (T/NK, B/Plasmablasts and Myeloid) from the granulocyte depleted samples were further clustered to identify finer subpopulations using a resolution of 225 (dimensions 15x15) for T/NK and 144 (12x12) for the other two. The metacluster merging was performed at a resolution of 50 clusters for the T/NK population, 25 for the B/Plasmablast and 30 clusters for the Myeloid Page 35 of 92 population. Clusters were manually annotated and further managed to provide cell type resolution at two different depths -one to define major cell types, and a finer one to define subpopulations. The cell type assigned to each cell ID in each of the three major population was then used to reconstitute the frequency of each population inside the granulocyte depleted sample. The finer annotation was used to evaluate the frequency of cell subtypes.

Differential abundance analysis was performed using code from the diffcyt package (version 1.8.8) with the option testDA_edgeR. The analysis was adjusted for the confounding variables batch, age and sex. For the analysis of the non-depleted samples the normalization for composition bias was deemed necessary to account for the pronounced differences in neutrophil proportions. For all other analyses no normalization was performed. The subpopulations of the granulocyte depleted samples were analyzed independently of each other to avoid composition effects.

The full results of the clustering and differential abundance analysis results are shown for nondepleted whole blood together with principal components analysis (Data S3); and clustering and differential abundance analysis for granulocyte depleted whole blood broad populations, myeloid, B cell and T/NK populations (Data S3).

Two subpopulations of cells were more precisely defined by manual gating rather than by clustering. To do that, the cluster of the parent population (cMono or MAIT) was exported from the R environment and gated in FlowJo. Also, the median intensity of CD38 expression were identified in different cell types within the granulocyte depleted samples.

The density plots were made by down-sampling to have the same number of events across all different conditions and using the geom_pointdensity function of the ggpointdensity package, with an adjust of 4, to generate the ggplot.

The two datasets were aligned using Seurat (v3.9.9.9010) using only samples that had been analyzed on both technologies (115 samples, after exclusion of one sample due to low number of cells in CITE-Seq). The integration was based on an overlapping set of 38 markers that were common between the CyTOF and ADT antibody sets. Both datasets were filtered to exclude cells with unclassified or uncertain annotations. The CITE-Seq dataset was further filtered to exclude potentially low-quality cells with i) number of genes, ii) number of features of ADT, iii) number of total UMI (RNA) or iv) number of total UMI of ADT lower than 0.001 of their relevant distributions. The two datasets were then down-sampled to 1,000 cells per sample and the integration was performed on 230,000 cells in total.

The anchors between the query and reference datasets were found based on a CCA with 30 dimensions and the labels were transferred using a PCA for the weight reduction of the query dataset. The analysis was performed with either the CITE-seq or the CyTOF as reference datasets in turn and all cells were given a predicted annotation from their counterpart dataset. Then the predicted annotation of the cells was compared to their original annotation and all the common major cell types were validated with very high accuracy ( Figure S2B) . Finally, the visualization of the two datasets was performed with the CITE-seq data as a reference. The anchors were used to impute the ADT markers for the CyTOF cells which were then merged with the CITE-seq data and centered in order to run a UMAP analysis to visualize all the cells together ( Figure S2A ).

Data were analyzed with FlowJo version 10. Cells were gated on live leukocytes and recompensated in FlowJo. Frequencies of individually gated populations were exported and plotted in Prism. Alternatively, 15,000 CD3 + cells per sample were concatenated and exported for subsequent analysis in R. Supporting data for analysis of memory CD4 + T cell subsets is shown (Data S3).

We trimmed adaptor sequences using TrimGalore (v0.6.2, https://github.com/FelixKrueger/TrimGalore), and aligned reads to the reference genome (GRCh38.100) using multi-sample 2-pass mapping with STAR v2.7.3 (Dobin et al., 2013) (ENCODE Best Practices recommended parameters). We quantified gene expression using featureCounts (v1.6.4) (Liao et al., 2014) and annotations from Ensembl (v100). We calculated and checked QC metrics from FastQC, mapping metrics from STAR, duplication rates and other RNASeq metrics from Picard (v2.23) (http://broadinstitute.github.io/picard/) and RNASeQC (v2.3.5) (DeLuca et al., 2012) and checked for outliers in principal component analysis. We removed 1 sample with high proportions of duplicates and short reads, low exonic rate, number of mapped reads and number of genes identified, and which was also clearly outlying on PCA. We also confirmed that there were no sex mismatches that could indicate sample mix-ups based on sex chromosome gene expression using QTLTools (Delaneau et al., 2017) . This resulted in a data set of 143 samples from 123 patients. We filtered out features that did not have at least 10 reads in at least 10 samples, retaining 23,063 features for downstream analysis. We then normalized the data using the trimmed mean of M-values method R package from edgeR; https://doi.org/10.1093/bioinformatics/btp616 and log2-transformation of counts-per-million. Code is available on GitHub (https://github.com/COMBATOxford/bulkrnaseq/mapping/).

We carried out principal component analysis (PCA) on the normalized filtered data for 123 patients using prcomp (R v 3.6.2) with default parameters. We additionally performed PCA on just the samples from hospitalized COVID-19 patients (n=64). We selected informative PCs by taking the first n PCs that cumulatively described at least 80% of variance in the data. We explored the relationship between these PCs and the clinical variables available (transformed and filtered as described above) using Spearman correlation ( Figure S3A) . Additionally, we investigated the biological relevance of the PCs through pathway enrichment analysis with the R package XGR (Fang et al., 2016b) and Reactome pathways. For each PC, we took the 500 genes with the highest absolute loading scores as input and compared to a background of all detected genes.

For PCA plots, to aid visualization 95% data ellipses were generated assuming a multivariate tdistribution. These were drawn using the stat_ellipse method from the R package ggplot2, using default parameters. Data ellipses were prepared in the same fashion for CITE-seq composition and expression analysis, and proteomics.

We performed differential expression analysis on the normalized data with one sample per patient using the limma R package (https://doi.org/10.1093/nar/gkv007). We included age, age 2 and sex as fixed effects and compared each patient subgroup with all others, as well as investigating Page 37 of 92 COVID-19 severity as a quantitative trait (mild=1, severe=2, critical=3). We defined significance for downstream analysis as FDR <0.05 and fold change >1.5. We investigated the impact of variation in cell proportions across samples calculated from hospital measurements for neutrophils, monocytes, and lymphocytes, by adding neutrophil and monocyte proportions to the model and comparing the estimated fold changes to the results to the more basic model (correlation plot showing the influence of cell proportion on detection of differentially expressed genes in Figure S3C ). Pathway enrichment analysis was performed using Reactome pathways via the XGR R package (Fang et al., 2016b) , with Fisher's exact test and filtering of redundant terms by the xEnrichConciser function.

We applied weighted gene correlation network analysis (WGCNA) to describe modules of highly correlated genes within the whole blood total RNA-sequencing data (i.e. 143 samples from 123 patients) (Langfelder and Horvath, 2008) . We utilized the "cornet" pipeline to wrap the WGCNA R package and perform gene set enrichment analyses (https://github.com/sansomlab/cornet.git). In brief, a stepwise approach of correlation network construction and module detection was adopted using log2-transformed counts-per-million RNA-sequencing data. Using the soft thresholding power of 4, a signed-hybrid network was built implementing the biweight midcorrelation as the adjacency function. The adjacency matrix was transformed into a topological overlap matrix in order to calculate the dissimilarity and a dissimilarity threshold of 0.3 was used to merge modules with very similar expression profiles. Module eigengene values (module first principal component) were used to summarize modules and perform module-trait correlation analyses (Pearson correlation). Gene set over representation analyses were performed using the default settings of the "cornet" pipeline (https://github.com/sansomlab/cornet.git).

Pre-processing of 10X Libraries (Gene expression, ADT, TCR, BCR) As described above a total of n=140 PBMC samples from COVID-19, sepsis, influenza and healthy volunteers were mixed into n=10 pools. Each pool comprised of n=14 samples (each from a different individual) and cells from each pool were captured using n=7 10X channels ( Figure S1D) . From each channel, four libraries were generated, Gene expression (GEX), Surface proteome (ADT), TCR and BCR repertoires. The libraries were sequenced on an Illumina NovaSeq 6000 sequencer across nine S4 flowcells (4 lanes per flow cell). For each of the sample pools, the n=7 GEX and ADT libraries were sequenced on 3 dedicated lanes (to prevent index hopping between sample pools). The BCR/TCR libraries were sequenced on the remaining 6 lanes. Mapping indexes for the bulk and single cell data were built using GRCh38 genome sequences and Ensembl version 100 annotations. For construction of the Cellranger reference transcriptome we included genes with the following Ensembl biotypes: protein_coding, IG_V_gene, IG_V_pseudogene, IG_D_gene, IG_J_gene, IG_J_pseudogene, IG_C_gene, IG_C_pseudogene, TR_V_gene, TR_V_pseudogene, TR_D_gene, TR_J_gene, TR_J_pseudogene, TR_C_gene (n=20,615 genes; lncRNA not included). FASTQ files were generated using Cellranger (v3.1.0) mkfastq and individually QC'ed using FastQC (Andrews, 2010) . To mitigate index hopping between channel libraries within sample pool FASTQ files were filtered per pool using the 10x Genomics index-hopping-filter (https://support.10xgenomics.com/docs/index-hopping-filter) (v1.0.1). Per-channel (n=70) GEX and CITE-seq matrices were prepared using Cellranger (v3.1.0) count. Cell identification was performed on the GEX modality with both Cellranger (v3.1.0) and EmptyDroplets (Lun et al., 2019) (1.8.0) and the union of the calls from both algorithms taken forward as the set of identified cells.

The GATK variant calling pipeline (https://github.com/gatk-workflows/gatk4-rnaseq-germlinesnps-indels) (v4.1.7.0) was applied to the bulk RNA-seq data to identify sequence variants for the genetic demultiplexing of the single cell data. Minor modifications were made to the variant calling pipeline to allow execution from FASTQ files, incorporate metadata into the BAM alignment files, and produce a block gVCF file. We applied both Demuxlet (Kang et al., 2018 ) (v2; https://github.com/statgen/popscle/commit/038044660337b50ec89b0b355493ceb707f18ecd) and

Vireo (Huang et al., 2019) (v0.4.0; https://github.com/single-cellgenetics/vireo/commit/4aecb54a4f2fa3a3413e2cd3e9f3515744385cba), proceeding with the demultiplexing calls from Vireo as they were more robust to variation in read depth. The dataset was demultiplexed using full genotypes (per-pool gVCF files) from the bulk sequencing. For each of the n=70 10x channel sequence libraries, Vireo consistently demultiplexed 60 -75% of cells. A total of n=884,587 singlet cells were demultiplexed.

The pre-processing workflow was encapsulated in a set of CGAT python pipelines (Cribbs et al., 2019) and version controlled. Cell quality was assessed with a suite of metrics that included total UMI (GEX), number of genes detected, percent mitochondrial gene expression, percent ribosomal gene expression, percent IgG expression, Scrublet doublet score (Wolock et al., 2019) and total UMI (ADT). Ambient RNA was assessed. Cell QC statistics, demultiplexing assignments and cell and patient metadata were centrally warehoused in an sqlite database. Following inspection of the QC metrics, the dataset was filtered to retain cells with ngenes > 300 and pct_mitochondrial < 10%. In total, n=836,148 cells were selected for downstream analysis. Expression data for selected cells was extracted from the per-channel count matrices using R and combined into a single market matrix with AWK. RNA velocity matrices were computed for each of the channels with Velocyto (http://velocyto.org/) (La Manno et al., 2018) (v0.17.5) . Velocity data for selected cells was extracted and combined into a single matrix using the Loompy Python library (http://loompy.org).

We used expert immunological knowledge to guide a curated integration of the data from the different modalities (GEX, ADT and VDJ) to identify and label the cell sub-populations present in the CITE-seq dataset ( Figure 1B) . As detailed below, we first performed separate clustering of gene expression, clustering of surface protein expression, and analyses of T and B cell receptor V(D)J sequences. Next, led by expert understanding of the three feature spaces we prioritized use of ADT surface phenotype for definition of major cell lineages and subsets where definitive marker expression was available. Cell types and subsets were further refined using information from the repertoire and GEX layers, or in the absence of definitive ADT information were identified by GEX cluster phenotype. Finally, the identified cell types and subsets were further divided by inferred functional state based on targeted assessment of information from all three modalities. For example, cell cycle phase was determined by GEX phenotype, T cell memory vs effector status was distinguished using information from both the GEX and ADT layers, while assignment of B cell maturation status involved use of information from all three modalities (including BCR mutational status). Information from all three modalities was used to identify and exclude doublets from downstream analysis.

Prior to alignment data was normalized and log1p transformed using Scanpy (Wolf et al., 2018) . The top n=3000 highly variable genes (HVG) were selected (Seurat flavor with n_top_genes set). Effects associated with total number of UMIs were regressed out and PCA components computed. We evaluated alignment with Harmony (Korsunsky et al., 2019) , Scanorama (Hie et al., 2019) and BBKNN (Polanski et al., 2020) and chose to proceed with Harmony following inspection of UMAP projections of the alignment results (data not shown). Following inspection of the PCA scree (knee/elbow) plot, Harmony alignment of samples (n=140 levels) was performed in Python using the top n=65 PCs. Leiden clustering, marker discovery (wilcox tests) and visualization (with violin plots, UMAPs, volcano plots, MA plots and heatmaps), automatic cell type identification (with singleR (Aran et al., 2019) and via overrepresentation analysis of xCell (Aran et al., 2017) genesets), basic composition analysis, geneset over-representation analysis (including of GO, KEGG and Biocarta genesets) and visualization of the aligned datasets was performed using pipeline_scxl.py (https://github.com/sansomlab/tenx). Neighbor graphs were built using Euclidean distance and the "HNSW" algorithm (as implemented in ScVelo (Bergen et al., 2020) , following the example of Pegasus ). Clustering was performed using Scanpy, marker discovery with Seurat (Stuart et al., 2019) and pathway overrepresentation analysis with gsfisher (Croft et al., 2019) (https://github.com/sansomlab/gsfisher).

After alignment of the full manifold, we iteratively divided, re-aligned and re-clustered the dataset to achieve high-resolution clustering of different cell subsets. Variable gene identification was performed separately within each subset as described above. In the first step the major cell types -T/NK, myeloid and B/plasmablast cells were extracted into n=3 separate subsets. Each of these three subsets was then separately aligned and clustered (n=45 PCs for n=482k T/NK cells; n=45 PCs for n=279k myeloid cells; n=40 PCs for n=66k B/plasmablast cells). This process was then repeated with a final ( Figure  S1E . Choice of principal component number was guided by inspection of the scree plots.

Each of the cell six subsets was subject to Leiden clustering with a range of resolutions (1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 6, 8) and cluster phenotypes assessed using pipeline_scxl.py (https://github.com/sansomlab/tenx) as described above, and ADT surface phenotype. For each subset a base resolution was selected and sub-populations further refined by manual splitting or merging of clusters based on alternate cluster resolutions guided by expert knowledge of cell identities and phenotypes to define a total of n=131 GEX clusters. These clusters were intersected with information from the ADT and repertoire modalities for the final multimodal annotation of the dataset as described below.

For each channel the ADT signal was normalized independently to its background, using the DSB algorithm (Mulè et al., 2020) taking the set of cells identified from the GEX data as the foreground and unassigned GEMs with log10(GEX_nUMI +1) >=1.5 as the background. The workflow was written in R and Python, encapsulated in CGAT pipelines and version controlled. The normalized ADT data was subject to hierarchical stochastic neighbor embedding (hSNE) and unsupervised Gaussian mean shift (GMS) clustering as implemented in Cytosplore (van Unen et al., 2017) (www.cytosplore.org). This analysis, which was based on the expression of all n=192 ADT tags, partitioned the dataset into 66 discrete clusters. Based on expression of known lineage markers (CD3, CD4, CD8, CD56, CD19, CD20, CD27, CD38, CD14, CD123, CD1c, CD33, CD235, CD34, Vd2, Va72 and CD161) these clusters were apportioned to NK cells, B cells, Plasmablasts, Vd2+

plasmablasts were defined by re-clustering these cells based upon expression of CD10, CD24, CD25, CD19, CD20, CD22, HLADR, CD21, CD23, IgD, CD1c, IgM, CD38, CD29, CD71, CD39, CD27, HLA-ABC and IgA. In total, this analysis of the ADT data identified 12 T-cell subpopulations, 7 B cell/plasmablast subpopulations and five distinct sets of doublets (with markers of more than one major immune lineage). These clusters were intersected with information from the GEX and repertoire modalities for the final multimodal annotation of the dataset as described below.

T and B cell V(D)J gene usage and receptor sequences were quantified using Cellranger VDJ (v3.1.0) with reference sequences from the IMGT. For TCR chain usage we recorded additional information for each cell including (i) "clone proportion", computed as the proportion of the entire T cell repertoire occupied by the clone in the sample, (ii) whether the cell was a doublet (defined as a singlet clone with TRA-TRA-TRB-TRB chains), (iii) whether the cell possessed an iNKT receptor phenotype (TRAV10 with TRAJ18), and (iv) whether the cell carried a MAIT receptor phenotype (TRAV1-2, with either TRAJ12/TRAJ20/TRAJ33). Immunoglobulin sequences were further analyzed using IMGT High V-Quest. B cell/plasmablast doublets were identified as cells with multiple heavy chain (HC) contigs where the ratio of the number of UMIs for the top two ranked HC contigs was >0.125. Light chain (LC) information was not used for doublet identification as it is known that 1% of B cells may express dual light chains. The TCR and immunoglobulin chain usage information was intersected with information from the GEX and ADT modalities for the final annotation of the dataset as described below.

As described below, we performed curated intersections of the information from the different modalities to discern the cellular identities and functional phenotypes of the different PBMC populations. Information from all three modalities was used to screen and exclude doublets: in addition to inspection of Scrublet scores (computed from the GEX layer), lineage according to surface phenotype, immune receptor sequence and gene expression profile was required to be congruent. In total, following curated intersection of GEX, ADT and repertoire information we obtained n=853 clusters. Of these we retained n=346 clusters (n=787,928 cells) for downstream analysis. The retained cells comprised of n=343 clusters (n=782,245 cells) that were assessed to represent genuine singlet cell clusters based on expert knowledge. We assigned each of these clusters to a "minor subset", "major subset" and "cell type" category in order to accommodate downstream analyses at different levels of granularity as summarized in Figure 1C . In addition, we retained two clusters that appeared to comprise monocyte-platelet doublets (n=5562 cells) and a minor, lower confidence cluster of AXL+ DCs (n=121 cells) for inclusion in cluster-level analyses. We excluded n=92 clusters (n=20,391 barcodes) that likely represented cell-cell doublets/multiplets and n=286 clusters (n=19,858 cells) which appeared to be comprised of mixtures of cell types. A further n=129 clusters (n=7,971 cells) of uncertain phenotype were also excluded. Scripts were written in Python and R and version controlled.

Known T and NK cell subsets were identified by surface protein (ADT) and immune receptor (TCR) phenotype as recorded in Table S2 . These identity assignments were then intersected with the gene expression clusters (GEX) to identify subsets of the different lineages that displayed distinct transcriptional phenotypes (e.g. naïve vs memory T cells). These sub-populations are shown in Data S4 for CD4 T-cells, for CD8 T-cells, for double positive (DP) T-cells, for double negative (DN)

B cell and plasmablast cell subsets were defined using information from the GEX, ADT and immunoglobin repertoire modalities as detailed and shown in Data S4. Sub populations of innate immune cells and other blood cells were delineated according to gene expression clusters and annotated based on their GEX and ADT phenotypes. The identified sub-populations are shown in Data S4 for mononuclear phagocytes (MNPs), for hematopoietic stem (and progenitor) cells (HSCs) and for platelets and erythrocytes.

To assess the relative immune receptor diversity of the different T, B and plasmablast populations we computed repertoire Shannon and Gini indices (detailed in the Repertoire analysis section, and as shown in Data S4). Because of the large variance in cell number between clusters we estimated the indices by bootstrapping with a sample size of n=30, n=1000 times. Individual bootstrap samples were drawn without replacement to avoid introduction of artificial clonality. The indices were not estimated for clusters with fewer than n=30 cells.

Three samples with <500 cells in total, and three samples with confirmed or suspected malignancy, were excluded from all composition analyses. For the hospitalized COVID-19 and sepsis clinical categories, only samples closest to maximum severity and only one sample per individual were included, such that for each category the following numbers of samples were analyzed: healthy volunteers, n=10; COVID-19 acute in-patient mild (OUH), n=12; COVID-19 acute in-patient severe (OUH), n=20; COVID-19 acute in-patient critical (OUH), n=18; COVID-19 community COVID-19, n=12; influenza, acute in-patient, n=10; and sepsis acute in-patient, n=15. Scripts for all analysis were written in Python and R and version controlled.

Composition analysis was performed for the different levels of cellular granularity summarized in Data S5, including for cell types; major cell subsets; minor cell subsets; T and natural killer (NK) cell clusters; B and plasmablast (PB) cell clusters; and mononuclear phagocyte (MNPs) clusters. To initially inspect cluster abundance across clinical categories, the percentage frequencies of cell subpopulations were quantified per sample and visualized as boxplots. For cell types, major cell subsets and minor cell subsets (Data S5), percentages were calculated out of total PBMCs. For higher resolution immune cell clusters (Data S5), frequencies were calculated out of total T and NK cells, total B and PB cells, or total MNPs. Selected panels are reproduced in the main and supplementary figures.

Principal component analysis (PCA) was used for exploratory analysis, with pre-filtering to remove clusters with a count of n<10 cells in <10 individuals. Cluster counts were converted to centered log-ratios (CLRs) using the ALDEx2 R package, v1.18.0. Briefly, 1,000 Monte Carlo samples of the Dirichlet distribution were generated from the cluster counts for each sample and converted to CLRs. The median CLR for each cluster-sample combination was used for PCA which was performed using prcomp (R version 3.6.2) with default parameters. Association tests for clinical variables were performed using an omnibus analysis of variance (ANOVA) to test for association between the top 15 PCs (collectively explaining >80% of the total variance) and clinical category (source), age, sex, and sample pool, with the Benjamini-Hochberg procedure for multiple testing correction and a significance threshold of Pc<0.05. For the hospitalized COVID-19 cases, association tests with additional variables that were related to acute disease or were secondary to infection were also performed. All variables had one or no missing values. These additional variables included: ethnicity; weight; days from symptom to sample and from symptom to admission; maximum temperature in the 24 hours preceding sampling; persistent fever in the 24 hours preceding sampling; fever, WHO ordinal, ventilation status, SaO2/FiO2 ratio, and SOFA oxygenation score at the time of sampling; quantile normalized white cell, neutrophil, lymphocyte, monocyte and platelet counts in the 24 hours before or after sampling; quantile normalized highest concentration of C-reactive protein in the 24 hours before or after sampling; clinical and/or radiological evidence of thromboembolism during hospitalization; length of hospital stay; OpenSAFELY COVID-19 mortality propensity score; and death in the hospital.

The full results of the PCA and association test results are shown in Data S5 for the cell types; for the major cell subsets for the minor cell subsets; for the T and NK clusters; for the B and PB clusters; and for the MNP clusters with selected panels reproduced in the main and supplementary figures.

To compare cell subset/cluster abundance across clinical categories, the percentage frequencies of cells were quantified per sample and visualized as boxplots. For major and minor cell subsets, percentages were calculated out of total PBMCs. For higher resolution immune cell clusters, frequencies were calculated out of total myeloid cells, total B and plasmablast cells, or total T and natural killer (NK) cells. To determine the statistical significance of differences in cell subset/cluster abundance between groups, differential abundance analysis was performed using edgeR, v3.28.1 (Amezquita et al., 2020) . Cell subset/cluster counts were modelled adjusting for age, sex and sample pool using a quasi-likelihood negative binomial generalized log-linear model (glmQLFit function), with pre-filtering to remove clusters with a count of n<10 cells in <10 individuals. Counts were normalized by the total number of cells in each sample for the major and minor cell subset analyses, or by the total myeloid, the total B and plasmablast, or the total T and NK cells for higher resolution immune cell cluster analyses. Differential abundance testing was performed using the quasi-likelihood F-test. The Benjamini-Hochberg procedure was implemented to correct for the number of clusters and the number of pairwise clinical category comparisons, and a significance threshold of Pc<0.05 was used. Testing for composition effects did not provide evidence for biases in cluster abundance. In addition, for the hospitalized COVID-19 patients, edgeR analysis was performed using ANOVA to identify statistically significant associations between cluster abundance and patient characteristics/clinical variables (as detailed for the PCA association tests).

The full results of the differential abundance and covariate analysis results are shown in Data S5 for the cell types; for the major cell subsets; for the minor cell subsets; for the T and NK clusters; for the B and PB clusters; and for the MNP clusters.

Data pre-processing Pseudobulk counts were generated for each combination of gene and sample at minor subset, major subset and cell type level by summing together the within-group gene counts. These were then converted to reads-per-million (RPM) by normalizing by the total count across all genes for each combination of sample and cell type. Finally, residuals were calculated by taking the log(1 + count) and subtracting the predicted value from a linear model with pool as the independent variable. The six poorly performing samples mentioned in the Composition section above were removed from all further processing.

PCA was carried out on the residuals using the prcomp function in R. We generated PCs separately at minor subset, major subset and cell type level using residual values from genes with a mean RPM > 1 in among all sources, and with six poorly performing samples removed. Outlying samples (more than 5 SDs away from the mean across the top 10 PCs) were then removed and PCs recalculated. Plots are shown in Data S6 for the first two principal components of pseudobulk gene expression for each of the major subsets. To test differences in major subset PCs across categories, we used an omnibus ANOVA (i.e. a test of a linear model including all diagnoses/severity categories as dummy variables, against a model with no difference between these groups) to test for association between the first 10 PCs and diagnosis/severity (across all samples) and disease severity (measured by WHO criteria, within COVID-19 samples) (Data S6).

The linear model also includes terms for age and sex. Omnibus p-values were corrected for multiple testing using the Benjamini-Hochberg procedure. We also repeated the principal component analysis only within hospitalized COVID-19 samples and tested their association with clinical covariates. Only two major subsets showed adjusted p-values < 0.01, classical monocytes and DCs, with association patterns shown in Data S6 for clinical covariates in hospitalized patients.

We generated UMAPs from the top n=10 PCs for each level of cell clustering (using the R package "umap"). For this analysis we excluded clusters with n=0 cells for 10 samples and excluded samples with n=0 cells in remaining clusters. The resulting UMAP plots are shown in Data S6 combining PCs across all clusters at the minor subset, major subset or cell type levels.

We carried out differential expression tests for a range of contrasts (including subgroups of COVID-19 against healthy volunteers, sepsis and influenza vs healthy volunteers, COVID-19 vs sepsis, COVID-19 vs influenza and COVID-19 subgroups against one another. These tests were performed using the R edgeR library (McCarthy et al., 2012) . We included age, sex and pool as covariates, and filtered out samples with <5 cells and gene with mean RPM < 1, with significant genes selected based on log fold change (logFC) > 2 and false discovery rate (FDR) < 0.01 (calculated by Benjamini-Hochberg). To test for a set of genes that differed across COVID-19 severity categories, we also carried out a "COVID-19 omnibus" likelihood ratio test, with dummy variables for WHO severity and a separate variable for mild recording healthcare workers (in this case, genes were selected if they had logFC > 2 between any pair of categories).

The number of differentially expressed genes (FDR < 0.01, absolute log-fold change > 2) in each cell cluster for each contrast is shown in Data S6, volcano plots of differentially expressed genes for critical COVID-19 patients vs healthy volunteers for each major subset are shown in Data S6, and the top 10 most differentially expressed genes for each contrast (across all cell clusters) in Table S3 . We also investigated the relationship between the frequency of a cell subset and the number of differentially expressed genes (Data S6), and the relationship between differential composition and number of differentially expressed (Data S6).

Gene-set enrichment analysis GSEA of differentially expressed genes were performed using the FGSEA algorithm (Korotkevich et al., 2021) . We performed GSEA for biological pathways from the MSigDB database (including for KEGG, GO, canonical pathways, regulatory sets, immune signatures and Hallmark genesets). Enrichment analysis was carried out separately for each pair of cell cluster and contrast, with genes ranked by p-value. For the Hallmark genesets, we also carried out a signed enrichment test (ranking by log(pvalue) x sign(logFC)), shown as a heatmap for major subsets and a range of contrasts in Figures S4D, S5E and S6H . The top 10 pairs of cell cluster and GO term, canonical pathway and immunologic signature for critical COVID-19 vs healthy volunteers are shown in Table S3 .

In addition, we applied GSEA to a set of experimentally derived interferon-stimulated genes (ISGs) (Schoggins et al., 2011) . Genes within this gene set were differentially expressed across a wide range of minor subsets and contrasts (Figure S4E) , with the most significant enrichment being seen in HSCs in critical COVID-19 patients ( Figure S4F) . We focused in on classical dendritic cells, as this minor subset showed enrichment of differential expression in ISGs (p < 1e-5) in 9 out of 13 contrasts tested. A subset of ISGs were identified both as driving the ISG enrichment signature (leading edge genes) and showing individually strong differential expression (FDR < 0.001, absolute fold change > 3), shown in a hierarchically clustered heatmap of gene expression in classical dendritic cells of highly differentially expressed genes from the leading edges of interferon stimulated gene sets ( Figure S4G ).

We performed separate WGCNA (Langfelder and Horvath, 2008) analyses of n=7 selected "major cell types" consisting of cell populations that were annotated as "cell types" (B cells, Plasmablasts, NK cells) or "major subsets" (cMono, ncMono, CD4 T cells, CD8 T cells). For this analysis perpatient pseudobulk-summarized RPM-normalized counts from the prioritized sample set were used as input. Patients with low total cell numbers (<500 cells; n=3), with confirmed or suspected malignancies (n=3), or from the London COVID-19 cohort (n=2) were excluded from the analysis.

Input gene expression matrices were filtered to retain genes with expression above a certain threshold in a given minimum number of samples. Based on inspection of expression distributions RPM thresholds of 1 (cMono, NK cells, CD4 T cells and CD8 T) and 3 (ncMono, PB and B cells) were chosen for the indicated "major cell types". The number samples that constituted the smallest patient group was used as the minimum number of samples. After filtering we retained 9,000-12,000 genes per "major cell type". Log2 RPM+1 values were batch corrected using the ComBat algorithm (sva R package v 3.36.0) (Johnson et al., 2007; Leek et al., 2012) specifying the multiplexing sample pool as the adjustment variable (together with an intercept term).

The WGCNA analysis was performed using pipeline_wgcna.py from the https://github.com/sansomlab/cornet repository, and geneset over-representation analysis was carried out using the gsfisher R library https://github.com/sansomlab/gsfisher. The parameters used for WGCNA runs are given in Data S6. We excluded 3 "ambient_rna" modules whose eigengene gene loadings correlated strongly with ambient RNA species abundance (as quantitated in "empty" droplets) (Spearman's rho 0.83, p < 2.2 x 10 -16 ; data not shown) and unassigned grey modules from downstream analysis.

Modules were characterized and named by inspection of their gene membership, overrepresentation of biological pathways (GP Biological Process, Go Cellular Component, KEGG, MSigDB REACTOME and MSigDB transcription factor motifs) and correlation with AUCell (Aibar et al., 2017) (v 1.12.0) expression scores for specific sets of genes (Data S6) (identified WGCNA modules). The "curated type I IFN response" geneset comprised of GBP2 , IFI27, IFI44L, IFIT1,  IFITM1, IFITM3, IFNA1, IFNA10, IFNA13, IFNA14, IFNA16, IFNA17, IFNA2, IFNA21, IFNA4,  IFNA5, IFNA6, IFNA7, IFNA8, IFNB1, IRF1, IRF3, IRF7, IRF9, ISG15, MX1 MX2, OASL, RSAD2 , SIGLEC1, TRIM56, USP18 and XAF1.The "curated AP1 TF family" geneset comprised of the genes shown in Figure 3C . The custom "zf-C2H2 (PF00096) domain" geneset comprised of genes containing the PF00096 domain (Pfam database; Figure 3B ). For geneset over-representation analysis only genesets with a BH corrected p-value < 0.05 are reported (separate corrections performed within major cell type and ontology).

Module eigengene correlations with disease were computed by numerically coding the given x_vs_y disease groups as x=1 and y=0 respectively ( Figure 3A) . Module eigengene associations with variance between patient groups were assessed by ANOVA (COVID+HV_ANOVA, All_Groups_ANOVOA, Figure 3A) . Module eigengene correlations with clinical variables were computed using pairwise-complete observations. Where appropriate clinical variables were quantile normalized as described below for analysis of clinical phenotype. Stars in Figure 3A indicate a significant association between the contrast, clinical variable or geneset score and the module eigengene (BH corrected p-value < 0.05, p-value correction performed separately for each contrast, variable or score).

Here, we used both single cell (sc) and bulk cell VDJ sequencing to decipher the immune responses in COVID-19 patients. The benefit of scVDJ-seq is the joint analysis of the VDJ BCR or TCR sequence along with the gene expression and CITE-seq for each individual cell, allowing for a detailed linking of antigen receptor function or type with cellular phenotype or function. However, this is limited by the number of cells that can be captured in a single experiment, usually between a few 10's to 1,000's. Therefore, bulk VDJ-seq can also be used to capture a more holistic view of the immune cell repertoire. While this does not capture VDJ sequences on a single cell level, many aspects of the repertoire may be investigated, such as changes in clonality, dynamics, selection and tolerance. Because bulk VDJ-seq is able to capture the immune receptors of 1,000's to 100,000's of cells, this means that the resulting repertoire is more representative of the individual that single cell VDJ-seq, and has greater power to detect low frequency clones or repertoire features. Finally, through capturing a higher magnitude of BCR/TCR sequencing, imputing the germline IGHV/J germline alleles becomes possible, and therefore allowing for analyses of differential genetics between groups.

Additional pre-processing of scBCR-seq data for repertoire analysis The BCR outputs from CellRanger were run through IMGT V-QUEST (Giudicelli et al., 2004) to determine the TCR chain types and annotation. In addition, only BCRs from droplets that were confidently annotated as (a) singlets, (b) within the B cell annotated clusters in the GEX data, and (c) were confidently genetically demultiplexed were included. Processed FASTA sequences, corresponding to high confidence VDJ contigs from 10X Genomic's Cell Ranger v4.0.0 pipeline, were annotated using IMGT HIGHV-QUEST. To ensure high quality contigs, non-productive rearranged sequences were removed and only sequences corresponding to cell barcodes that past QC cut-offs for gene expression data were analyzed. For B cell VDJ sequences where multichains were detected, e.g. two or more heavy or light chain sequences, only contigs ranked above the first derivative of log10 ranked contig ratios were selected for downstream analysis. Contig ratios were defined by:

Where multi-chains passed this quality threshold, the contig corresponding to the lowest UMIs detected was considered ambient RNA contamination and removed. Where multi-chains did not pass this threshold, we considered these to be low confidence of being true single cell droplets and likely to represent homotypic doublets; thus, they were removed from further analysis.

Additional pre-processing of scTCR-seq data for repertoire analysis The TCR outputs from CellRanger were filtered based on called productivity, and chain identity (TRA or TRB). Only T cells that contained either (a) 1 beta chain, (b) 1 alpha and 1 beta, (c) 2 alpha chains and 1 beta were retained. Given that bona fide T cells with UMI counts of 1 have been validated in previous datasets (data not shown), the minimum number of UMIs required per cell to accept was 1. In addition, only TCRs from droplets that were confidently annotated as (a) singlets, (b) within the T cell annotated clusters in the GEX data, and (c) were confidently genetically demultiplexed were included.

The IMGT annotation of cell types demonstrated that the majority of droplets for which TCRs were captured were able to pass the chain type filters (Table S2) . After removal of low-quality droplets, there remain 94 samples corresponding to the maximal disease severity per patient. Four samples were excluded due to low T cell capture (<200 cells). Age was used as a co-variate in downstream analyses.

Clonal assignment was performed using a custom procedure. First, concatenated VH/VJ nucleotide sequences were clustered using an adapted R implementation of CD-HIT (https://github.com/thomasp85/FindMyFriends). Subsequently, within-cluster CDR3 normalized Levenshtein distances were generated using the "stringdist" R package. The clonal threshold was set as the local minimum of the CDR3 distance distribution. Convergent clones were assigned using the same procedure without constraints accounting for biological replicates. The presence of published SARS-CoV-2 binding antibody sequences from CoV-AbDab within the dataset was identified by first requiring identical VH and VJ gene pairings, followed by CDR3 distance thresholding (as described above).

The clonal expansion index (CEI) was calculated as the Gini index (unevenness) of the number of total BCRs per clone. The clonal diversification index (CDI) was calculated as the Gini index of the number of unique BCRs per clone. This is a measure of unevenness based on how many nonidentical members of each clone are diversified from their inferred germline ancestor.

A clone ID was defined by a concatenation of amino acid CDR3 chains present, with cells sharing identical clone IDs classed as members of the same clonotype. Clonal proportions were calculated by dividing the number of cells in each clone per sample by the total of number of cells per sample. Shannon diversity was calculated from count data of clones per sample using the R package entropy (Hausser, 2009 ). Statistical analysis was performed using a linear model with covariates for age and sample size.

For the scTCR analysis the cluster annotations were merged to create simplified clusters for analysis. The following reannotated clusters were used in TCR clonality analysis, with constituent "minor subsets" shown in brackets: CD8 T effector memory (CD8.TEMRA + CD8.TEM); CD8 T central memory (CD8.TCM + CD8.TCM.CCL5); CD8 T effector (CD8.TEFF + CD8.TEFF.prolif); CD8 naïve (CD8.NAIVE); CD4 T effector (CD4.TEFF + CD4.TEFF.prolif); MAIT (MAIT). Where comparisons of clone size across samples has been performed in terms of absolute number of member cells each sample was randomly down-sampled without replacement (n=1250 for CD4, n=250 for CD8) to generate repertoires of identical size across all samples. The mean size of clone was then calculated from the down-sampled repertoires. This process was bootstrapped 100 times and the mean size from these iterations is presented. This approach allows for accurate comparisons between samples containing different cell numbers and controls for sample size.

Cytotoxicity score was calculated using the AddModuleScore function in Seurat (Butler et al., 2018) with a gene-set of the top 50 genes that significantly correlated with IFNG expression in an independent single cell dataset (Watson et al., 2020a) and were identified as variable features in Seurat. The gene-set was then used as an input to generate phenotype scores.

CDR3b sequences of T cells identified as CD4+ or CD8+ T cells were broken down into 4-amino acid length sequences (Kmer). A Fisher's exact test was performed per Kmer across each group of patients (HV, COVID-19 or sepsis patients) and p values adjusted using a Bonferroni's correction, with a leave-one-out re-sampling of individuals employed to ensure inter-individual reliability. Kmer sequences significantly enriched in COVID-19 patients over HV or sepsis patients across more than 95% of the re-samples were identified as COVID-enriched Kmers. The proportion and cellular phenotype (i.e. subset and cytotoxicity) of T cells with a CDR3b sequence containing a COVID-enriched Kmer was subsequently compared between COMBAT clinical groups. The method used to identify CDR3 Kmers associated with COVID-19 is illustrated (Data S3).

Putative viral antigen specificity of clonotypes within COMBAT cohort was determined using publicly available databases of viral antigen specific CDR3 sequences. At the time, two databases containing TCR sequences with reported binding to SARS-CoV-2 epitopes were available: VDJdb (Bagaev et al., 2019) and ImmuneCODE™ (Nolan et al., 2020) . From VDJdb, all unique CDR3 amino acid sequences (alpha as well as beta chain of human TCR) of all lengths and MHCrestrictions with reported binding to known SARS-CoV2 associated antigens were selected. Sequences with low confidence scores assigned to TCR:peptide:MHC complexes or missing sequencing/ specificity validation data were then filtered out. ImmuneCODE™ database includes deeply sampled TCRb repertoires from over 1,400 subjects exposed to or infected with the SARS-CoV-2 virus. Within these, over 135,000 TCRs were deemed to be SARS-CoV-2-specific with high-confidence using Multiplex Identification of Antigen-Specific T-Cell Receptors Assay (MIRA) and were included (Klinger et al., 2015) . Unique CDR3 amino acid sequences thus identified were collated and interrogated against TCR sequences within COMBAT cohort for matches within CDR3a, CDR3a2 and CDR3b regions. A similar strategy was used to identify CMV, EBV and Influenza specific clones. Finally, proportions of individual repertoires occupied by viral clones were compared between COMBAT clinical groups.

Raw sequencing reads were filtered for base quality (median Phred score >32) using QUASR (Watson et al., 2013) . Forward and reverse reads were merged if they contained identical overlapping region of >50bp, or otherwise discarded. Universal barcoded regions were identified in reads and orientated to read from V-primer to constant region primer. The barcoded region within each primer was identified and checked for conserved bases. Primers and constant regions were trimmed from each sequence, and sequences were retained only if there was >80% per base sequence similarity between all sequences obtained with the same barcode, otherwise discarded. The constant region allele with highest sequence similarity was identified by 10-mer matching to the reference constant region genes from the IMGT database (Lefranc, 2011) , and sequences were trimmed to give only the region of the sequence corresponding to the variable (VDJ) regions.

Isotype usage information for each BCR was retained throughout the analysis hereafter. Sequences without complete reading frames and non-immunoglobulin/TCR sequences were removed and only reads with significant similarity to reference IGHTCR V and J genes from the IMGT database using BLAST (Altschul et al., 1990) were retained. Ig/TCR gene usages and sequence annotation were performed in IMGT V-QUEST, where repertoire differences were performed by custom scripts in Python.

Bulk BCR sequencing resulted in 2,356,813 BCRs consisting of 1,905,867 unique BCR sequences, resulting in 79 samples passing all QC measures and with the minimum number of BCR sequences per sample at 1000. Bulk TCR sequencing resulted in 1,200,656 TCRs consisting of 1,159,363 unique TCR sequences, resulting in 77 samples passing all QC measures and with the minimum number of TCR sequences per sample at 1,000.

Analysis methods are based on . To account for the greater numbers of BCR RNA molecules per plasmablast compared to other B cell subsets, the normalized isotype usages, defined as the percentage unique VDJ sequences per isotype, thus controlling for differential RNA per cell and reducing potential biases from differential RNA per cell. Similarly, mean somatic hypermutation levels and CDR3 lengths were calculated per unique VDJ region sequence to reduce potential biases from differential RNA per cell. IGHV gene usages were determined using IMGT, and proportions were calculated per unique VDJ region sequence. The representation of IGHV genes in the BCR repertoire reflects their presence in the germline, the naïve repertoire and their expansion after antigenic exposure. We therefore compared the frequency of IGHV gene use in PBMC-derived BCRs identified by sequence as being enriched for naive (IgM + D + SHM -: >78% naïve B cells by flow cytometry) and antigen-experienced B cells (including both unswitched (IgM + D + SHM + ) and class-switched memory (IgA + /G + /E + ) subsets) as shown in (Bashford-Rogers et al., 2019) .

Relative class-switch event frequency was the frequency of unique VDJ regions expressed as two isotypes (i.e. from more than one B cell, where one has undergone class-switch recombination). This was determined as proportion of unique BCRs present as both isotypes IgX and IgY within a random subsample of 8000 BCRs, where the mean of 1000 repeats was generated. This provides information on the frequency of BCRs observed associated with any two isotypes (class-switching events) while accounting for total read depth, but not accounting for differences in the relative frequencies of BCRs per isotype.

The per-isotype normalized class-switch event frequencies determine frequency of unique VDJ regions expressed as two isotypes whilst normalizing for differences in isotype frequencies.

To account for differences in isotype proportions, BCRs from each isotype were randomly subsampled to a fixed depth of 100 BCRs, and the proportion of unique VDJ sequences present between each pair of isotypes was counted. The mean of 1,000 repeats was generated.

The single-cell RNA-seq data sets were subjected to the standard RNA-velocity pipeline, and the trailing analyses were performed using the scVelo package (v0.1.24). The scVelo package was used to normalize the counts and select highly variable genes based on spliced counts. Following this, the dynamical model implemented in scVelo was used to estimate the RNA velocity for the cells. The estimated velocities were then visualized using Partition-based graph abstraction (PAGA) plots using the previously computed UMAP embeddings.

Supporting plots are provided in Data S3 for B cell cluster proportions; plasmablast diversity indices, clonal expansion index and clonal diversification index; correlation analysis of IGHV genes by single cell VDI and RNA VDJ; and IGHV gene usage of IGHD/M unmutated sequences, IGHD/M mutated sequences and class-switched sequences; constant region genes across B cell clusters; clonal overlap across comparator groups; and frequencies of known SARS-CoV-2 binding BCRs.

Raw data pre-processing was performed with Cell Ranger ATAC (10X Genomics). 'cellrangeratac count' pipeline was used to align reads and generate single-cell accessibility counts for the cells. The reference genome assembly file was Ensembl GRCh38 v100 Primary Assembly corresponding to hg38, downloaded from http://ftp.ensembl.org/pub/release-100/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz. This file was used as the reference genome file for alignment and generation of single-cell accessibility counts. Annotations were from Ensembl GRCh38 v100 Gene Annotation downloaded from http://ftp.ensembl.org/pub/release-100/gtf/homo_sapiens/Homo_sapiens.GRCh38.100.gtf.gz. VIREO (Huang et al., 2019) and demuxlet (Kang et al., 2018) packages were used to demultiplex patient samples within a channel and identify inter-sample nuclei doublets. Concordantly filtered barcodes from individualized samples were then used to create individual fragment files for each patient using HTSlib -c function. Following the creation of fragment files, downstream analysis of the scATAC-seq data was performed using the ArchR v0.9.3 R package (Granja et al., 2021) . Fragment files were first checked using the reformatFragmentFiles function. Arrow files were then generated by reading each sample's specific fragment files and a tile matrix was created using 500-bp bins. Cells with a transcription start site enrichment score < 4, with fewer than 1,000 detected fragments or containing intra-sample nuclei doublets were removed, resulting in ~46,000 cells, which were subjected to dimensionality reduction with iterative Latent Semantic Indexing (LSI) and Singular Value Decomposition (SVD), followed by Uniform Manifold Approximation and Projection (UMAP) embedding calculation to visualize the data structure in the two-dimensions. ~4,000 cells were manually removed at this stage, as we identified putative batch effects. The resulting ~42,000 cells were subjected to UMAP and were clustered using the implementation from Seurat R package (Stuart et al., 2019) . Cluster-specific gene activity scores were identified based on the local chromatin state, and marker genes were identified (FDR ≤ 0.01 & log2FC ≥ 1.25). Unconstrained integration with cognate scRNA-seq profiles was performed using the addGeneIntegrationMatrix (ArchR) method and scRNA-seq cell type annotations were used to label scATAC-seq clusters. We performed peak calling using MACS2 with the addReproduciblePeakSet (ArchR) function using pseudo-bulk replicates. Such replicates were grouped with different variables, such as cell-type, condition and patient, as well as a combination of these variables. Differentially accessible peaks (FDR <= 0.1 & log2FC >= 0.5) were identified between pairwise comparisons, and peak-to-gene linkages were calculated using the addPeak2GeneLinks (ArchR) method using a correlation cut-off of 0.45 and resolution = 1,000. We then used the 'cisbp' motif set to annotate motifs in accessible peaks using the addMotifAnnotations (ArchR) function. Motif enrichments in differentially accessible peaks were calculated using the peakAnnoEnrichment (ArchR) method. Finally, motif footprinting was performed by measuring Tn5 insertions in genome-wide motifs and normalized by subtracting the Tn5 bias from the footprinting signal. The concentrations of 51 circulating proteins in plasma were presented as mean fluorescence intensity (FI)  SEM and compared between HV and each disease severity group of COVID-19 using unpaired student's t-test, and then depicted with Prism (Data S3).

PCA on 171 plasma samples was conducted with the scater package in R. The values of FI were normalized with logNormCounts function, and then calculated with the RunPCA function. The loadings were generated by the value of first two components of eigenvectors multiplied by 10.

The heatmap was colored by the log10 of the fold change in the natural log of the FI, normalized against the mean value of HV for the plasma cohort and HS for the serum cohort. When comparing mortality in severe and critical COVID-19 patients, the data is normalized to the survivor group. The color-scale is bounded at ±5-fold change (0.7 in log10), with an increased FI shaded red; decreased FI shaded green; unchanged FI shaded yellow (Data S3).

The p-values in the Volcano plots are calculated using a two-tailed two-sample unpaired t-test (ttest2, MATLAB). The t-test was taken for the natural log of the FI of the test and control conditions. The p-values are plotted against the log2 of the fold change in the natural log of the FI between the test and control conditions.

The 

The correlation coefficient, r, was calculated using a Pearson correlation coefficient (corcoeff, MATLAB). All quoted r values have an associated p-value (also computed by corcoeff) of less than 0.05.

All values used for the correlation network analysis were the FI of Luminex results and/or the clinical readout. Severities were scored as HV=0, CH=1, CM=2, CS=3 and CC=4. Quality control (QC) was conducted on the matrix of expression values. Seven out of total 51 proteins were removed from the correlation analysis due to their low FI (<10 after subtracting FI in blank) and small SD (<10). Pearson correlation coefficient was performed using pairwise-complete correlation. Correlation matrix plots were generated using a modified package corrplot. Correlation matrix summary plots were made manually by R. Network plots were created by ggraph, and nodes were selected based on |r|>0.5 and p<0.05.

Data was analyzed by the Fragpipe pipeline consisting of Fragpipe 13.0 (Kong et al., 2017) , MSFragger 3.0, Philosopher 3.2.9 (da Veiga Leprevost et al., 2020) and Python 3.8.2. Blank runs were excluded and each file defined as experiment to facilitate LFQ. Data were searched against a fused target/decoy database generated by Philosopher and consisting of human UniProt SwissProt sequences and UniProt SARS-nCov02 (retrieved 17/07/2020), plus common contaminants. The database had 40,860 entries (including 50% decoy entries). MSFragger parameters were set to allow a precursor mass tolerance of plus/minus 10ppm and a fragment tolerance of 20ppm. Isotope error was left at 0/1/2 and masses were set to re-calibrate. Protein digestion was set to semi specific trypsin with up to 2 allowed missed cleavage sites, allowing peptides between 7 and 50 residues and mass range 500 to 5,000 Da. N-terminal protein acetylation and Methionine oxidation were set as variable modifications. ID validation was done with PeptideProphet and ProteinProphet (Nesvizhskii et al., 2003) with default settings.

Label free quantitation was conducted with IonQuant (Yu et al., 2020) and Match-Between-Runs enabled (with default parameters) and using Top-3 quantitation. Feature detection tolerance was set to 10ppm and RT Window to 0.6 minutes with an IM Window of 0.05 1/k0. For matching, ion, peptide and protein FDRs were relaxed to 0.1 and min correlation set to 0 in order to allow prefractionated library samples to be included. MBR top runs was set to 600.

The proteomics dataset was processed as follows: (1) Protein filtering such that proteins with at least 50% of valid values in one group were kept; (2) Sample filtering such that samples with more than 50% of missing values were removed from the dataset; (3) Data normalization with log2 transformation and median-centering of the dataset. Imputation of missing values was performed using a mixed model that combines a K-Nearest Neighbour approach (KNN) when at least 60% of valid values are present, otherwise a Minimum probability approach is used where missing values are randomly drawn from a Gaussian distribution (shift=1.8, nstd=0.3). The resulting data matrix contains 353 samples and 105 proteins. Thirteen samples were further excluded from analysis for malignancy, immunosuppression, or being alternative samples.

Unsupervised hierarchical clustering was performed based on Euclidean distance and Ward's method for calculating linkage (Data S3). Differential abundance analysis was performed by fitting protein abundance in linear models with the limma package, using only one sample at the maximal severity of the patient and including age and sex as covariates. The Benjamini-Hochberg procedure was applied to correct for multiple comparisons. FDR < 0.05 and fold change > 1.5 was taken as statistical significance. Pathway enrichment analysis was performed using the XGR package with annotations either from Gene Ontology Biological Process or the Reactome pathway database. Significantly enriched terms were defined by FDR < 0.05 in hypergeometric tests with all proteins detected in plasma (including in library samples) as the background. Statistical analysis was performed in R.

Protein-protein interaction data was retrieved from STRING v11 database with a confidence score cut-off of 0.7 and zero additional interactors. The network was visualized through Cytoscape v3.8.0 platform (Shannon et al., 2003) using perfuse force directed layout, and divided into clusters with the Markov cluster algorithm applied in the "clusterMaker" plugin. Node color was mapped to Pearson's correlation coefficients of PC1 scores and the protein level across samples, as lower PC1 score was shown to correlate with higher disease severity.

Page 52 of 92

The MS-based proteomics data were also analyzed using the Clinical Knowledge Graph (CKG) (Santos et al., 2020) . CKG provides a Python framework for downstream analysis and visualization of proteomics data: protein ranking, dimensionality reduction, functional Principal Component Analysis (PCA), Analysis of Variance (ANOVA), protein-clinical variable correlation analysis and network summarization.

CKG runs 3 feature reduction algorithms: Principal Component Analysis (PCA), Uniform Manifold Approximation and Projection (UMAP) and functional PCA. The functional PCA is based on the results of the method single-sample Gene Set Enrichment Analysis (Subramanian et al., 2005) , which identifies enrichment of Biological processes (Gene Ontology) (GOBP) in single samples derived from the ranked intensities of the identified proteins. This method generates a vector of biological processes enrichment scores for each sample. Loadings of the top 15 proteins and GOBP driving the separation of the conditions studied are included in the PCA and functional PCA respectively. In this analysis, the drivers are biological processes such as acute-phase response and inflammation and retinoid and lipoprotein metabolic processes and cholesterol transport.

We performed ANOVA analysis to identify differentially regulated proteins across conditions. Further, we run posthoc analysis (pairwise t-test) to show specific differences when comparing disease conditions to healthy volunteers or community COVID-19 and also between severity levels.

We performed functional enrichment analysis (Fisher's exact test) to identify enriched GO among the up-and down-regulated proteins in each pairwise comparison. Enrichments are plotted as scatterplots showing up-or down-regulated enriched GOBP and their adjusted p-values (Benjamini-Hochberg FDR: cutoff 0.05).

CKG performed a Spearman correlation analysis using the clinical metadata and the proteomics dataset. This Protein-clinical correlation is shown in a network where nodes are either clinical variables (diamond shape) or proteins (circle shape), and edges represent correlation (red=positive correlation, blue=negative correlation; Spearman correlation cutoff>=0.5; Benjamini-Hochberg FDR: cutoff 0.05). CKG applied a clustering algorithm (Louvain community detection) to identify clusters of highly connected nodes (nodes colored by cluster).

The results from all CKG analyses are summarized in a single visualization all the findings in the different analyses and all relevant biomedical context associated (diseases, drugs, biological processes, pathways, protein complexes, publications). The summarization algorithm prioritizes what nodes are shown in the network based on betweenness centrality. The top 15 central nodes are shown for each node type.

Plasma protein abundance for specific proteins is summarized (Data S3).

The similarity network fusion (SNF) analysis was performed using the function in CKG that makes use of the python library pySNF. The method was used to analyse the proteomics datasets (MSbased proteomics and Luminex) in combination to identify in an unsupervised manner clusters of similar patients. The number of clusters was not defined initially and optimized using the eigengap method and the clusters identified using Spectral clustering. The SNF analysis used Euclidean distance to calculate similarity (k_affinity=20, mu_affinity=0.6). The function returns the clusters and a mutual information score for each feature included in the analysis (MIscore). We used this score to prioritize a reduced number of features mainly driving the separation of the clusters (11 features) (MIscore>=0.15). The clusters are visualized using PCA plots.

In order to validate the identified clusters using SNF on the proteomics data, we used an independent cohort studied using a different technology, namely targeted proteomics by Olink (Filbin et al., 2020) . Access to the dataset was granted by Olink https://info.olink.com/mgh-covidstudy-overview-page. The processing of the data was done using CKG analytics core functions to map protein identifiers to names, transform the data into wide format and impute missing values using a mixed model as previously described. The clustering on these data followed a similar approach to how the SNF clusters were calculated (optimal clusters and Spectral clustering) but using only the selected features in the SNF analysis (7/11 features).

To evaluate the clinical relevance of the identified clusters in COMBAT, we performed a survival analysis and plotted the Kaplan-Meier curve using R packages survival (Therneau and Grambsch, 2000) and survminer. The input data is a data frame specifying the time to event, the event (death or end of observation) and the groups (SNF clusters). The comparison of the survival distributions between clusters was performed and the p value given using log-rank test. The hazard ratio was calculated using Cox proportional hazard model. In the Olink dataset survival status is only available at 4 timepoints: 0, 3, 7 and 28 days. Deaths at these time points were collected according to WHO category 1, defined as death, in these timepoints (WHO 0, WHO 3, WHO 7 and WHO 28). We compared the 28-day mortality between the two SNF clusters by chi-squared test.

To evaluate the robustness of the identified clusters between technologies/studies and to eliminate the possibility that only the 7 selected features used for the clustering of the Olink data were similar, we calculated the correlation (Pearson correlation) of fold-changes between clusters for all the common proteins in these studies (n=43).

A tensor and matrix decomposition method, Sparse Decomposition of Arrays (SDA), as defined in (Hore et al., 2016) , was used to integrate 152 samples from the different 'omics datasets defined above, allowing that some samples were missing certain 'omics data.

The whole blood total RNA-seq (9 missing samples) and the pseudobulk from 10X CITE-seq scRNA-seq (22 missing samples) were combined into a three-dimensional tensor consisting of 152 samples by 9 tissue types by 14,989 genes (which passed QC in both datasets). This expression was normalized by sample in each tissue type by log2-transformation of counts-permillion + 1.

The number of cells per cell type was included as defined by 10X CITE-seq (a two-dimensional matrix of 152 samples by 97 cell types, with 22 samples missing); and mass cytometry (CyTOF) (One two-dimensional matrix of 152 samples by 10 cell types for the all cells dataset, with 21 samples missing, and a further two-dimensional matrix of 152 samples by 51 cell types for the granulocytes-depleted cells dataset, with 20 samples missing). We filtered out any samples with fewer than 500 cells in any matrix. The data in each matrix was normalized by a log2 transformation of counts-per-million + 1.

The proteomics data from Luminex (in a two-dimensional matrix of 152 samples by 51 proteins, with 20 samples missing) and mass spectrometry (Tims-TOF) (in a two-dimensional matrix of 152 samples by 105 proteins, with 17 samples missing) were used, with the data normalized as described in the Luminex and Tims-TOF sections.

As described in (Hore et al., 2016) , to find robust components we ran the tensor and matrix decomposition ten times for 1,000 components. Each time, around 290 components were estimated to be zero. Once again, similar to the (Hore et al., 2016) method the absolute correlation (r) was calculated for the sample scores for each pair of components, clustered using hierarchical clustering on 1-r (dissimilarity measure) and formed flat clusters in which the components in each flat cluster have no greater a cophenetic distance than 0.4. We chose the flat clusters that had components from at least 5 of the 10 runs. The final sample, tissue and gene or protein or cell score was the mean of all the components within the chosen clusters. This resulted in 381 clusters.

Components were identified as associated with COVID-19 if they (1) showed significant variation (BH adjusted p < 0.01) in an analysis of variance between the COVID-19 categories and healthy volunteers and (2) showed a significant |spearman's rho|>= 0.5 (and Benjamini/Hochberg adjusted p < 0.01) in at least one of the contrasts between the COVID-19 groups vs healthy volunteers (in total we found n=130 such components) ( Figure S11C . To identify COVID-19 specific components the median loadings of the components for the different comparator (source) groups (i.e. including influenza and all-cause sepsis) were clustered (k-means). Individual component associations, for example with comparator group, severity or clinical features were assessed with Spearman correlation (between component loadings and numerical variables) or ANOVA (between component loadings and categorical variables). The overview heatmap ( Figure  S11A ) was generated by combining the top significant (Benjamini/Hochberg adjusted p < 0.05, abs(r) > 0.5, max 10) components from each of the individual analyses. Pathway enrichment of gene expression (those with posterior inclusion probability > 0.5, weighted by their ranking in loading score magnitude) for individual components was done using gene set enrichment analysis as implemented in Pi's xPierGSEA (Fang et al., 2019) .

We used supervised machine learning (from sklearn) to classify samples according to their WHO severity based on PCs obtained from data across the different modalities (timsTOF, Luminex, and total RNA-seq) ( Figure S10A) . We performed permutation feature scoring to find the most important PCs to predict severity. After that, we extracted the most important features of the most important PCs and reran the algorithm directly on these features, again ranking them according to their importance.

The pre-processed data from individual COMBAT assay datasets were automatically processed using the standard extract-transform-load (ETL) procedure to generate an integrated dataset. Datasets were merged using shared variables, the COMBAT sample ID and assay-specific sample IDs. Next, we generated novel features, such as name of specific assays to indicate if a sample was analysed in a specific assay (yes/no), ordering of samples shared between assays (first_sample_across_assays), day when sample was obtained from the maximum disease (day_sampling_from_max_disease), state of the disease when sample was acquired (recovered/ongoing), hospitalized, ventilation, oxygenation, sampling (<=10 days as early or >10 days as late phase of the disease), samples acquired before or after maximum disease (sampling_from_max_disease) and disease progress for longitudinal samples (deterioration or recovery) . Donors that were excluded or frail were not included in the integrated dataset. We standardized names of the immune cells subsets and analytes to reflect the measurement, such as frequency (freq_cell subsets) or intensity (Luminex parameter_intens). In total, the integrated COMBAT dataset contained information on 428 samples from 268 donors on more than a million parameters.

For the multi-omics data integration, after filtering for samples not analyzed across all assay modalities and restricting to the first available sample after admission from sepsis and hospitalized COVID-19 patients, the final dataset included 15 sepsis patients and 53 hospitalized COVID-19 patients with 184 features analyzed using CyTOF, 79 using flow cytometry, 8 using GSA, 105 mass spectrometry, 102 Luminex and 23,063 features from whole blood total RNA-seq. The data for each assay was centered and scaled, missing values were imputed, features with zerovariance and near-zero-variance were removed and finally, highly correlated features with cut-off 0.85 were also removed.

To avoid 'curse of dimensionality', reduce overfitting and improve accuracy, we implemented a wrapper approach as described (Kohavi and John, 1997) . Briefly, the initial dataset containing all 68 donors was partitioned into training (52 donors) and test set (16 donors) with balanced class distribution of sepsis and COVID-19 patients using the data partitioning function (Kuhn, 2008) as described previously (Tomic et al., 2019) . The same training/testing dataset split was used for each assay. The reduction of features using the recursive feature elimination (RFE) algorithm was performed on the training set and the model was evaluated using the 10-fold cross-validation repeated 5 times. For whole blood total RNA-seq, prior to training the model using RFE algorithm, we have performed analysis of differentially expressed genes between sepsis and hospitalized COVID-19 patients, reducing the number of genes to 2,989 based on the FDR <0.05 and fold change >1.5, and accounted for age and sex in the model. The RFE method removes the features that do not contribute to the final model, while it keeps the features that contribute the most to the final model as evaluated using the variable importance score (Kuhn and Johnson, 2019) . Finally, after the RFE analysis, the selected features from each assay (36 features from CyTOF, 20 from flow cytometry, 20 from Luminex, 32 from mass spectroscopy, 28 from whole blood total RNA-seq and 1 from GSA dataset) were merged in the final training dataset containing 52 donors and 137 features.

To identify immunological and molecular signatures that can discriminate between sepsis and COVID-19 patients, we used SIMON (Sequential Iterative Modeling "Over Night") (Tomic et al., 2019) . SIMON is a free and open-source software that provides a standardized ML method for data pre-processing, data partitioning, building predictive models, evaluation of model performance and selection of features. For the analysis, we applied four machine learning algorithms, naïve Bayes, Shrinkage discriminant analysis (sda), Support vector machine with linear kernel ('svmLinear2') and C5.0 decision tree. Since the entire ML process in SIMON is unified, resulting models built with different algorithms can be compared and the best performing models can be selected. First, models are built on training set and the performance is evaluated using a 10-fold cross-validation repeated five times and cumulative error rate is calculated. To prevent overfitting, in the second step, each model is evaluated on the withheld test set. The performance of classification models was determined by calculating the area under the receiver operating characteristic curve (AUROC) for test set (test AUROC). The best performing model was built using the naïve Bayes (testing AUROC 0.85, 95% CI 0.59-1.00). In the final step, SIMON calculated the contribution of each feature to the model as variable importance score (scaled to maximum value of 100).

To support consistent and coherent communication of data and metadata within the project, a unified identifier system for all samples was implemented. The COMBAT sample identifier system encodes information regarding the sample providence in terms of cohort, de-identified patient ID, location where the sample was taken, the time point of the sample relative to the initial collection and details regarding the processing of the sample itself.

Datasets from each modality were stored within the consortium via the COMBAT datawarehouse, consisting of over 100TB of fast storage connected to a research computing cluster. This enabled data processing to occur within the datawarehouse, reducing the risk of duplication of datasets and the possibility of uncontrolled changes. Once datasets were ready to be shared within the consortium, to support (for example) data integration work, they were formally given a unique identifier and placed in a dedicated dataset directory. The existence of this deposition and its associated metadata, including information regarding associated samples and status was then made available via a web application which captures this information in a back-end database. This application allows consortium members to search by modality and status, providing information about the purpose of each dataset and its location in the datawarehouse.

The governance of data management was supported by the existence of a short but well-defined Data Access Agreement, which all consortium members were required to sign before gaining access to the datawarehouse. Furthermore, granular permissions within the datawarehouse enabled careful access controls to be applied to particularly sensitive data (such as rich clinical data). Web applications supporting the consortium are all protected by a federated Shibbolethbased authentication approach, allowing collaborators from outside of Oxford to gain access as required.

In order to visualize the data from the different modalities (experiments), a module, Multi Experiment Viewer (MEV https://github.com/Hughes-Genome-Group/MEV) was built for Multi Locus View (MLV DOI 10.1101 such that the data was pivoted on sample rather than genome location. The data was loaded in from each modality and for the CITE-seq data, pseudobulk values (on a sample/cell type basis) were used. In order to compare radically different datatypes (read counts, fluorescence, % cell type etc.), percentiles for each observation were calculated from the 10th to the 90th, in steps of 10. Values were then placed into 10 bins, e.g., if a value was <=10th percentile it would be given a value of 1 and if >2nd and <=3rd, it would get a value of 2 etc. This was used as the default value, although depending on the modality other values e.g., raw counts, log transformed values were also added and can be selected by the user. Users can create their own views by searching for genes or loading in specific data sets and then combine them with a limited set of clinical data. Charts such as histograms, heatmaps and scatter plots can then be added and cross-filtered to identify samples or features of interest. An instance can be found at https://mlv.combat.ac.uk, with links to predefined views of the data.

Page 57 of 92 Shiny apps (https://shiny.combat.ox.ac.uk) were developed using the R package shiny to display the results of the whole blood total RNAseq differential expression analysis and the principal component analysis. In each, derived data are loaded (limma fitted models and pre-calculated principal components respectively) together with limited metadata, and plots are generated within the app using ggplot2. Figure  1 and Data S4 NLRP1  MAP3K14  OGFRL1  RPS6KA5  FOXP1  RTN1  HGSNAT  DNAJC10  NAP1L1  COQ8A  ZNF652  TRIO  OGT  ATG16L2  DPEP2  CMTM3  SNX29  TBC1D9  ZDHHC7  CSF3R  MS4A4A  AP1S1  NKG7  IDH1  RNASE2  NDUFB3  PSMD1  H2AJ  CCND3  MCEMP1  UBE2F  RPN2  LMAN2  S100A8  COPG1  CYP19A1  FDPS  SLC39A8  SQOR 

Spearman correlation (between eigengene and dummy-coded contrast or continuous variable/score) or ANOVA (variance of eigengene between groups; scaled F-value 0-1). *within-column BH adjusted P < 0.05. Highlights  Blood atlas delineating innate and adaptive immune dysregulation in COVID-19  Shared and specific immune signatures of COVID-19, influenza and all cause sepsis  Multi-omic immune profiling differentiates hospitalised patient severity in COVID-19  Immune activation and proliferation involving AP-1/p38MAPK associated with COVID-19

In Brief A multi-omic analysis of patient blood samples reveals both similarities and specific features of COVID-19 when compared with samples obtained from sepsis or influenza patients which could yield better targetted therapies for severe COVID-19.

SCENIC: single-cell regulatory network inference and clustering

Basic local alignment search tool

Orchestrating single-cell analysis with Bioconductor

FastQC: A Quality Control Tool for High Throughput Sequence Data

xCell: digitally portraying the tissue cellular heterogeneity landscape

Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage

Systems biological assessment of immunity to mild versus severe COVID-19 infection in humans

A Novel LC System Embeds Analytes in Preformed Gradients for Rapid, Ultra-robust Proteomics

VDJdb in 2019: database extension, new analysis infrastructure and a T-cell receptor motif compendium

Analysis of the B cell receptor repertoire in six immune-mediated diseases

Sepsis and Coronavirus Disease 2019: Common Features and Anti-Inflammatory Therapeutic Approaches

Generalizing RNA velocity to transient cell states through dynamical modeling

Longitudinal Multi-omics Analyses Identify Responses of Megakaryocytes, Erythroid Cells, and Plasmablasts as Hallmarks of Severe COVID-19

A prospective analysis of lymphocyte phenotype and function over the course of acute sepsis

Deciphering the state of immune silence in fatal COVID-19 patients

Potent neutralizing antibodies from COVID-19 patients define multiple targets of vulnerability

Integrating single-cell transcriptomic data across different conditions, technologies, and species

Topological methods for genomics: present and future directions

Gene-set integrative analysis of multi-omics data using tensorbased association test

T cell responses in patients with COVID-19

CGAT-core: a python framework for building scalable, reproducible computational biology workflows. F1000Research

Distinct fibroblast subsets drive inflammation and damage in arthritis

Philosopher: a versatile toolkit for shotgun proteomics data analysis

Immune therapy in sepsis: Are we ready to try again?

The Genomic and Immune Landscapes of

A complete tool set for molecular QTL discovery and analysis

RNA-SeQC: RNA-seq metrics for quality control and process optimization

Reduction and Functional Exhaustion of T Cells in Patients With Coronavirus Disease 2019 (COVID-19)

STAR: ultrafast universal RNA-seq aligner

Developmental Analysis of Bone Marrow Neutrophils Reveals Populations Specialized in Expansion, Trafficking, and Effector Functions

Peripheral CD8(+) T cell characteristics associated with durable responses to immune checkpoint blockade in patients with metastatic melanoma

Multi-insight visualization of multi-omics data via ensemble dimension reduction and tensor factorization

A genetics-led approach defines the drug target landscape of 30 immune-related traits

Pi: Leveraging genetic evidence to prioritise drug targets at the gene, pathway and network level

XGR software for enhanced interpretation of genomic summary data

ANOVA-like differential expression (ALDEx) analysis for mixed population RNA-Seq

Serial evaluation of the SOFA score to predict outcome in critically ill patients

Plasma proteomics reveals tissuespecific cell death and mediators of cell-cell interactions in severe COVID-19 patients

IMGT/V-QUEST, an integrated software program for immunoglobulin and T cell receptor V-J and V-D-J rearrangement analysis

Interleukin-6 Receptor Antagonists in Critically Ill Patients with Covid-19

Targeting Macrophages as a Therapeutic Option in Coronavirus Disease

ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis

p38 MAPK inhibition: A promising therapeutic approach for COVID-19

Algebraic Systems Biology: A Case Study for the Wnt Pathway

Impaired type I interferon activity and inflammatory responses in severe COVID-19 patients

Entropy Inference and the James-Stein Estimator, with Application to Nonlinear Gene Association Networks

Efficient integration of heterogeneous single-cell transcriptomes using Scanorama

Baricitinib treatment resolves lower-airway macrophage inflammation and neutrophil recruitment in SARS-CoV-2-infected rhesus macaques

The role of growth factor receptors in viral infections: An opportunity for drug repurposing against emerging viral diseases such as COVID-19?

Dexamethasone in Hospitalized Patients with Covid-19

Tocilizumab in patients admitted to hospital with COVID-19 (RECOVERY): preliminary results of a randomised, controlled, open-label, platform trial

Tensor decomposition for multiple-tissue gene expression experiments

Vireo: Bayesian demultiplexing of pooled single-cell RNA-seq data without genotype reference

Matplotlib: A 2D Graphics Environment

Adjusting batch effects in microarray expression data using empirical Bayes methods

Phenotypical and functional alteration of unconventional T cells in severe COVID-19 patients

Baricitinib plus Remdesivir for Hospitalized Adults with Covid-19

Multiplexed droplet single-cell RNAsequencing using natural genetic variation

Empirical statistical model to estimate the accuracy of peptide identifications made by MS/MS and database search

Multilayer networks

Multiplex Identification of Antigen-Specific T Cell Receptors Using a Combination of Immune Assays and Immune Receptor Sequencing

Wrappers for feature subset selection

MSFragger: ultrafast and comprehensive peptide identification in mass spectrometry-based proteomics

Fast gene set enrichment analysis

Fast, sensitive and accurate integration of singlecell data with Harmony

Building predictive models in R using the caret package

Feature engineering and selection. A practical approach for predictive models

Comprehensive mapping of immune perturbations associated with severe COVID-19

Combinatorial Single-Cell Analyses of Granulocyte-Monocyte Progenitor Heterogeneity Reveals an Early Uni-potent Neutrophil Progenitor

RNA velocity of single cells

WGCNA: an R package for weighted correlation network analysis

Immunophenotyping of COVID-19 and influenza highlights the role of type I interferons in development of severe COVID-19

The sva package for removing batch effects and other unwanted variation in high-throughput experiments

IMGT unique numbering for the variable (V), constant (C), and groove (G) domains of IG

Activation and evasion of type I interferon responses by SARS-CoV-2

Cumulus provides cloud-based data analysis for large-scale single-cell and single-nucleus RNA-seq

featureCounts: an efficient general purpose program for assigning sequence reads to genomic features

The Molecular Signatures Database (MSigDB) hallmark gene set collection

EmptyDrops: distinguishing cells from empty droplets in droplet-based single-cell RNA sequencing data

c-Jun overexpression in CAR T cells induces exhaustion resistance

Longitudinal immune profiling reveals key myeloid signatures associated with COVID-19

Mature CD10(+) and immature CD10(-) neutrophils present in G-CSF-treated donors display opposite effects on T cells

Why have clinical trials in sepsis failed?

AP-1 imprints a reversible transcriptional programme of senescent cells

Identification and characterization of VEGF-A-responsive neutrophils expressing CD49d, VEGFR1, and CXCR4 in mice and humans

Deep immune profiling of COVID-19 patients reveals distinct immunotypes with therapeutic implications

Scater: pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R

Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation

Immunoglobulin heavy chain expression shapes the B cell receptor repertoire in human B cell development

Consensus Clustering: A Resampling-Based Method for Class Discovery and Visualization of Gene Expression Microarray Data

Cytokine release syndrome in severe COVID-19

Distinct inflammatory profiles distinguish COVID-19 from influenza with limited contributions from cytokine storm

Normalizing and denoising protein expression data from droplet-based single cell profiling

fastcluster: Fast Hierarchical, Agglomerative Clustering Routines for R and Python

A statistical model for identifying proteins by tandem mass spectrometry

A large-scale database of T-cell receptor beta (TCRβ) sequences and binding associations from natural and synthetic exposure to SARS-CoV

Instant Clue: A Software Suite for Interactive Data Visualization and Analysis

CyTOF workflow: differential discovery in highthroughput high-dimensional cytometry datasets

Parallels in Sepsis and COVID-19 Conditions: Implications for Managing Severe COVID-19

Coronavirus Disease 2019-Associated Thrombosis and Coagulopathy: Review of the Pathophysiological Characteristics and Implications for Antithrombotic Management

MAIT cell activation and dynamics associated with COVID-19 disease severity

Nucleotide sequence analysis of the V regions of two IgM cold agglutinins. Evidence that the VH4-21 gene segment is responsible for the major cross-reactive idiotype

Scikit-learn: Machine Learning in Python

The PRIDE database and related tools and resources in 2019: improving support for quantification data

Immune checkpoint inhibition in COVID-19: risks and benefits

BBKNN: fast batch alignment of single cell transcriptomes

Heparin reduces nonspecific eosinophil staining artifacts in mass cytometry experiments

CoV-AbDab: the coronavirus antibody database

limma powers differential expression analyses for RNA-sequencing and microarray studies

edgeR: a Bioconductor package for differential expression analysis of digital gene expression data

Clinical Knowledge Graph Integrates Proteomics Data into Clinical Decision-Making. bioRxiv

A diverse range of gene products are effectors of the type I interferon antiviral response

Severe COVID-19 Is Marked by a Dysregulated Myeloid Cell Compartment

Tensor clustering with algebraic constraints gives interpretable groups of crosstalk mechanisms in breast cancer

Derivation, Validation, and Potential Treatment Implications of Novel Clinical Phenotypes for Sepsis

Cytoscape: a software environment for integrated models of biomolecular interaction networks

Use of machine learning to identify a T cell response to SARS-CoV-2

Elevated Calprotectin and Abnormal Myeloid Cell Subsets Discriminate Severe from Mild COVID-19

The cellular immune response to COVID-19 deciphered by single cell multi-omics across three UK centres

Comprehensive Integration of Single-Cell Data

Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles

One-class Differential Expression Analysis using Tensor Decompositionbased Unsupervised Feature Extraction Applied to Integrated Analysis of Multiple Omics Data from 26 Lung Adenocarcinoma Cell Lines

The trinity of COVID-19: immunity, inflammation and intervention

Modeling Survival Data: Extending the Cox Model

Automated Machine Learning System, Reveals Immune Signatures of Influenza Vaccine Responses

SIMON: Open-Source Knowledge Discovery Platform. Patterns (N Y)

Perseus: A Bioinformatics Platform for Integrative Analysis of Proteomics Data in Cancer Research

Genomics in the Clud: Using Docker, GATK, and WDL in

The stringdist Package for Approximate String Matching

CytoNorm: A Normalization Algorithm for Cytometry Data

Visual analysis of mass cytometry data by hierarchical stochastic neighbour embedding reveals rare cell types

Kepler Mapper: A flexible Python implementation of the Mapper algorithm

SciPy 1.0: fundamental algorithms for scientific computing in Python

Similarity network fusion for aggregating data types on a genomic scale

Diverse Functional Autoantibodies in Patients with COVID-19

Immune Checkpoint Blockade induces peripheral cytotoxicity and persistence of large effector CD8<sup>+</supT cell clones

Immune Checkpoint Blockade induces peripheral cytotoxicity and persistence of large effector CD8<sup>+</sup> T cell clones. bioRxiv

Viral population analysis and minority-variant detection using short read nextgeneration sequencing

diffcyt: Differential discovery in high-dimensional cytometry via high-resolution clustering

Factors associated with COVID-19-related death using OpenSAFELY

The single-cell epigenomic and transcriptional landscape of immunity to influenza vaccination

SCANPY: large-scale single-cell gene expression data analysis

Scrublet: Computational Identification of Cell Doublets in Single-Cell Transcriptomic Data

The functional paradox of CD43 in leukocyte recruitment: a study using CD43-deficient mice

Prediction models for diagnosis and prognosis of covid-19 infection: systematic review and critical appraisal

Single-cell epigenomic landscape of peripheral immune cells reveals establishment of trained immunity in individuals convalescing from COVID-19

Identification of immune correlates of fatal outcomes in critically ill COVID-19 patients

IonQuant enables accurate and sensitive label-free quantification with FDR-controlled match-between-runs

Fast Quantitative Analysis of timsTOF PASEF Data with MSFragger and IonQuant

NK.CD16hi.IFN.resp NK.CD16hi.cyc

Cat# 99814 anti-human CD273 (B7-DC, PD-L2) ( 

Contains two anndata objects:(1) "COMBAT-CITESeq-DATA": the raw and normalized gene expression and ADT data, cluster annotations, summary repertoire information and detailed metadata as reported and analyzed in this study.(2) "COMBAT-CITESeq-EXPRESSION-ATLAS": raw and normalized gene expression data from an alternative mapping of the data to a transcriptome index that additionally included genes belonging to the "lncRNA" biotype. With ADT data, cluster annotations and a subset of clinical metadata. Suitable for visualization with cellxgene. This paper 

Anonymized patient metadata, data dictionary. This paper EGAD00001007931 CBD-RAW-CYTOF-MYELOID Single cell resolution mass cytometry (CyTOF) data generated from whole blood and granulocyte depleted whole blood, with a focus on myeloid populations. Dataset contains ungated, raw data files in fcs format and the list of antibodies used. This paper 10.5281/zenodo.513 9561 CBD-RAW-CYTOF-WB Single cell resolution mass cytometry (CyTOF) data generated from whole blood. Dataset contains ungated, raw data files in fcs format and the list of antibodies used. ( Software and Algorithms ALDEx2 (Fernandes et al., 2013) v1.18.0ArchR (Granja et al., 2021) v0.9.3 https://www.archrp roject.com/ AUCell (Aibar et al., 2017) v1.12.0 BBKNN (Polanski et al., 2020) https://github.com/ Teichlab/bbknn BLAST (Altschul et al., 1990) v0.0.5Cytoscape (Shannon et al., 2003) v3.8.0Cytosplore (van Unen et al., 2017) www.cytosplore.or g demuxlet (Kang et al., 2018) v2 https://github.com/ statgen/demuxlet diffcyt v1.8.8 edgeR (Robinson et al., 2010) v3.28.1 https://bioconducto r.org/packages/rel ease/bioc/html/edg eR.html EmptyDroplets (Lun et al., 2019) v1.8.0 Entropy (Hausser, 2009) http://www.strimme rlab.org/software/e ntropy/ eXploring Genomic Relations (XGR) (Fang et al., 2016b) http://galahad.well. ox.ac.uk/XGR Fastcluster (Mulner, 2013) v1.1.25 FastQC (Andrews, 2010) v0.11.9 https://github.com/ s-andrews/FastQC featureCounts (Liao et al., 2014) v1.6.4 fgsea (Korotkevich et al., 2021) https://bioconducto r.org/packages/rel ease/bioc/html/fgs ea.html J o u r n a l P r e -p r o o f (Lefranc, 2011) http://www.imgt.or g/ IMGT V-QUEST (Giudicelli et al., 2004) http://www.imgt.or g/IMGTindex/V-QUEST.php InstantClue (Nolte et al., 2018) v0.5.3 http://www.instantc lue.uni-koeln.de/ IonQuant (Yu et al., 2020) https://ionquant.ne svilab.org/ KeplerMapper (van Veen et al., 2019) https://reposhub.co m/python/datavalidation/scikittda-keplermapper.html Limma (Ritchie et al., 2015) https://bioconducto r.org/packages/rel ease/bioc/html/lim ma.html MATLAB https://uk.mathwor ks.com/help/matla b/ Matplotlib (Hunter, 2007) https://matplotlib.or g/ MSFragger (Kong et al., 2017) v3.0 MSigDB (Subramanian et al., 2005)