Virus-induced genetics revealed by multidimensional precision medicine transcriptional workflow applicable to COVID-19.

Precision medicine requires the translation of basic biological understanding to medical insights, mainly applied to characterization of each unique patient. In many clinical settings, this requires tools that can be broadly used to identify pathology and risks. Patients often present to the intensive care unit with broad phenotypes, including multiple organ dysfunction syndrome (MODS) resulting from infection, trauma, or other disease processes. Etiology and outcomes are unique to individuals, making it difficult to cohort patients with MODS, but presenting a prime target for testing/developing tools for precision medicine. Using multitime point whole blood (cellular/acellular) total transcriptomics in 27 patients, we highlight the promise of simultaneously mapping viral/bacterial load, cell composition, tissue damage biomarkers, balance between syndromic biology versus environmental response, and unique biological insights in each patient using a single platform measurement. Integration of a transcriptome workflow yielded unexpected insights into the complex interplay between host genetics and viral/bacterial specific mechanisms, highlighted by a unique case of virally induced genetics (VIG) within one of these 27 patients. The power of RNA-Seq to study unique patient biology while investigating environmental contributions can be a critical tool moving forward for translational sciences applied to precision medicine.

INTRODUCTION

Precision/personalized medicine represents the concept of rapidly identifying altered biology within an individual and using the findings to guide therapy. Much of the promise of this new strategy is obscured by the inadequacy of available technology to discern altered biological pathways in a noninvasive fashion and concomitantly view multiple different facets of altered biology. Whenever a disease process is multisystemic such as infection, neoplasm, or toxin-mediated, blood has potential to provide early detection to characterize injury. RNA-Seq has been implemented in many trait mapping projects, using large cohorts of specific phenotypes/trait groups relative to sedation controls to identify biomarkers. It remains to be shown how the millions of dollars of integrated RNA-Seq work can be translated into individual patients of broad phenotypes to guide precision medicine. The power of RNA-generated data likely is in its potential to effect therapeutic changes in individual patients with diverse clinical presentations, such as multiple organ dysfunction syndrome (MODS). Coronavirus disease 2019 (COVID-19) infections highlight this point well, where not every patient with COVID-19 needs to be treated for the MODS that results in lethality, just those that advance into sepsis, where biomarkers of these pathways are critical for earlier treatment of patients.

The pediatric intensive care unit (PICU) sees diverse phenotypes involving environmental induction (virus, bacteria, trauma), resulting in individual outcomes that could benefit from an expeditious diagnostic and therapeutic timeline. This is ideal for the development and implementation of rapid precision/personalized medicine. Multitime point analyses that can dissect disease by using omic technologies in n = 1 biology is currently lacking in the PICU. Patients with MODS have high rates of morbidity and mortality and represent some of the most complex cases within the PICU. Using a combination of PAXgene tubes with ribosomal reduced RNA-Seq, we present here a dissection of both cellular and acellular RNA signatures of PICU MODS patients, suggesting utility for transcriptomics in precision medicine.

Over the last decade, our ability to extend life in patients with MODS through interventions like dialysis and extracorporeal membrane oxygenation (ECMO) has improved, but there remains an inability to serially track the dyshomeostasis induced by disease from trigger to organ injury, an end point for a variety of infectious and noninfectious agents. To halt injury to organs and promote recovery, the clinical team needs rapid, actionable clinical information within the short treatment window of any ICU. That information needs to include simultaneously the etiology, i.e., was it infectious or noninfectious, and the extent of homeostatic disruption. Bacterial cultures take at least 24 h to grow and 72 h to speciate. Viral panels check only for common viruses, requiring insights of what to test. Both viral and bacterial analyses are often ordered only when the phenotype is suspected of matching, meaning unique responses are often missed for proper diagnosis. At the same time, the treating physicians are often unsure of the patient’s clinical state, such as whether the patient is immunosuppressed or experiencing an overwhelming immune and inflammatory response (12, 30, 31). We are seeing that even in epidemics from an infectious agent like COVID-19, there are differing clinical phenotypes that lead to MODS, requiring changes in therapeutic course (14, 15). In any critically ill patient with MODS, the first 48 h are spent by the clinical team supporting organ function with little ability to make effective changes to the clinical course. At present, our perception of organ damage is limited to a few organs that have validated, but imperfect, means of measurement such as the heart, lung, kidney, liver, and brain, missing markers for many other tissues. The transcriptome is an attractive tool, offering a snapshot of the interaction between affecting pathology and patients’ biology. Thus, transcriptomics is a practical, real-time approach to diagnosis and management, which, with decreasing time of processing and improved computational capacity to provide insight, could be a powerful tool to predict disease etiology and clinical course. Doing so reduces the diagnostic odyssey in a population that can ill afford delay, providing an agnostic guide to patient-specific disease. The individual responses seen in each MODS patient make this group a prime cohort for personalized/precision medicine strategy development (5), as we lay forth below.

METHODS

Sample Collection

Institutional Review Board approval [2016-062-SH/Helen DeVos Children’s Hospital (HDVCH)] was obtained for a precision medicine initiative to characterize patients admitted to our PICU with multiorgan dysfunction. Blood samples were collected at up to three independent time points, at baseline (day 0), greater than 48 h to less than 72 h (day 3), and greater than 7 days (day 8). These time points were collected to reflect the trajectory of critical illness with measurements in the acute, stabilization, and recovery phase of illness. The clinical team extracted from the clinical notes, diagnosis, transfer status from other hospital, genomic diagnosis/comorbidity for patient (if available), 2 yr history within the hospital system, pediatric logistic organ dysfunction (PELOD) score, race, age, sex, clinically identified infections, organs experiencing failure, and several clinical measures when available.

Patients

The sedation patients were less than 18 yr of age presenting for routine sedation in an outpatient clinic without an inflammatory disease. The definition for MODS was taken from the original definition (47) in which, among other criteria, MODS patients were clinically identified with two or more organs in failure. ECMO is defined as MODS patients that required ECMO therapy.

PAXgene RNA

PaxGene RNA tubes were stored at room temperature for at least 2 h following the addition of blood, overnight at −20°C, and at −80°C for long-term storage. PAXgene tubes were thawed at room temperature for 2 h and processed using Qiagen’s QIAsymphony PAXgene Blood RNA Kit.

RNA-Seq

RNA-Seq libraries were generated using RiboZero-Globin rRNA kit (Illumina), KAPA RNA HyperPrep Kit (v1.16) (Roche Life Sciences), and Bio Scientific NEXTflex Adapters (Bioo Scientific). Libraries were pooled and run over seven flowcells using 150 bp sequencing kit (v2) in paired-end on an Illumina NextSeq 500.

Human Bioinformatics

Paired-end reads (fastq) were aligned in salmon_0.14.1 (45) using the Homo sapiens Gencode 31 for transcripts or our custom single gene to sequence transcript library for genes. Mapped transcript per million (TPM) for all samples were processed through NetworkAnalyst3.0 (66) using DESeq2 (40) to identify gene-level annotations that differ between sedation control, MODS, and ECMO groups. Pathways identified in Network Analyst to have variable gene expressions with an adjusted p value (< 0.05) were extracted for Gene Ontology- and KEGG-associated genes. Pooling all day 0 TPM values for transcripts, we calculated z-scores to identify day 0 samples with >2 standard deviations above the mean of all and a value >1 TPM. The biotype data for the transcripts in each of these samples >2 standard deviations were compared with the percentage of biotype annotations in the entire transcriptome to calculate transcript enrichment. The genes corresponding to the transcripts >2 standard deviations were assessed with STRING (51) for pathway enrichment. For patients with transcriptomes at day 0 and day 8, Log2 fold change was calculated for day 0 relative to day 8 on that patient and compared with day 0 of the patient relative to day 0 of all other patients. Neutrophil-specific genes were taken from the LM22 matrix of CIBERSORTx (43, 44), and cytokines identified from protein annotations on UniProt for cytokines. Heatmaps were created with the Broads Morpheus tool with each gene normalized over all samples, and clustering performed by one minus Pearson correlation. CIBERSORTx digital cytometry cell fraction imputation (43, 44) was performed using the gene mapping on the LM22 signature matrix without batch corrections and 500 permutations.

Infection Bioinformatics

Using only the first of the paired-end reads, we performed bacterial and viral read mapping with kraken2-2.0.7-beta (61) used against the full transcriptome database of the tool downloaded on June 24, 2019. Samples were normalized to every million mapped reads to human and clustered based on patient ID. Z-scores were calculated for each viral or bacterial species annotations across all time points for all patients to determine outliers. When needing a secondary read alignment, we obtained Burrows-Wheeler aligner (BWA) (36) alignment to individually downloaded transcriptomes.

Organ Biomarker Database

Using GTEx (28) tissue annotated expression data sets, we filtered for genes that are unique to each tissue/organ. All sex-specific organs were removed from the database to begin. Genes had to have an expression value >5 and an enrichment of 10-fold in one tissue over the average of all. Each annotated tissue needed to have at least five genes identified to be included in the search. From the patient samples included genes for annotation had to have a 20-fold elevation in one of the samples relative to the average of all samples. The TPM values for each tissues were added together for sample and the percent of genes with a value above 0 determined.

RNASEH2B Patient Assessment

A detailed analysis of patient 24 was performed. A 2 yr follow-up PAXgene tube sample was isolated from the patient, and RNA extracted as before. Illumina RNA-Seq was performed on the isolated RNA as before. Leftover RNA was run on Nanopore RNA-Seq using PCR-cDNA Sequencing kit and MinION flow cell, to compare values with RNA obtained from a control sample. We compared Illumina and Nanopore differentially expressed transcripts with isolate overlapping data, running the gene list through NetworkAnalyst3.0 for enriched pathways. The last bit of isolated RNA was sent to iRep for Immune Repertoire analysis of seven chains of CDR3 of the immune system. Reads from day 0 were reassessed for viral mapping using Taxonomer (https://www.taxonomer.com) and the Epstein-Barr virus (EBV) transcriptome aligned with BWA (36). Rapid whole-exome sequencing on the patient and both parents was performed clinically through GeneDx. Assessment of variant outcomes in RNASEH2B was done according to our previously published workflow (46). Reads from all of the patient’s RNA-Seq were aligned onto the RNASEH2B transcript using BWA (36) followed by calculation of allelic bias using the variant inherited from the father. Reads were then manually viewed in Sequencer software to interpret the splice variant site for protein prediction.

RESULTS

Patient Cohort Demographics

The cohort consists of four sedation controls (Sedation), 17 MODS, and 6 MODS patients that needed extracorporeal membrane oxygenation (defined as ECMO throughout). Diverse ethnicities are represented in the data: Caucasian, 16; Hispanic/Caucasian, 6; Black, 4; East Asian, 1; and Hispanic/Black, 1. While there were mostly even numbers of male/female in control and MODS, ECMO patients tend to be male. Renal failure (89% MODS, 100% ECMO) and liver failure (30% MODS, 50% ECMO) are higher in ECMO than MODS group. Eight of the patients were transferred from local hospitals to HDVCH to receive more advanced PICU support at a quaternary children’s hospital. A total of 12 patients had been diagnosed with additional comorbidities, including multiple rare diseases diagnosed via genomics. Blood draws were done on each patient followed by RNA-Seq, with six patients having only a day 0 transcriptome, 8 with day 0/3, and 13 with day 0/3/8 transcriptomes (Fig. 1) based on length of stay in the PICU. Throughout the paper and figures each patient is referred to by deidentified number, with details of grouping added (sample:group:day:sex).

Fig. 1.

Transcriptomics workflow used to study pediatric intensive care unit (PICU) patient precision medicine. Shown in red are returnable data sets from the RNA-Seq. ECMO, extracorporeal membrane oxygenation; MODS, multiple organ dysfunction syndrome.

Transcriptome Comparisons with Known Markers

Going from large cohort studies to individual patients, RNA-Seq analysis has been very challenging to date. While large sample sizes are used to derive significance and powerful marker genes, it impedes understanding of RNA-Seq signatures in each patient, mainly when diverse phenotypes are present. The collection of 61 samples from three time points of 27 different PICU patients highlights this. Unsupervised clustering of sample gene-level annotations (27,854) over multiple days shows minimal bulk overlap of the samples (Fig. 2). Using the day 0 RNA-Seq of the 27 patients, we find some clustering of the sedation controls relative to the MODS/ECMO patients (Fig. 2).

Fig. 2.

61 transcriptomes analyzed from 27 patients. Principal component analysis of transcriptomes annotated at the gene level for all 3 days (day 0 green, day 3 blue, day 8 magenta) (A) or on day 0 for the various groups (Sedation control blue, MODS green, MODS with ECMO red) (B). C: using the five PERSERVERE genes, we quantify the total transcripts per million (TPM, x-axis) and the percent of the five genes expressed (y-axis). Two outliers are labeled (sample:group:day:sex). Expression of various neutrophil genes (D) or cytokines (E) shown on the left as the sum of their TPM (x-axis) and the % of genes in the group expressed (y-axis). On the right is a heat map of each sample and the various genes. Clustering is based on one minus Pearson’s correlation. ECMO, extracorporeal membrane oxygenation; F, female; M, male; MODS, multiple organ dysfunction syndrome.

Validated biomarkers for pediatric sepsis severity (60) were examined and found in two samples to have significant elevation (Fig. 2), suggesting that sepsis severity is not always the driver of MODS. The two patients with elevated sepsis severity markers had severe outcomes, as seen in patient case descriptions. Another set of genes thought to contribute to MODS are those associated with neutrophil biology (9), where we show the elevation of several neutrophil-specific genes (43) within our cohort (Fig. 2). Patient 11 day 0 has the highest neutrophil expression coming from elevation of S100A12, FPR1, VNN2, LST1, AQP9, FPR2, TNFRSF1, MNDA, FCGR3B, which are connected particularly with neutrophil motility (38). Patient 8 day 8 (passed away from unidentified infection) has the second-highest neutrophil profile but utilizing a completely different subset of genes than patient 11, including DPEP2, APOBEC3A, MEFV, FFAR2, MAK, P2RY13, which can be found linked to sepsis activation and tissue recruitment (20, 27). Cytokines are known to contribute to MODS (58), where 19 blood cell transcribed cytokines show elevation in patients (Fig. 2). Patient 8 day 8, in addition to the unique neutrophil activation, has the highest cytokine activation, including IFNK, LTB, CD70, IL18, CCL5, BMP7, GPI, GDF5, IL10, CCL22, CXCL16, TNFSF14, IL11. Patient 20 day 0 has the second-highest cytokine activation, including CCL5, BMP7, GPI, GDF5, IL10, CCL2, CXCL, TNFS, IL11.

Cell Composition

From the gene level annotation, we assessed cell composition using CIBERSORTx digital cytometry (43, 44). Of the 22 cell types annotated, there is a broad spectrum of composition predictions (Fig. 3). Monocytes and neutrophil genes make up the most substantial portion in most samples, with other cell types having outliers from the cohort averages. Two patients on day 3 are outliers for monocyte gene levels (patients 26 and 27). Two patients on day 8 have continued elevated levels of neutrophil genes (patients 21 and 20). Memory B cell genes have two samples as outliers on day 0 (patients 21 and 5). Several samples show elevated CD8 T cell genes, including patient 24 day 0. One patient on day 0 shows elevation of dendritic cell genes (patient 20) and one on day 8 for Tregs (patient 4). Eosinophil-elevated genes are seen in patient 25 on day 0 and 3, with additional elevation of patient 6 on day 0 and patient 2 on day 3. Four of the patients show elevation of day 0 CD4 memory T cell genes, including patient 24 and 19 and two of the sedation controls (patients 17, 16). Two samples have elevation of mast cell genes, day 3 patient 24, and day 0 patient 1. Day 0 samples show an elevation of M0 macrophage genes with highest levels seen in patients 4, 1, and 7 with elevation in patient 22 on day 3 and patients 5 and 15 on day 8.

Fig. 3.

Blood cell composition of each sample. CIBERSORTx plots of each sample deconvoluted cell identifies shown as box and whisker plots for 16 of the cell types for day 0 (red), day 3 (orange), and day 8 (gray). Outliers are labeled (sample:group:day:sex). ECMO, extracorporeal membrane oxygenation; F, female; M, male; MODS, multiple organ dysfunction syndrome.

Bacterial/Viral Mapping

Many of the patients were known to have clinical infections, and the blood draws for RNA-Seq were performed in PAXgene tubes known to support isolation of bacterial and viral RNA. Therefore, we assessed RNA reads associated with bacteria and viruses using Kraken2, with 20,266 annotations. Mapped reads were normalized to human reads with values around 1,000–2,000 reads for bacteria and viruses (Fig. 4). Patients 5 and 11 had the highest bacterial read mapping.

Fig. 4.

Bacterial/viral mapping from RNA-Seq. A: box and whisker plot for domain annotations normalized to every million reads of human for day 0 (red), day 3 (orange), and day 8 (gray). Box and whisker plots for the highest mapped bacteria (B) and samples with a significant outlier (C) reads mapped relative to the number of human mapped reads. Outliers are labeled (sample:group:day:sex). ECMO, extracorporeal membrane oxygenation; F, female; M, male; MODS, multiple organ dysfunction syndrome.

For bacteria (Fig. 4), this workflow identified a Staphylococcus aureus infection clinically confirmed that showed up at day 8 in patient 7 and day 0 in patients 2 and 12, which are some of the highest values mapped for any species. Several species show up in all samples, including Mycobacterium kansasii, Staphylococcus simulans, Rathayibacter rathayi, Stenotrophomonas maltophilia, Bacillus cereus, Clostridium, Escherichia coli, Klebsiella pneumoniae, Streptomyces rubrolavendulae, and Cupriavidus taiwanensis. Additional insights from the top mapped bacteria include normal flora elevation of Polynucleobacter necessarius and Bordetella parapertussis in patient 24 suggested to have issues in antigen processing and presentation (case study presented below for RNASEH2B), multiple Streptococcus strains (including pyogenes) identified in patients 11 and 5 that were culture positive, Pandoraea faecigallinarum in patient 10, and Cryobacterium arcticum in patient 10 day 0.

Focusing on bacteria and viruses that are extreme outliers within one or more samples, we identify Human gammaherpesvirus 4 (Epstein-Barr, EBV) in patient 24 that was clinically confirmed with viral PCR (Fig. 4). In certain patients, infections that were unsuspected such as Measles morbillivirus (patient 6), EBV (patient 9), and Human betaherpesvirus 7 (patient 20 and 19) that may have been associated with disease outcomes were identified. The M. morbillivirus was the most surprising outcome, with bioinformatics confirming the etiology as that from vaccination and the clinical teams’ investigation of clinical notes revealed a measles vaccine was administered 4 days before enrollment in the study. Patient 8, a patient with developmental delay and limitation of support, had a suspected infection that clinicians were unable to identify, ultimately resulting in the patient passing away from MODS several days following the study. Using the bacterial mapping, we isolated unique signatures of Cloacibacillus in all 3 days with a further elevation at day 8, a bacteria incredibly challenging to grow and shown to associate with multiorgan pathology similar to this patient (24). While reads align to some of these transcriptomes, reads may belong to closely related species not within our transcriptome database.

Organ Damage Biomarkers

While assessing the transcripts within the RNA-Seq samples, we noticed the presence of multiple albumin (ALB) transcripts in patients 24 and 27 at day 0 (Fig. 5). The ALB gene is known to be highly specific to the liver, and its presence was a surprise. Both patients were annotated with liver failure, marked by bilirubin and AST blood elevation. Reads in the ALB gene at day 0 of patient 27 cover 89% of the gene transcript with 10 or more reads, showing these read annotations are not erroneous (Fig. 5). Interestingly, both of these patients had a day 3 spike in AST blood levels, not at day 0 when ALB is high (Fig. 5). This suggests that ALB RNA in the blood could have potential to mark liver damage. Moreover, from the Human Protein Atlas, we know that 15 genes are unique to bone marrow. These genes show that all three time points of patient 8 are elevated for blood markers of bone marrow cells, with day 8 the most elevated, and 2 days elevated in patient 12 (Fig. 5).

Fig. 5.

Organ-specific biomarker transcripts. A: mapping of each RNA-Seq to transcripts of ALB, shown as the summed TPM count (x-axis) and the percent of ALB transcripts identified (y-axis). B: BWA alignment of reads from patient 27 day 0 to ALB (NM_000477.7). C: correlation of ALB mapped reads to blood measured AST at day 0 (black) and day 3 (red). The black line connects the two points for patient 24 and red line for patient 27. D: Human Protein Atlas annotated bone marrow-specific genes assessed in each patient for the summed TPM and the % of genes identified above 0. E: GTEx tissue-specific genes annotated for 12 tissues. Shown for each tissue is the number of genes used for annotations next to the tissue label. The x-axis shows the summed transcript per million (TPM), and the y-axis shows the percent of tissue-specific genes with expression >0. Genes, listed in order of expression value, are shown in red for the top nonpatient 20 sample. Genes in blue are shown for small intestine, kidney, nerve, and stomach, where a small number of genes have very high expression. ECMO, extracorporeal membrane oxygenation; F, female; M, male; MODS, multiple organ dysfunction syndrome.

To begin understanding the patient-specific tissue responses in each of the transcriptomes, we developed a list of genes that are unique to tissues from GTEx (39) and found only a small subset of patient RNA-Seq to have detectable transcripts for these tissues (Fig. 5). Three patients are identified as extreme outliers at day 0 (patient 20, 8, and 19) in all of the tissues. At day 0, these are the same three patients that show elevation of both Shewanella baltica and Shigella boydii bacteria in Fig. 4. Top genes identified within these three patients are labeled on the figure. Tissues with additional samples identified include liver, pituitary, small intestine, kidney, nerve, and stomach. Five genes outside of ALB are seen elevated in patients 27 and 24 (APOA1, APOH, SARDH, NPC1L1, CPS1) for the liver. For small intestine, two genes (DEFA5, HRASLS2) are seen elevated at day 8 in patient 8 and day 0 in sedation control patient 17. Three patients at day 0 have elevation of the kidney tubule gene NAT8 (patient 10, 18, 24). Patient 9, who had necrotizing fasciitis, had the elevation of RNA for PGA3 that is highly associated with stomach expression.

Transcript Biotype Annotations

A strategy of mapping biological outcomes within each patient was to map transcripts that are two standard deviations above other samples annotated to 226,030 transcripts of GENCODE followed by the assessment of transcript biotypes (Fig. 6) and enriched pathways. Again patient 20 day 0 RNA-Seq is an outlier with 18,062 transcripts elevated, followed by patient 6, 8, and 9 as outliers. Annotation of the biotypes identifies several samples enriched, suggestive of altered biology, including patient 20 with retained introns, patient 7 for rRNA, and patient 12 for scaRNA, snoRNA, and mitochondrial tRNA. Utilizing the list of genes elevated for each patient with gene ontology (GO) enrichment suggests multiple overlaps with clinical findings. Enriched genes associated with biological pathways are listed in the case presentations below.

Fig. 6.

Transcript biomarker outliers. Far right panel: the number of transcripts within a transcriptome that are more than two standard deviations from the mean with outlier samples identified. The remaining panels show biotype annotations of elevated transcripts for each patient with the red bar represented the expected frequency of the biotype based on total presence in the transcriptome. ECMO, extracorporeal membrane oxygenation; F, female; M, male; MODS, multiple organ dysfunction syndrome.

Building Precision Medicine Dynamic Insights from Transcriptomics

A total of 13 of the patients had RNA-Seq for three time points, allowing for comparison of day 0 versus day 8 transcriptomes and average of all three time points relative to other patients (Fig. 7). Several of the patients, including patient 20, have an elevation or suppression of many transcripts in day 0 relative to day 8 (Fig. 7, left) and transcripts elevated in all time points of the patient relative to others (Fig. 7, right). GO enrichment again revealed a remarkable overlap to patient clinical records. Of note, patient 24 had an elevation of 67 transcripts, including multiple miRNA, snoRNA, and snRNA at day 0 that were unique to the patient, the only time this occurred, suggesting a unique biological outcome of potential syndromic biology. Integration of all data for each patient case understanding is presented below.

Fig. 7.

Patient unique transcripts for day 0–8 and patient to all genes. Left: genes up (gray) or down (red) in patients at day 0 RNA-Seq relative to day 8. Right: genes up (gray) or down (red) in all 3 days RNA-Seq of patient relative to all other patients. Shown below is the plot of both day 0/8 (x-axis) and patient to all (y-axis). Labeled for each patient are the primary diagnoses and the GO enriched terms that overlap to phenotype (red and blue). The most surprising observation was the HHV-7 detected infection in patient 20, who had sarcoma, that aligns to genes decreased at day 0 that correlate to herpes-associated sarcoma. ECMO, extracorporeal membrane oxygenation; EBV, Epstein-Barr virus; FDR, false discovery rate; HLH, hemophagocytic lymphohistiocytosis; PPI, protein-protein interaction; SRP, signal-recognition particle. On the bar graph, A, B, and C are shown below for each patient.

Patient Precision Medicine Insights

Using the data for each patient, we detail the 23 MODS/ECMO cases, integrating the transcriptome data with notes and primary diagnosis. Organs in failure are L for liver, R for renal, and C for cardiovascular with MV for mechanical ventilation.

1 (ECMO: sepsis/acute myocarditis; R, MV, C; PELOD-12).

Clinical enterovirus and Staph infection; high mapping Staphylococcus simulans and Staphylococcus aureus day 0; 1,485 day 0 elevated transcripts [GO: ribosomal protein genes (32)]; elevation of mast cells at day 0; high levels of M0 macrophages at day 0.

2 (MODS: sepsis/dilated cardiomyopathy; L, MV, C; PELOD-11).

Clinical rhino-/enterovirus; high detected Staphylococcus aureus at day 0; 383 day 0 elevated transcripts (GO: viral carcinogenesis); elevation of eosinophils and naïve CD4 T cells at day 3.

3 (MODS: sepsis; L, R, MV; PELOD-12).

Clinical ornithine transcarbamylase deficiency, newborn seizures, no clinically identified pathogens; high day 0 mapping of Pasteurella multocida and Tsukamurella paurometabola; P. multocida has been observed to cause newborn seizures and sepsis (50); T. paurometabola has been observed to cause sepsis (33); 37 elevated day 0 transcripts.

4 (MODS: sepsis; MV, C; PELOD-11).

Died; Sturge Weber syndrome; clinical Staph infection; detected S. aureus; enrichment of known chronic based infectious agent (22) Pandoraea apista and the highest level of K. pneumoniae at day 3; 601 up (GO: myeloid leukocyte activation) and 465 down (GO: mRNA metabolic process) transcripts day 0 to day 8; patient elevation of 16 down transcripts in all three days (GO: defensin/corticostatin family); day 8 spike in Tregs, which are associated with systemic inflammatory response (48); highest level of M0 macrophages at day 0.

5 (MODS: sepsis; R, MV; PELOD-22).

Clinical group A beta hemolytic Streptococcus; enrichment of Streptococcus pyogenes; 287 up (GO: neutrophil degranulation) and 1,549 down (GO: metabolic process) transcripts day 0 to 8; elevation of memory B cells at day 0; high levels of M0 macrophages at day 8; highest bacterial mapped reads.

6 (MODS: sepsis; MV; PELOD-11).

Schaaf Yang syndrome (MAGEL2 identified variants); clinical Pseudomonas; only patient with detectable M. morbillivirus reads (received vaccine before admission); detection of Pseudomonas phage PaP1 and Pseudomonas entomophila a known hemolytic soil bacterium (53); 7,466 day 0 elevated transcripts [GO: ribosomal protein genes (32)]; elevation of eosinophils at day 0.

7 (ECMO: sepsis; R, MV; PELOD-20).

Vascular sling; clinical Staph infection; day 8 highest S. aureus detection; 199 up (GO: immune system) and 443 down (GO: cellular process) transcripts at day 0 to 8; high level of M0 macrophages at day 0.

8 (MODS: sepsis; MV; PELOD-20).

Died; seizure disorder with identified variants in CRELD1 (p.Q320RfsX25, p.C192Y); clinically confirmed Metapneumo; Raoultella and Cloacibacillus enriched in all 3 days of RNA-Seq with Cloacibacillus highest at day 8 [Cloacibacillus has been associated with multiorgan pathology in multiple cases (24)]; elevation of 72 transcripts at all days relative to other patients suggestive of neutrophil degranulation (1e-12), leukocyte activation (FDR 6e-13), and linked with known neutrophil protease associated organ pathologies (3, 6, 55, 62). Second highest PERSEVERE gene elevation. Neutrophil activation of DPEP2, APOBEC3A, MEFV, FFAR2, MAK, P2RY13. Cytokine activation of IFNK, LTB, CD70, IL18, CCL5, BMP7, GPI, GDF5, IL10, CCL22, CXCL16, TNFSF14, IL11; day 0 elevation of both Shewanella baltica and Shigella boydii; spike in bone marrow-specific genes at all 3 days with the highest at day 8; elevation of tissue damage biomarkers at day 0.

9 (MODS: sepsis; MV; PELOD-12).

Necrotizing fasciitis; clinical Streptococcus pyogenes infection, 45 units/L aspartate aminotransferase, 2 g/dL albumin, indicating possible liver damage; high mapping of Clostridium perfringens on day 0 and 3, C. perfringens is associated with necrotizing fasciitis and hypoxia (10); elevation of PGA3 associated with stomach at day 0; 5605 elevated transcripts of which some indicate HIF-1 signaling pathway enrichment involved in metabolic adaptation to hypoxia; abnormally large proportion of elevated transcripts are retained introns; spike in mapping of human gammaherpesvirus 4 on day 3 and observed at day 0 at lower levels.

10 (MODS: sepsis; MV; PELOD-5).

Clinically identified coronavirus; day 0 high mapping of Pandoraea faecigallinarum known to be a chronic lung colonizer (23) and several Staphylococcus infections; high level of mapped reads to Cryobacterium arcticum at day 0; high level of NAT8 at day 0 associated with kidney tubule cells; spike in COL28A1 at day 0 associated with nerve tissue.

11 (MODS: sepsis; R, MV; PELOD-11).

Clinical Streptococcus pyogenes infection; day 0 high mapping of S. pyogenes as well as other Streptococcus infections; day 0 high mapping of Listeria monocytogenes and S. aureus; Neutrophil activation of S100A12, FPR1, VNN2, LST1, AQP9, FPR2, TNFRSF1, MNDA, FCGR3B; second-highest day 0 bacterial mapping.

12 (ECMO: sepsis; R, MV; PELOD-11).

Clinical human orthopneumovirus; day 0 high mapping of Staphylococcus aureus and Burkholderia; spike in bone marrow-specific genes at days 0 and 3; 23 decreased transcripts compared with other patients (GO: Staphylococcus aureus infection); 210 elevated day 0 transcripts of which snoRNA, scaRNA, and mt tRNA are abnormally high.

15 (ECMO: sepsis; L, R, MV; PELOD-11).

Clinical Haemophilus infection and human orthopneumovirus infection; day 0 high mapping of Haemophilus influenzae, Acinetobacter baylyi [known to be opportunistic with other pathogens (18)], and Pasteurella multocida [known infection from animal scratches (4)]; high level of M0 macrophages at day 8.

18 (MODS: sepsis; R, MV; PELOD-20).

Clinical Serratia marcescens, rhino-/enterovirus, coronavirus; all 3 day high mapping of Dichelobacter nodosus; day 0 high mapping of Paenibacillus sp. FSL R5-0912 and Lysinibacillus fusiformis, which have been observed to lead to sepsis (59); high mapping of Staphylococcus aureus; high level of NAT8 at day 0 associated with kidney tubule cells; 419 transcripts increased from day 0 to 8 (GO: immune system process).

19 (MODS: sepsis; R, MV, C; PELOD-11).

Clinical methicillin-resistant Staphylococcus aureus (MRSA), Staphylococcus aureus, rhino-/enterovirus infection; day 0 high mapping of Staphylococcus simulans and human betaherpesvirus 7; day 0 elevation of both Shewanella baltica and Shigella boydii; elevation of tissue damage biomarkers at day 0.

20 (MODS: acute kidney injury; R, MV; PELOD-11).

Clinical alveolar rhabdomyosarcoma and common variable immunodeficiency with predominant immunoregulatory T cell disorders; high mapping human betaherpesvirus 7 (HHV-7) day 0 known to infect CD4+ T cells (41); day 8 elevation of Kocuria palustris, which is known to be associated with catheter-related bacteremia (35); 18,062 day 0 elevated transcripts, a high proportion of which are retained introns suggesting cancer (65); 700 transcripts down day 0–8 (GO: Kaposi's sarcoma-associated herpesvirus infection); HHV-7 z-score day 0/3/8: 6.2/−0.2/−0.2; Cytokine activation of CCL5, BMP7, GPI, GDF5, IL10, CCL2, CXCL, TNFS, IL11. Day 8 continued elevation of neutrophils; elevation of dendritic cells at day 0; day 0 elevation of both Shewanella baltica and Shigella boydii; elevation of tissue damage biomarkers at day 0.

21 (MODS: sepsis; MV; PELOD-11).

MECP2 duplication at Xq28 region; clinical rhino-/enterovirus; 56 units/L clinical aspartate aminotransferase, 2 g/dL albumin, indicating possible liver damage; MECP2 duplication has been linked to immunodeficiencies (7); day 8 continued elevation of neutrophils; elevation of memory B cells at day 0.

22 (MODS: sepsis; MV; PELOD-11).

Clinical Staph infection; 25% burns that developed sepsis; 58 units/L clinical aspartate aminotransferase, 1.1 g/dL albumin, 38*109/L white blood cells, Halothece mapped Cyanobacteria that can have cyanotoxins that yield similar outcomes observed in patient namely hepatocyte necrosis (11, 21); known 1,468 day 0 elevated transcripts [GO renal allograft rejection (49)]; high level of M0 macrophages at day 3.

23 (MODS: CO toxicity; MV; PELOD-10).

Clinical rhino-/enterovirus; high mapping of Acinetobacter and Burkholderia infections on day 0 and 3; 44 transcripts low at day 0 (GO: histone H3-R26 citrullination). H3-R26 citrullination is involved in increasing immune response (42); 180 transcripts increased from day 0 to 8 (GO: immune system process).

24 (ECMO: hemophagocytic lymphohistiocytosis/EBV; L, R, MV; PELOD-10).

Hemophagocytic lymphohistiocytosis (HLH), EBV detected by PCR and seen in RNA-Seq of patient at day 0; elevation of Bordetella parapertussis at day 0 and 8; enriched for genes involved in T cell activation; highest PERSEVERE gene elevation at time point 0 in cohort; highest level of CD8 T cells and memory CD4 T cells; elevation of mast cell genes at day 3; elevation of blood ALB transcripts at day 0 with five additional liver genes; high level of NAT8 at day 0 associated with kidney tubule cells; patient case write-up described below.

25 (MODS: sepsis; L, R, MV; PELOD-22).

Clinical hemophagocytosis in bone marrow, clinical rhino-/enterovirus; high mapping of Methylobacterium (25), Streptomyces (26) on day 0 and 3; high mapping of Bordetella parapertussis (57) on day 0; 201 transcripts low on day 0 (GO: myeloid leukocyte activation); 265 transcripts increased from day 0 to 8 (GO: translation, cotranslational protein targeting to membrane, GTP hydrolysis and joining of the 60S ribosomal subunit); elevation of eosinophils at day 0 and 3.

26 (MODS: febrile infection-related epilepsy; L, MV; PELOD-32).

POLQ and SASH3 variants; clinical febrile infection-related epilepsy; liver failure; high mapping of Mycobacterium tuberculosis and Pasteurella multocida; both pathogens have been observed to cause seizures (34, 56); day 3 elevated monocytes.

27 (ECMO: sepsis/myocarditis; L, R, MV, C; PELOD-12).

Died; clinical myocarditis and culture-negative endocarditis; clinical Staph infection; liver failure; renal failure; brain injury; respiratory failure; high mapping of Pasteurella multocida (common infection from cat scratches) and Gordonia polyisoprenivorans VH2, P. multocida has been observed to cause endocarditis as well as respiratory tract symptoms (8, 64), and G. polyisoprenivorans VH2 has been observed to cause endocarditis (54); family had two cats as pets; day 3 elevated monocytes; ALB RNA detected in blood at day 0 with five additional liver genes; several pituitary associated genes elevated.

A Case of Virally Induced Genetics, Patient 24

Patient 24 was suggested to have unique potentially syndromic biology based on data presented, agreeing with clinicians who ordered rapid whole-exome sequencing. A previously healthy, active 16 yr old female was admitted with fever, malaise weakness, and upper respiratory symptoms that progressed rapidly into multiorgan failure (Fig. 8). The symptoms started about 10 days before hospital presentation and were mainly respiratory. She was initially treated at home with a 10-day course of amoxicillin. Over the next week, she continued to deteriorate, with on-and-off febrile episodes thought to be “flu-like” with significant myalgia, eventually becoming bedridden with altered mental status leading to an emergency department visit. She continued to deteriorate despite cardiac support and required ECMO within 24 h of being admitted to HDVCH PICU. The rapid deterioration was suggestive of viral illness (identifying EBV by PCR), a ferritin level >10,000 ng/L, and a bone marrow exam suggesting hemophagocytic lymphohistiocytosis (HLH), a diagnosis that typically takes several days to make but is visible in day 0 RNA-Seq. The patient samples cluster relative to control and a 2 yr follow-up transcriptome (Fig. 8). A second tool confirmed EBV reads from RNA-Seq, including the gene EBER1 (Fig. 8). Rapid whole-exome sequencing identified compound heterozygous variants in RNASEH2B that the reference laboratory annotated as likely pathogenic for Aicardi-Goutières, c.58 G>C (V20L amino acid, paternal), and a potential splicing variant c.245-8A>G (unknown protein change, maternal). RNASEH2B is associated with Aicardi-Goutières syndrome, an early-onset autoinflammatory disorder that can present as a monogenic form of systemic lupus erythematosus (SLE) (13). Acquired or secondary HLH may occur in the setting of SLE or other autoinflammatory conditions with or without an infectious trigger. The V20L variant, based on our computational tools, has no alteration of the protein (Fig. 8, D and E), suggesting it to be benign, and therefore, the patient does not have homozygous associated Aicardi-Goutières. However, the V20L variant was used to show an allelic imbalance from the RNA-Seq that was pronounced on all but day 0 (Fig. 8), suggesting that a unique form of Aicardi-Goutières may be present.

Fig. 8.

Patient with active Epstein-Barr virus (EBV) and RNASEH2B splice variant. A: patient timeline: * marks multi-omic data collection points, red marks interventions, black marks major events, and blue marks phenotypes observed. B: clustering of control (red) and patient samples (green), including a 22 mo follow-up. C: second viral read caller to confirm EBV infection with EBER1 (gray) and EBER2 (red) reads from EBV detected at day 0 shown below. D: known structure of RNASEH2A (gray), RNASEH2C (magenta), and RNASEH2B (cyan). E: zoom-in view of V20 not falling at a critical site suggesting L would be functional. F: percent of reads from the RNASEH2B alleles at multiple RNA-Seq measurements suggesting allelic bias to the allele inherited from the father. G: reads aligned near the potential splice site variant inherited from mom that confirms a frameshift and early truncation. H: total RNASEH2B levels over the multiple days suggesting reduced levels relative to controls (red). I: pathways activated in the patient’s 2 yr follow-up, suggesting immune system dysfunction. J: immune system repertoire analysis performed on the patient follow-up for seven chains (TRB, T cell receptor beta; TRA, T cell receptor alpha; IGH, immunoglobulin heavy chain; TRD, T cell receptor delta; TRG, T cell receptor gamma; IGK, immunoglobulin light chain kappa; IGL, immunoglobulin light chain lambda). Shown is the percent of CDR3 identified in each group with a callout for further dissection of the IGH subgroups. K: diversity index (DI) for each of the sequenced chains. DI corresponds to elevated recombination sites, with lower levels suggesting an enrichment of certain sequences within the chain. The IGL is shown for diversity, suggesting the elevation of epitope recognition. This epitope is seen in 27 of the top 35 CDR3 sequences for IGL with a significant enrichment (8.9e-82) and composing 21.7% of all the CDR3 IGL sequences. CRP, C-reactive protein; ECMO, extracorporeal membrane oxygenation; ER, emergency room; HDVCH, Helen DeVos Children’s Hospital; HLH, hemophagocytic lymphohistiocytosis; PICU, pediatric intensive care unit; WES, whole exome sequencing.

Interestingly, the splice variant inherited maternally (c. 245-8A>G) was expressed only in the day 0 samples, yielding a frame-shift mutation (Fig. 8). Levels of RNASEH2B were around 50% in later days (Fig. 8), suggesting the splice variant had a dominant-negative effect, which would be normally suppressed by nonsense-mediated decay (NMD). This loss of NMD in day 0 sample is interesting and correlates with the increasing copies of serum EBV noted clinically and in the RNA-Seq. Interestingly, it is known that EBV-coded proteins can suppress NMD (29). With the patient’s potential syndromic biology identified, the patient had a 2 yr follow-up transcriptome analysis. This revealed excessive inflammation and continued alterations in the immune system (Fig. 8). The immune system was profiled by sequencing the CDR3 region of seven chains, identifying a high level of IgD (Fig. 8) and low diversity of immunoglobulin kappa and lambda chains with enrichment of a functional motif (Fig. 8). Hyper-IgD syndrome can be caused by mutations in MVK, and this gene was reanalyzed in the exome and showed 100% coverage at a minimum of 10X with no variants identified. This all suggests that the patient has extensive immune system dysfunction that combined with an active EBV infection turning on the suppressed RNASEH2B splice variant to cause a dominant-negative HLH that recovered as the EBV was cleared. This patient represents an example of the complex interplay between environment and genetics that we coin as virally induced dominance and highlights the power and need for precision transcriptomic analysis performed on more patients, with these cases missed in large cohort-style study design.

Data and Materials Availability

All of the mapping data for the RNA-Seq are available for gene and transcript levels in the supplemental Excel file at https://doi.org/10.6084/m9.figshare.12252197.v1.

DISCUSSION

Transcriptomics has been speculated to have significant utility for personalized/precision medicine (16, 37). While fields like oncology have embraced these tools, their value remains undetermined for most other areas in health care. Patients admitted to the intensive care units have altered homeostatic pathways that result from unclear triggers. Transcriptomics has the potential not only to identify the degree of disruption but also to identify etiology in patients that often don’t have the luxury of time for a long diagnostic odyssey. Thus far for PICU patients, strategies have focused on developing larger cohorts based on similarity of immune responses, typically utilizing microarray strategies focusing on those that have sepsis or MODS (2, 63). In these models of research, clusters of genes are examined for similarity of response and highlighted, where the differences between individual patients are either ignored or not studied.

In the clinical arena, multiple aspects of disease etiology and outcomes exist, and as the technologies for understanding sepsis and MODS biomarkers have evolved, there is much space for developing tools that have potential to actually impact therapy (52). A typical cohort of PICU patients presenting MODS, as shown in this paper, is composed of a broad etiology of disease, including diverse environmental contributions and various infections. As the case reviews of these patients show, multiple factors are also identified in many of these patients, from genetics and infection to multiple infections sometimes in each patient. Even within an extensive cohort collection, the combinations that give rise to pathologies in each patient are likely such that no two patients present the same. In this work, we show how the collection of blood samples into the PAXgene tube system, followed by reduced ribosomal RNA-Seq, can elucidate MODS understanding via outlier mapping strategies. Collecting information from transcriptomics reveals a remarkable ability to map pathologies and understand the diverse biological responses seen within each patient in the PICU and gains us more insights as the cohorts of information grow. In the age of expanding sequencing technologies and precision medicine, this work highlights that it is important to consider how each patient is unique and not to lose these data by overly investing into cohort-based patient assessments that tend to obscure individual unique outcomes and mechanisms of pathology. This is further exemplified by viral outbreaks, such as COVID-19, that show selectivity in age- and sex-related vulnerability and inflammatory phenotypes leading to differences in presentation, severity, and management (14, 17). These differences may result from common factors or could be contributed by unique biology that is challenging to capture in large cohorts, especially if there is no in-depth analysis in a patient-by-patient manner. Further work is needed to identify correlations of our precision medicine transcriptomics data to additional comorbidities within our cohort.

One of the most promising parts of the PAXgene collection and its analysis could be in pathogen detection, as the sequencing uncovered microbial/viral RNA that were in the sample and correlate well to clinically confirmed pathogens. This emerging field, labeled metagenomics, is still in its fledgling phase as we are trying to understand its clinical implications (19). The potential for clinical relevance is undeniable, as we showed clinical overlap in multiple patients, not only detecting viruses such as measles and EBV, but also bacteria, confirmed by clinical tests and those that cannot be cultured well in clinical laboratories. The power of such findings is in the potential to guide clinical decisions on day 1 as shown by others through discovering unexpected infections (1) and antimicrobial choice (19).

As we sequence more patients, we will have more cases that align phenotype with correlating transcriptomic insights. This will require blood draws from larger numbers of patients collected by similar protocols, but with a focus on understanding each patient instead of the entire cohort. Investments into larger data sets should lead to the scaling of bioinformatic tools that will delineate initiation events, patient responses, and determine when and why patients respond adversely. The current computational workflows have significantly decreased in compute time with the quasi-based alignment strategies, ranging from a few minutes to an hour, depending on available computing resources. One limiting factor that needs to be surmounted is the current bottleneck in the time it takes to generate RNA-Seq data by existing platforms. Newer tools such as Nanopore hold a promise for speed, for clinical care to transition from observational science, as in this study, into a tool that generates personalized/precision data in real time, leading to the appropriate confirmatory clinical tests to guide decision-making in the PICU.

GRANTS

Spectrum Health Office of Research funding initiative for precision medicine (S. Rajasekaran, M. L. Leimanis) and National Institutes of Health Office of the Director and National Institute of Environmental Health Sciences Grant K01ES025435 (J. W. Prokop).

DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

AUTHOR CONTRIBUTIONS

J.W.P., M.L.L., M.A., A.S.B., J.H., D.J.S., C. Bupp, and S.R. conceived and designed research; J.W.P., M.W., M.A., M.R., W.P., M.L.B.-S., D.C.S., M.D., B.E.B., M.E., J.H., and S.R. performed experiments; J.W.P., R.S., R.G., D.N., K.U., N.L.H., J.V.V., J.S.S., O.S., T.A.L., C.C.K., V.S., J.M., and S.R. analyzed data; J.W.P., R.S., R.G., M.L.L., B.C., N.L.H., B.B., A.H., C. Braunreiter, A.S.B., M.H., C.M., D.J.S., L.H.S., C. Bupp, and S.R. interpreted results of experiments; J.W.P., R.G., and S.R. prepared figures; J.W.P., D.N., M.A., C. Bupp, and S.R. drafted manuscript; J.W.P., R.S., R.G., M.L.L., D.N., K.U., N.L.H., L.H.S., C. Bupp, and S.R. edited and revised manuscript; J.W.P., R.S., R.G., M.L.L., D.N., K.U., B.C., N.L.H., J.V.V., J.S.S., O.S., T.A.L., B.B., A.H., C. Braunreiter, C.C.K., V.S., J.M., M.W., M.A., M.R., A.S.B., W.P., M.L.B.-S., D.C.S., M.D., B.E.B., M.E., J.H., M.H., C.M., D.J.S., L.H.S., C. Bupp, and S.R. approved final version of manuscript.