{'fraction_non_alphanumeric': 0.029998240572400445, 'fraction_numerical': 0.015072429769514985, 'mean_word_length': 5.139153915391539, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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 RNAgenerated 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.

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.

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 Network-Analyst3.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.

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 hospi-tals 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).

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. 2A). 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. 2B).
Validated biomarkers for pediatric sepsis severity (60) were examined and found in two samples to have significant elevation ( Fig. 2C), 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. 2D). 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 (

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).

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. 4A). Patients 5 and 11 had the highest bacterial read mapping.
For bacteria (Fig. 4B), 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. 4C). 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. 5A). 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. 5B). Interestingly, both of these patients had a day 3 spike in AST blood levels, not at day 0 when ALB is high (Fig. 5C). 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. 5D).
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. 5E). Three patients are identified as extreme  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 GEN-CODE 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

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   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.

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. 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. 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. 8A). 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. 8B). A second tool confirmed EBV reads from RNA-Seq, including the gene EBER1 (Fig. 8C). 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. 8F), suggesting that a unique form of Aicardi-Goutières may be present.
Interestingly, the splice variant inherited maternally (c. 245-8AϾG) was expressed only in the day 0 samples, yielding a frame-shift mutation (Fig. 8G). Levels of RNASEH2B were around 50% in later days (Fig. 8H), 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. 8I). The immune system was profiled by sequencing the CDR3 region of seven chains, identifying a high level of IgD (Fig.  8J) and low diversity of immunoglobulin kappa and lambda chains with enrichment of a functional motif (Fig. 8K). 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 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. 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.