Plasma cell-free RNA signatures of inflammatory syndromes in children

We collected and analyzed 370 plasma samples from pediatric patients with inflammatory and infectious conditions at four hospitals in the United States, including Rady Children’s Hospital San Diego (UCSD), Emory, Children’s National Hospital (CNH), and University of California San Francisco (UCSF) (Fig. 1A). This sample set included patients diagnosed with KD, MIS-C, viral infection, bacterial infection, and other hospitalized pediatric controls, some but not all with inflammatory disease (for example, arthralgia, Crohn’s disease flare, parenchymal lung disease, chronic lung disease, toxic shock syndrome, and postvaccine myocarditis), as well as healthy children (Fig. 1A and SI Appendix, Table S1). Included cases of bacterial and viral infections were heterogeneous with respect to infection site and pathogen (Dataset S1). Age, sex, race, severity, and inflammatory status were recorded for each patient (Fig. 1 B–E). All patient samples were collected within 4 d of hospitalization during the acute phase of disease (Materials and Methods). RNA next-generation sequencing was used to quantify the genes and cell types of origin of cfRNA in each sample, with an average of 20.7 million reads sequenced per sample (Materials and Methods). No batch effect was observed between samples from different hospitals (Materials and Methods and SI Appendix, Fig. S1 A and B). Machine learning models were built to identify biosignatures associated with different diseases.

Our initial focus was on characterizing changes in cfRNA profiles in plasma that are common across the various inflammatory conditions. We performed pairwise differential gene abundance analysis between the healthy control group and each disease group (Materials and Methods, Fig. 2A, and Dataset S2). This analysis identified differentially abundant genes (DAGs) for all comparisons [Benjamini–Hochberg (BH) adjusted P-value < 0.05]. The smallest number of DAGs (n = 2,686) was identified when comparing healthy controls to patients in the other hospitalized control group, likely due to the smaller sample size and higher heterogeneity of this sample group. The greatest number of DAGs (n = 6,606) was identified when comparing healthy controls to patients with KD (Fig. 2A). Further analysis of the DAGs for all conditions revealed a significant number of shared genes (1,871 DAGs, Fig. 2A). Chief among these were histone protein coding genes that were elevated for all disease groups, indicating that an increase in histone transcripts is a universal indicator of inflammation (SI Appendix, Fig. S2A). The analysis also identified immune related transcripts Myeloperoxidase (MPO), Elastase, Neutrophil Expressed (ELANE), CD53, and C-X-C Motif Chemokine Receptor 2 (CXCR2) elevated in each disease group, homeostasis related transcripts CDK19, ANAPC5, and 32 mitochondrial protein coding RNAs in the control group (SI Appendix, Fig. S2B). Pathway analysis of the 1,871 overlapping DAGs revealed an enrichment of neutrophil and cell replication transcripts related to inflammation (Fig. 2B and Dataset S3). These findings point to shared signatures of inflammation among disease groups and emphasize the need for intergroup comparisons to develop disease-specific biosignatures.

KD and MIS-C share many clinical characteristics: They are highly inflammatory, present with endothelial dysfunction, and have multiple overlapping signs, including skin rashes, mucosal involvement, and fever. There are currently no molecular tests to distinguish between KD and MIS-C. We therefore assessed whether cfRNA could be used to differentiate these severe pediatric inflammatory syndromes. We divided the KD (n = 101) and MIS-C (n = 97) samples into training (60%), validation (20%), and test (20%) sets, ensuring a roughly equal representation of hospital of origin and disease subclassification across all three sets (Fig. 3A). We used the training data for feature selection and to train machine learning models. We then used the validation set to select the final model, based on the model with the highest receiving operator characteristic area under the curve (AUC). Last, we evaluated the performance of the final model using the test set. To prevent influence of the test set on the training data and to ensure unbiased results, we normalized each set separately using a variance stabilizing transformation (Materials and Methods).

In our initial analysis using the training set, we identified 1,242 DAGs between KD and MIS-C (DESeq2, BH adjusted P-value < 0.01 and base mean > 10; Fig. 3B and Dataset S4). We refined this gene list based on adjusted P-value, base mean, individual gene AUC, and fold change, resulting in a final tally of 132 genes for model input (Fig. 3C). We then trained 14 machine learning classification models, evaluating their performance on the validation set (Materials and Methods). The generalized linear model with elastic net regularization (GLMNET) model with least absolute shrinkage and selection operator (LASSO) regression exhibited the highest validation set AUC and was selected as our final model (Fig. 3D). A unique feature of the GLMNET LASSO algorithm is the feature selection step, which selected 25 genes for the modeling (SI Appendix, Table S2). Using this trained model, we tested the classification performance using the test set [receiver operating characteristic (ROC)-AUC train = 1.00, validation = 0.98, and test = 0.98] (Fig. 3E). We also observed similar distributions in the classification scores across sample sets, further confirming that there was little to no overfitting of the training and test sets (Fig. 3F).

We compared the KD vs. MIS-C gene panel to previously reported differentially expressed genes (DEGs) from a study of cultured endothelial cells treated with KD and MIS-C sera (24). We found significant overlap in genes; nine of the 25 genes in our panel were also identified as DEGs in the previously reported study (nominal P-value < 0.05). Of the overlapping genes, those elevated in acute KD are associated with autophagy (ANXA6, CD44, and PDCD4) and endothelial–mesenchymal transition (EndoMT) (STK17B and VIM), processes linked to endothelial cell function in inflammatory conditions (25–29). This provides evidence that the key differences in endothelial dysfunction and cardiac outcomes between MIS-C and KD may be linked to differences in endothelial autophagy and activation of EndoMT pathways.

We next asked whether cfRNA could be used in the more challenging scenario of multidisease classification. For this, we developed a machine learning framework that combines the outputs of one-vs.-one models using a random forest (RF) multiclass algorithm (Materials and Methods and Fig. 4A). We first split our dataset into training (70%) and test (30%) sets, with roughly comparable proportions of samples from patients with KD, MIS-C, viral infection, or bacterial infection, while also stratifying the groups evenly with respect to hospital of origin and disease subclassification. We performed differential abundance analysis for each pairwise comparison using only the training data (Dataset S5). Next, we trained individual one-vs.-one GLMNET models with LASSO regression for each sample group (MIS-C, KD, viral infection, or bacterial infection), using the top 100 genes identified in the differential abundance analysis. These genes were selected based on gene abundance AUC and level of significance (Materials and Methods). Each of the one-vs.-one models demonstrated high performance with individual model scores classifying samples with high accuracy (test AUC: min = 0.86, max = 0.99; test accuracy: min = 0.83, max = 0.92; Fig. 4B and SI Appendix, Table S3). The union of genes used by each GLMNET model with LASSO regression generated a 109-gene panel (Dataset S6). To combine the outputs of the individual models, we trained a multiclass RF classifier using the classification scores from the one-vs.-one models. Using this framework, our multiclass machine learning model achieved high accuracy in both the training and test sets (accuracy, train = 100% and test = 80%). The bacterial infection group had the highest rate of misclassified samples, likely due to the smaller sample size and heterogeneity of both pathogen and site of infection (SI Appendix, Table S3 and Dataset S1). The high performance of the multiclass machine learning model on a relatively small number of genes points to the potential utility of cfRNA in differentiating complex inflammatory conditions in a clinical setting.

We further investigated whether cfRNA could be employed not only for classification but also for disease characterization. Since inflammation and/or infection can impact multiple organ systems, understanding organ-specific damage and function is crucial for guiding clinical management. Previous work from our group and others has shown that quantifying the cfRNA cell-type-of-origin (CTO) is a viable noninvasive method to assess cell, tissue, and organ-specific injury (17, 30). Here, we employed deconvolution of the cell types of origin of cfRNA using BayesPrism and the Tabula Sapiens human cell atlas as a reference. The CTO estimates were compared with known biomarkers and other indicators of specific organ injury or dysfunction (Materials and Methods and Fig. 5A). We first compared hepatocyte and intrahepatic cholangiocyte contributions to the cfRNA in plasma to alanine transaminase measurements (ALT, n = 141). We observed significantly elevated levels of cfRNA from hepatocytes and intrahepatic cholangiocytes in patients with high ALT (ALT > 100 IU/L) compared to patients with normal ALT (ALT < 40 IU/L, Fig. 5A). To explore cardiac function and damage, samples were categorized as having either normal or abnormal cardiac function (Materials and Methods). In the abnormal group, we observed increased levels of cfRNA from cardiac muscle cells and pericytes, likely indicative of increased cardiac cell injury or death. Interestingly, this group also exhibited elevated cfRNA levels from kidney epithelial cells, intrahepatic cholangiocytes, club cells and type I pneumocytes, and bronchial epithelium, suggesting other end organ damage associated with impaired cardiac function (Fig. 5A). Of note, patients with abnormal cardiac function did not have elevated levels of smooth or skeletal muscle cell derived cfRNA (SI Appendix, Fig. S3A). Last, to evaluate lung function and damage we focused on samples from healthy individuals and viral infection patients with COVID-19, further stratifying the COVID-19 samples by disease severity (Materials and Methods). We observed elevated levels of cfRNA from club cells and type I pneumocytes in moderate COVID-19 cases compared to healthy individuals, with a greater increase for severe cases (Fig. 5A). We observed similar trends in the levels of bronchial epithelium derived cfRNA; however, the difference was not statistically significant.

We next assessed whether cfRNA could be used to further stratify sample groups by distinguishing between COVID-19 infection from Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and infections from other viral pathogens [influenza, respiratory syncytial virus (RSV), Epstein-Barr virus (EBV), adenovirus, etc.]. We separated viral infection samples into COVID-19 and non-COVID-19 viral infections and randomly split the data into training (70%) and testing (30%) datasets. Next, we trained a GLMNET with LASSO regression model to differentiate between these groups (Materials and Methods). The trained model had high performance on the training and test sets despite including only seven genes (train AUC = 0.99 and test AUC = 0.93), demonstrating the potential of cfRNA for respiratory viral species differentiation (SI Appendix, Fig. S3 B and C).

We next simulated a cfRNA “report” for each patient in the multiclass test set (Dataset S7) and discuss four patient cases in detail to illustrate how cfRNA can be integrated to support clinical decision-making. Each report includes diagnostic predictions from the multiclass classification algorithm and predicted organ involvement levels as z-scores of CTO fractions in samples from patients with disease compared to healthy patients.

The first case is a 13-y-old male with very-high-risk B cell acute lymphoblastic leukemia in delayed intensification chemotherapy who presented with a 1-d history of fever (Fig. 5B). His laboratory assessments were notable for leukopenia, anemia, thrombocytopenia, elevated C-reactive protein (CRP), mildly elevated ALT, normal kidney function, and developed hypotension following admission. He had known household contacts with COVID-19, and his nasopharyngeal PCR was positive for SARS-CoV-2, but he did not develop respiratory symptoms. The final diagnosis was COVID-19 in an immunocompromised host. Despite the immunocompromised state of the patient and lack of respiratory symptoms common in COVID-19, the cfRNA algorithm correctly identified the patient as having a viral infection. Furthermore, the elevated levels of liver and heart derived cfRNA are interesting given the clinical manifestations of mildly elevated ALT and hypotension.

In the second case, the cfRNA model correctly predicted the patient as having MIS-C, but the classifier score for MIS-C and viral infection were very comparable, producing what could be considered a “borderline” result (Fig. 5C). The patient is a 4-y-old male who presented with a 1-wk history of abdominal pain, vomiting, and watery diarrhea, with subsequent development of fever, conjunctival injection, and swelling of the hands and feet. The patient had normal white blood cell counts, with lymphopenia, thrombocytopenia, hyponatremia, acidosis, and elevated brain natriuretic peptide (BNP), CRP, and erythrocyte sedimentation rate, with normal liver and kidney function. He was SARS-CoV-2 PCR negative but SARS-CoV-2 immunoglobulin G positive. No additional viral testing was performed. His course was complicated by hypotension, but an echocardiogram was grossly normal. The final diagnosis was MIS-C, which was the top prediction from the cfRNA model. Interestingly, viral infection was the second most probable diagnosis from the cfRNA model. MIS-C is derived from a viral infection (SARS-CoV-2), and MIS-C and viral infection can be difficult to differentiate. There are at least three possible interpretations for these results: 1) the model is correctly identifying MIS-C along with a viral signature from SARS-CoV-2, 2) the patient is misdiagnosed with MIS-C and has only a viral infection, or 3) the detection of a viral signature is due to errors in the model. We also observed elevated liver cfRNA despite the patient having normal liver enzymes along with elevated heart tissue derived cfRNA, which is consistent with the elevated BNP.

In the third case, the model incorrectly predicted the diagnosis in an immunocompromised patient, but the correct diagnosis of bacterial infection was the second most likely prediction (Fig. 5D). The patient is a 16-mo-old male with a history of Wiskott–Aldrich syndrome requiring a bone marrow transplant. He was admitted with fever, fatigue, pancytopenia with severe thrombocytopenia, elevated CRP, and blood cultures positive for Klebsiella pneumoniae. The source of the infection is unknown but was believed to be secondary to bacterial gut translocation or central line infection. His pancytopenia and thrombocytopenia were attributed to his history of Wiskott–Aldrich syndrome. He was treated with cefepime and responded clinically to a full course of antibiotics. The final diagnosis was bacterial infection, which was the second ranked prediction by the cfRNA model. Interestingly, the patient presented with some symptoms characteristic of MIS-C, the apparently erroneous top prediction made by the cfRNA model; specifically with severe thrombocytopenia and elevated CRP, which are part of the MIS-C case definition (31).

In the last case, the model predicted an incorrect diagnosis, but the results of the multiclass model would have been informative in providing the correct diagnosis (Fig. 5E). The patient presented with fever during the COVID-19 pandemic and was initially diagnosed with MIS-C. However, a family member of the patient had recently been diagnosed with KD. Given the known genetic predisposition associated with KD, a reevaluation of the original case was conducted. The patient was subsequently found to meet the American Heart Association’s case definition for KD and had initial laboratory values and a subsequent clinical course that were consistent with KD (including periungual desquamation). The patient was reclassified accordingly. Although the multiclass model did not predict the proper diagnosis of KD, the classifier scores for KD, bacterial infection, and viral infection were very similar, while the classifier score for MIS-C was much lower. Furthermore, the MIS-C vs. KD pairwise classifier accurately predicted KD. While the final diagnosis of the multiclass model was incorrect, the multiclass model and MIS-C vs. KD classifier results could have aided the clinicians in determining that the patient did not have MIS-C upon initial admission.

The UCSF Institutional Review Board (IRB) (#21-33403), San Francisco, CA; Emory University IRB (STUDY00000723), Atlanta, GA; Children’s National Medical Center IRB (Pro00010632), Washington, DC; and Cornell University IRB for Human Participants (2012010003), New York, NY, each approved the protocols for this study. All samples and patient information were deidentified for analysis and shared with collaborating institutions. At Emory University, the IRB-approved protocol was a prospective enrollment study under which parents provided consent and children assent as appropriate for age. At Children’s National Medical Center and UCSF, the IRB protocols were “no subject contact” sample biobanking protocols under which consent was not obtained and data were extracted from medical charts. At UCSD, the IRB reviewed and approved collection and sharing of samples and data (IRB #140220). Signed consent and assent were provided by the parent(s) and pediatric patient, respectively.

Samples were acquired from UCSF as previously described (17). Briefly, hospitalized pediatric patients were identified as having COVID-19 by testing positive with SARS-CoV-2 real-time PCR (RT-PCR). Residual whole blood samples were collected in EDTA lavender top tubes and diluted 1:1: in DNA/RNA shield (Zymo Research). The remaining blood was centrifuged at 2,500 rpm for 15 min and the available plasma was retained. All samples were stored at −80 °C freezer until used. Samples were acquired from Emory and Children’s Healthcare of Atlanta as previously described (17). Briefly, pediatric patients were classified as having COVID-19 by SARS-CoV-2 RT-PCR and as having MIS-C if they met the CDC case definition. Controls were healthy outpatients with no known history of COVID-19 who volunteered for specimen collection. Whole blood was collected in EDTA lavender top tubes and aliquoted for plasma extraction via centrifugation at 2,500 rpm for 15 min. Samples were stored at −80 °C and shipped on dry ice to either UCSF or Cornell for analysis. Samples were acquired from UCSD prior to any treatment in all subjects in EDTA lavender top tubes and centrifuged at 2,000 g for 10 min. Plasma was collected and stored at −80 °C. Samples were acquired from CNH as previously described (17). Briefly, pediatric patients were classified as having MIS-C if they met the CDC case definition. Whole blood samples were collected and centrifuged at 1,300× G for 5 min at room temperature. Plasma was aliquoted into a cryovial and frozen at −80 °C.

Patients were stratified as previously described (17). For the purposes of this study, patients were classified as having MIS-C by multidisciplinary teams that adjudicated whether a patient met the CDC case definition of MIS-C. COVID-19 was defined as any patient with PCR-confirmed SARS-CoV-2 infection within the preceding 14 d who did not also meet the MIS-C case definition. KD patients at UCSD met the AHA definition for complete or incomplete KD. Viral and bacterial infection patients enrolled at UCSD were adjudicated and a final diagnosis assigned by the research team (one ED clinician and one pediatric infectious disease expert) 2 to 3 mo after the acute illness to allow time for serologies, recurrence, and clinical recovery to be assessed. Patients with a self-limited illness that resolved without treatment and for whom viral studies were either negative or not done were classified as having a “viral syndrome.” Clinical data were abstracted from the medical record and submitted into a REDCap databases housed at UCSF or UCSD.

Samples were processed as described previously (17). Briefly, samples were received on dry ice, RNA was extracted, and libraries were prepared and sequenced on a NextSeq or NovaSeq Illumina sequencer. Reads were trimmed to 61 bp, and sequencing data were processed using a custom bioinformatics pipeline which included quality filtering and trimming, alignment to the human GRCh38 reference genome, and counting of gene features.

Using the sequencing data, quality control was performed by analyzing DNA contamination, rRNA contamination, total counts, and RNA degradation. DNA contamination was estimated by calculating the ratio of reads mapping to introns and exons. rRNA contamination was measured using SAMtools (v1.14). Total counts were calculated using featureCounts³⁰ (v2.0.0). Degradation was estimated by calculating the 5′-3′ bias using Qualimap³¹ (v2.2.1). Samples were removed from analysis if either the intron to exon ratio was greater than 3, if a sample had less than 75,000 total feature counts, or if the 5′-3′ read alignment ratio bias was greater than 2.

Gene transcript abundances were compared using a negative binomial model implemented using the DESeq2 R package (45). Gene transcript base mean abundance, adjusted P-value, and log2 fold change were taken from the DESeq2 Results output. Gene transcript AUC was calculated using VST transformed counts and the pROC R package (46).

Samples were partitioned for machine learning applications taking into consideration diagnosis, hospital of origin, and disease subclassification. MIS-C and COVID-19 samples were subclassified by severity, as previously defined (17). KD samples were subclassified by phenotypic subclusters, as previously defined (47). Non-COVID-19 Viral and bacterial infection samples were evenly partitioned by diagnosis and hospital of origin only.

Samples were partitioned into training, validation, and test sets at a ratio of 60:20:20, partitioning evenly based on factors such as diagnosis, severity/subgroup, and hospital of origin. Subsequently, differential abundance analysis was conducted on the training data. Genes were filtered based on specific criteria (adjusted P-value < 0.01, base mean abundance > 100, gene transcript AUC > 0.65, and absolute log2 fold change > 0.25). Filtering criteria were chosen with the intention of selecting <150 genes. We tested multiple thresholds using our training and validation sample sets and found that the model results were not sensitive to the thresholds chosen. The test set was not used for determining thresholds for gene filtering.

Raw counts for the training, validation, and test sets were individually normalized using variance stabilizing transformation, and subsets were created based on the chosen transcript features. Variance stabilization transformation was performed using the DESeq2 package, and the dispersion function from the training set was used to transform the test and validation sets. Fourteen machine learning classification algorithms were employed using the R package Caret (https://doi.org/10.18637/jss.v028.i05), including generalized linear models (GLM) with Ridge and LASSO feature selection (GLMNETRIDGE and GLMNETLASSO), support vector machines with linear and radial basis function kernel (SVMLin and SVMRAD), RF, RF ExtraTrees (EXTRATREES), neural networks, linear discriminant analysis, nearest shrunken centroids (PAM), C5.0 (C5), k-nearest neighbors, naive bayes (NB), CART (RPART), and GLM. Training was performed using fivefold cross-validation and grid search hyperparameter tuning.

For each model, classification score thresholds were determined using Youden’s index on the training data. The trained models were then employed to predict labels for the validation set, and performance was assessed using the area under the ROC-AUC. The model achieving the highest AUC on the validation set was selected as the final model and subsequently applied to the test set, which had not been utilized in any phase of model training or selection. Prediction on both the validation and test sets utilized the Youden’s index threshold derived from the training set.

Samples were partitioned into train and test sets at a ratio of 70:30, considering factors such as diagnosis, severity/subgroup, and hospital of origin. Raw counts for the training and test sets were individually normalized using variance stabilizing transformation, as implemented using the DESeq2 package, and the dispersion function from the training set was used to transform the test set. One-vs.-one GLMNET LASSO models were trained for each pairwise comparison of samples groups (e.g., KD vs. Viral Infection) using the top 150 significant features (adjusted P-value <0.05, base mean abundance > 50, and absolute log2 fold change > 1) ordered by gene transcript AUC calculated using the training data. Trained models were used to calculate classification scores for all samples in the train and test dataset, resulting in six classification scores for each sample. Classification scores were used to train a multiclass RF algorithm which assigned probability scores for each condition. The final condition with the highest probability score from the RF was assigned as the prediction for each sample.

Cell type deconvolution was performed using BayesPrism (v1.1) (48) with the Tabula Sapiens single-cell RNA-seq atlas (Release 1) (33) as a reference. Cells from the Tabula Sapiens atlas were grouped as previously described by Vorperian et al. (30). Cell types with more than 100,000 unique molecular identifiers (UMIs) were included in the reference and subsampled to 300 cells using ScanPy (v1.8.1) (49). Deconvolution values were scaled from 0 to 1 for each cell type, and medians were calculated for plotting.

For liver damage assessment, samples were separated as having normal or high ALT measurements (normal: ALT < 40 IU/L, high: ALT > 100 IU/L). The ALT measurements were taken from the same blood draw as the plasma for cfRNA processing. For cardiac function, patients were categorized as abnormal if they had abnormal EKG/ECG and/or echocardiogram results. EKG/ECG and echocardiogram results were categorized as abnormal in the context of the patient narrative and final interpretation of the studies. For endothelial damage, samples were separated as either having KD/MIS-C or bacterial/viral infection. COVID-19 patients were determined to have moderate or severe disease using the following criteria:

Moderate: The patient must have been hospitalized due to COVID-19 respiratory disease and/or any systemic/nonrespiratory symptoms attributed to COVID-19 (e.g., neonatal fever, dehydration, new diagnosis diabetes, acute appendicitis, necrosis of extremities, diarrhea, encephalopathy, renal insufficiency, mild coagulation abnormalities, etc.).

Severe: The patient must have been hospitalized for COVID-19 with either high-flow oxygen requirement [high-flow nasal cannula, continuous positive airway pressure, bilevel positive airway pressure, intubation with mechanical ventilation, or extracorporeal membranous oxygenation] and/or evidence of end-organ failure [acute renal failure requiring dialysis, coagulation abnormalities resulting in bleeding or stroke, diabetic ketoacidosis, hemodynamic instability requiring vasopressors], and/or dying from COVID-19. These patients were almost always admitted to the intensive care unit (ICU).

Training was done using the same method as the KD vs. MIS-C model. Briefly, viral infection samples were partitioned into train and test sets at a ratio of 70:30, considering factors such as COVID-19 status, severity/subgroup, and hospital of origin. Subsequently, differential abundance analysis was conducted on the training data. Genes were filtered based on specific criteria (adjusted P-value < 0.01, base mean abundance > 100, gene transcript AUC > 0.65, and absolute log2 fold change > 0.25), and the top 150 genes, as ordered by gene transcript AUC, were selected as inputs for machine learning analysis.

Raw counts for the train and test sets were individually normalized using variance stabilizing transformation, and subsets were created based on the chosen transcript features. A GLMNET with LASSO regression was trained on the training set using fivefold cross-validation and grid search hyperparameter tuning. Classification score thresholds were determined using Youden’s index on the training data. The trained models were then employed to predict labels for the test set using the Youden’s index threshold derived from the training set.

The programming language R (v4.1.0) was utilized for all statistical analyses. Statistical significance was assessed through two-sided Wilcoxon signed-rank tests and Mann–Whitney U tests, unless specified otherwise. Machine learning algorithms were trained using the Caret R package, and pipelines were run using the Snakemake workflow management system (50, 51). In boxplots, boxes denote the 25th and 75th percentiles, the band within the box signifies the median, and whiskers extend to 1.5 times the interquartile range of the hinge. The alignment of all sequencing data was performed against the GRCh38 Gencode v38 Primary Assembly, with feature counting conducted using the GRCh38 Gencode v38 Primary Assembly Annotation (52).

C.J.L., M.D., R.L.D., C.A.R., J.C.B., C.Y.C., and I.D.V. designed research; C.J.L., A.B., J.L., E.B. processed samples; C.J.L., V.S., A.S.-G., A.B., W.S., J.N., M.E.W., M.O., M.A.G., J.-H.C., H.-M.H., H.W., J.K., C.S., A.H.T., J.C.B., C.Y.C., and I.D.V. analyzed data; W.S., J.N., M.E.W., M.O., M.A.G., P.E.M.K.D.R.G., C.H.A.R.M.S., J.-H.C., H.-M.H., C.S., A.H.T., M.D., R.L.D., C.A.R., J.C.B., and C.Y.C. identified and collected patient samples and clinical metadata; and C.J.L., J.C.B., C.Y.C., and I.D.V. wrote the paper.

C.J.L and I.D.V are inventors on submitted patents pertaining to cell-free nucleic acids (US patent applications 63/237,367 and 63/429,733). I.D.V. is listed as an inventor on submitted patents pertaining to cell-free nucleic acids (US patent applications 63/237,367, 63/056,249, 63/015,095, 16/500,929, 41614P-10551-01-US) and receives consulting fees from Eurofins Viracor. C.A.R. has received institutional support from ModernaTX, Inc., Pfizer Inc., BioFire Inc., GSK plc, MedImmune, Micron Technology Inc., Janssen Pharmaceuticals, Merck & Co., Inc., Novavax, PaxVax, Regeneron, and Sanofi Pasteur. She is co-inventor of patented RSV vaccine technology which has been licensed to Meissa Vaccines, Inc. C.Y.C. receives grant funding for research unrelated to this work from the Bay Area Lyme Disease Foundation, the Chan-Zuckerberg Biohub, and Abbott Laboratories, Inc.