
==== Front
Epigenomics
Epigenomics
Epigenomics
1750-1911
1750-192X
Taylor & Francis

39225561
10.1080/17501911.2024.2375187
2375187
Version of Record
Research Article
Research Article
Integrated epigenomic exposure signature discovery
Schuetter Jared a
Minard-Smith Angela a
Hill Brandon b
Beare Jennifer L a
Vornholt Alexandria d
Burke Thomas W g
https://orcid.org/0000-0002-7368-7742
Murugan Vel c
Smith Anthony K a
Chandrasekaran Thiruppavai c
Shamma Hiba J a
Kahaian Sarah C a
Fillinger Keegan L a
Amper Mary Anne S d
Cheng Wan-Sze d
Ge Yongchao d
George Mary Catherine d
Guevara Kristy d
Lovette-Okwara Nora d
Mahajan Avinash d
Marjanovic Nada d
Mendelev Natalia d
Fowler Vance G g
McClain Micah T g
Miller Clare M d
Mofsowitz Sagie d
Nair Venugopalan D d
Nudelman German d
Evans Thomas G h
Castellino Flora j
Ramos Irene d
Rirak Stas d
Ruf-Zamojski Frederique d
Seenarine Nitish d
Soares-Shanoski Alessandra d
Vangeti Sindhu d
Vasoya Mital d
Yu Xuechen d
Zaslavsky Elena d
Ndhlovu Lishomwa C e
Corley Michael J e
Bowler Scott e
Deeks Steven G f
Letizia Andrew G i
Sealfon Stuart C d
Woods Christopher W g
Spurbeck Rachel R * a
a Health Business Unit, Battelle Memorial Institute, Columbus, OH 43201, USA
b Nationwide Insurance, Columbus, OH 43215, USA
c Center for Personalized Diagnostics, Biodesign Institute at Arizona State University, Tempe, AZ 85281, 85281USA
d Icahn School of Medicine at Mount Sinai, New York City, NY 10029, 10029USA
e Division of Infectious Diseases, Department of Medicine, Weill Cornell Medicine, New York, NY 10021, USA 
f University of California San Francisco, San Francisco, CA 94143, 94143USA
g Division of Infectious Diseases, Duke University, Durham, NC 27710, USA 
h Barinthus Biotherapeutics, Harwell, Oxfordshire, UK
i Naval Medical Research Unit INDO PACIFIC, Singapore
j Biomedical Advanced Research & Development Authority-Administration for Strategic Preparedness & Response,Washington, DC 20201, USA
* CONTACT: spurbeck@battelle.org
3 9 2024
2024
3 9 2024
16 14 10131029
Aptara29 6 2024
01 8 2024
15 1 2024
28 6 2024
© 2024 Battelle Memorial Institute. Published by Informa UK Limited, trading as Taylor & Francis Group
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivatives License (http://creativecommons.org/licenses/by-nc-nd/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited, and is not altered, transformed, or built upon in any way. The terms on which this article has been published allow the posting of the Accepted Manuscript in a repository by the author(s) or with their consent.

Aim: The epigenome influences gene regulation and phenotypes in response to exposures. Epigenome assessment can determine exposure history aiding in diagnosis.

Materials & methods: Here we developed and implemented a machine learning algorithm, the exposure signature discovery algorithm (ESDA), to identify the most important features present in multiple epigenomic and transcriptomic datasets to produce an integrated exposure signature (ES).

Results: Signatures were developed for seven exposures including Staphylococcus aureus, human immunodeficiency virus, SARS-CoV-2, influenza A (H3N2) virus and Bacillus anthracis vaccinations. ESs differed in the assays and features selected and predictive value.

Conclusion: Integrated ESs can potentially be utilized for diagnosis or forensic attribution. The ESDA identifies the most distinguishing features enabling diagnostic panel development for future precision health deployment.

Plain Language Summary

This article introduces ESDA, a new analytic tool for integrating multiple data types to identify the most distinguishing features following an exposure. Using the ESDA, we were able to identify signatures of infectious diseases. The results of the study indicate that integration of multiple types of large datasets can be used to identify distinguishing features for infectious diseases. Understanding the changes from different exposures will enable development of diagnostic tests for infectious diseases that target responses from the patient. Using the ESDA, we will be able to build a database of human response signatures to different infections and simplify diagnostic testing in the future.

Tweetable Abstract

The exposure signature discovery algorithm simplifies diagnostic target discovery from integrating epigenomic and transcriptomic data sets enabling rapid diagnostics based on host response to exposure.

Article highlights

Integration and parsing of multi-omic data generates precise exposure signatures.

The number of features in the ESDA model is not strongly correlated with signature performance, but depends on the underlying biology.

The low sample size for some exposure types contributes to instability in some model predictions (BA), and therefore, high variability in accuracy metrics during LOOCV.

Ensemble results are robust for most exposure types.

The most prominent combination of epigenetic methods deployed in the ESDA models were EPIC array for DNA methylation, RNAseq, H3K4me1 and H3K4me3 histone modification assays.

Explosives (PETN) imprint the epigenome and result in a signature that can be detected 6–8 weeks after exposure.

ESDA-developed SARS-CoV-2 signatures track with acute SARS-CoV-2 with 89% performance. Individuals with convalescent SARS-CoV-2 have a distinct signature that indicates immune epigenetic remodeling.

ESDA-derived S. aureus signature can discriminate robustly between MRSA and MSSA in S. aureus infections.

An exposure signature database can enable a universal diagnostic test for exposure health that can distinguish infectious diseases and chemical exposures.

Keywords: 

diagnostics
epigenomics
exposure health
infection
machine learning
multi-omics
transcriptomics
Defense Health Agency 10.13039/100009898 9700140 Defense Advanced Research Programs Agency Biological Technology Office W911NF19C0041 This work was funded by the Defense Advanced Research Programs Agency Biological Technologies Office contract number W911NF19C0041 (RRS) and Defense Health Agency grant 9700140 through the Naval Medical Research Center (AGL).
==== Body
pmc1. Introduction

Traditional medical diagnosis depends on direct testing for the pathogen that is causing the patient's symptoms. However, testing for individual pathogens requires methods that are specific to each causative agent and often are fleeting targets. In time-critical cases, it is essential to identify the correct cause of disease to provide lifesaving treatment; however, tests for particular pathogens may not be able to determine the cause due to low detection limits or a transient presence that initiated the disorder. For example, the proper treatment of systemic inflammatory response syndrome (SIRS) depends on whether it is caused by an infection (sepsis); antibiotics are effective for microbial caused sepsis but have no effect on non-infectious SIRS [1]. For this reason, rapid differentiation of the cause of SIRS is necessary to enable an accurate and timely treatment regimen. Sepsis is one of the major causes of morbidity and early diagnosis is key to survival [2]. Blood cultures used to identify the bacterial etiology of sepsis take over 24 h and have a high incidence of false negatives, and molecular tests for specific pathogens are not sensitive enough to diagnose bloodstream infections. Another example relates to lower respiratory tract infections (LRTI), which have a high morbidity [3]. Since non-infectious respiratory syndromes often look symptomatically like LRTI, this can confound diagnosis when methods are used that search for specific pathogens. To improve recovery rates, it is therefore necessary to develop more sensitive and universal diagnostic tests for rapid identification and treatment of both infectious and non-infectious diseases which can be accomplished by targeting diagnostic testing to the host response.

Focusing on the host's response to the exposure provides the opportunity for a universal diagnostic by leveraging commonalities in each exposure response. While each infection stimulates unique responses, some factors, such as immune system responses, can be generalized to an infection by a virus versus a bacterium or a non-infectious disease. In all cases, whether symptoms are caused by an autoimmune deficiency or an infection, there is a human physiological response due to the combination of unique and universal changes in molecular regulatory mechanisms. Environmental stimuli can elicit a response through changes in methylation, regulatory RNAs, chromatin accessibility, and/or histone modifications that leads to changes in gene expression or post translational modification of proteins. By focusing on the effects of environmental changes on these regulatory mechanisms, a universal diagnostic platform can be developed focusing on a small, targeted panel of epigenomic features for distinguishing exposures that impact the epigenome. A targeted panel of identified host response features to exposure potentially could be developed into a rapid assay utilizing amplicon-based sequencing on an instrument such as the MinION, providing data in hours.

Building a host based diagnostic platform requires new methods to analyze and integrate data from large data sets. Here we developed and implemented a machine learning algorithm, the exposure signature discovery algorithm (ESDA), to identify the most important features present in multiple epigenomic and transcriptomic dataset to produce an integrated exposure signature (ES). Conventional research studies focus on the response of individual regulatory mechanisms to exposure, providing extensive data on differential exposure profiles, seldom attempting to integrate results from multi-omic data or distilling this information down to the most highly relevant features for distinguishing different exposures by a diagnostic test that can be conducted in a healthcare setting. The ESDA can used for diagnostic marker discovery to derive diagnostics that target these ESs using more common diagnostic methods such as quantitative polymerase chain reaction (qPCR) or targeted assays that are much more rapid.

The ESDA could also be used to reduce the time to answer for diagnosing rare infections or an autoimmune disorder, where the diagnostic process usually takes multiple assays. In these cases, the physician orders a set of tests, and if the results do not identify the cause, another set of tests is ordered. This process can be arduous with a high impact on the patient, as they are not being treated for the correct issue during the diagnostic period. A large omics-based assay analyzed by the ESDA can reduce the time to diagnosis in these cases by capturing and analyzing data to get to the answer from one assay. Furthermore, evolving techniques such as the MinION, which can produce and analyze data in real time, are rapidly reducing the sample to answer time for large omics experiments, making it more comparable to more traditional methods. Currently, multi-omic data analysis provides too much data for use as a clinical infectious disease diagnostic test, resulting in long sample to answer times and large computational effort. Here, the ESDA parses millions of features derived from multiple data sets and distills this information down to the most important features, ranging from two features to a few hundred. With the recent shift in the scientific field toward multi-omic studies, ESDA will output combined, holistic ESs focusing on the most important features for differentiation to enable development of host-based diagnostics.

2. Materials & methods

2.1. Human participants

The authors state that they have obtained appropriate institutional review board approval for all human subject studies outlined in this work. In addition, informed consent was obtained from the participants involved. The SARS-CoV-2 study protocol was approved by the Naval Medical Research Center institutional review board (protocol number NMRC.2020.0006) in compliance with all applicable federal regulations governing the protection of human subjects. The S. aureus sepsis study protocol was reviewed and approved by the Duke Medical School institutional review board (protocol number Pro00102421). The Bacillus anthracis vaccination study protocol was approved by the Battelle Memorial Institute institutional review board (protocol number 0708-100130312). HIV samples were derived from the iPrEX [4] and SCOPE studies (https://hividgm.ucsf.edu/scope-study). Influenza (H3N2) samples were derived from the Vaccitech H3N2 influenza vaccine study which was published [5]. Subjects provided informed consent prior to participation.

2.2. SARS-CoV-2 sample collection

PBMC samples were obtained from the SARS-CoV-2 Health Action Response for Marines (CHARM) cohort study, which has been previously described in Chen et al. 2023 [6] and Letizia et al. 2021 [7]. Samples analyzed were from participants who tested negative at least once 14 days prior to testing positive for SARS-CoV-2 (prediagnosis, N = 101), on the day of the first positive qPCR test (acute, N = 72) and tested negative at least 14 days after the last positive test (convalescent, N = 104).

2.3. S. aureus sepsis sample collection

These samples were obtained from patient donors who tested positive for either Methicillin-resistant Staphylococcus aureus (MRSA, N = 13) or Methicillin-sensitive Staphylococcus aureus (MSSA, N = 7). RNA and PBMC collection was previously described in Chen et al., 2023 [6]. Controls were PBMC samples obtained from a commercial vendor (Dx Biosamples LLC, N = 12).

2.4. HIV sample collections

Two studies were utilized in this work. The first was a retrospective analysis of PBMCs isolated on the iPrEx Study described in Grant et al. 2010 [4]. The iPrEx cohort (Iniciativa Profilaxis Pre-Exposición, or the PreP pre-exposure prophylaxis initiative) was a phase III clinical trial aimed at determining the efficacy of pre-exposure prophylaxis (PrEP), an antiretroviral prophylactic for preventing HIV infection. Ten iPrEx samples were analyzed from donors approximately ∼250 days prior to testing positive for HIV-1 (while HIV negative, ‘pre’), on the day of the HIV-1 positive blood draw (‘acute’, prior to onset of antiretroviral treatment) and approximately ∼200 days after treatment (‘chronic’, while on antiretroviral treatment). The second sample collection was derived from PBMCs collected and isolated on the SCOPE study (https://hividgm.ucsf.edu/scope-study).

2.5. Influenza A (H3N2) sample collection

Pre- and Post-challenge PBMC samples were analyzed from the BARDA-Vaccitech FLU010 Study with A/BELGIUM/4217/2015 (H3N2) challenge aimed to assess the protective capabilities of a vaccine candidate, VTP-100, as a standalone influenza vaccine [5]. Thirteen samples were analyzed from donors who received the placebo vaccine and were sampled before challenge dose and again 28 days after the challenge dose.

2.6. B. anthracis vaccination sample collection

Participants were individuals who were vaccinated against anthrax and worked with Bacillus anthracis in a controlled facility while wearing personal protective equipment. These individuals were trained scientists working in Biosafety Level 3 (BSL3) facilities and received vaccination against B. anthracis infection as part of the safety protocols of the facility. An inactivated, acellular vaccine primarily containing the non-pathogenic protective antigen (PA) protein, called BioThrax® (tradename) or Anthrax Vaccine Adsorbed was utilized. The donors also handled chemical weapon agents, precursor molecules to chemical agents, biological warfare agents, explosives, precursors to explosives, agricultural chemicals, or radioactive materials in a controlled environment. Blood collections were performed on participants after informed consent was given and a metadata questionnaire was completed. Controls were also collected for this cohort from individuals who worked in the same facility who were not vaccinated or worked with chemical agents, biological warfare agents, explosives, precursors to explosives, agricultural chemicals, or radioactive materials. From each participant, approximately 2.5 ml blood was collected in one BD PAXgene™ Blood RNA Tube (Fisher Scientific, catalog number 14-959-51D) and approximately 8 ml blood was collected in each of four BD Vacutainer™ Glass Mononuclear Cell Preparation (CPT) Tubes (32 ml total, Fisher Scientific, catalog number 14-959-51D). Within 90 min of collection, the PAXgene RNA tubes were placed into a -20°C freezer and the CPT tubes were processed for PBMCs and plasma. Nineteen B. anthracis vaccinated samples and 35 controls were assessed.

2.7. PBMC isolation from B. anthracis vaccination sample collection

CPT tubes were gently inverted and centrifuged at 1,500 × g for 10 min at room temperature (RT). After centrifugation, aliquots of plasma were taken from the top layer and stored at -80°C. For each sample, the PBMC layers were collected and combined in two (2) sterile 15 ml conical tubes. The combined PBMCs were then centrifuged at 800 × g for 3 min at RT. After centrifugation, the supernatant was discarded, the cells washed with 1 ml Dulbecco's Phosphate Buffered Saline (DPBS) and then combined into one (1) 15 ml conical tube. The cells were centrifuged again 800 × g for 3 min at RT and the supernatant discarded. The cells were then resuspended in Bambanker's medium to a final concentration of 1 × 106–2.5 × 106 cells/ml and split into 1 ml aliquots for freezing. The vials were placed into a Mr. Frosty filled with 100% isopropyl alcohol and stored at -80°C or in a pelican cooler full of dry ice overnight before transferring into liquid nitrogen for long term storage.

2.8. RNAseq Data Generation for B. anthracis vaccination & HIV (SCOPE dataset)

RNA sequencing was conducted from PAXgene RNA tubes for B. anthracis vaccination or PBMCs from the HIV studies at Ichan School of Medicine Mount Sinai. Briefly, total RNA was extracted from PAXgene blood or PBMCs using the RNAdvance Blood Kit and RNAdvance cell kit, respectively (Beckman Coulter) on a BioMek FXP Laboratory Automation Workstation (Beckman Coulter). Concentration and integrity (RIN) of isolated RNA were determined and RNAseq libraries constructed from total RNA using the Universal Plus mRNA-Seq kit (#9133, Tecan Genomics, San Carlos, CA, United States) as described in Chen et al. [6].

2.9. RNAseq data generation for influenza, HIV (iPrEx study), MRSA & MSSA

Total RNA extraction was performed from 100K PBMC cells using the RNAeasy Mini Kit (QIAGEN) according to manufacturer instructions. RNA quality and yield was assessed using TapeStation with RNA ScreenTape (Agilent) and Qubit RNA HS assays (Invitrogen). Stranded, poly-A-selected mRNA sequencing libraries were generated using KAPA mRNA HyperPrep Sequencing Library Kit (Roche Molecular Systems) using 100 ng total RNA input per sample. Sequencing was performed on an Illumina NovaSeq 6000 instrument with S4 flow cell using 100 bp paired-reads, and ∼50 M read clusters per sample.

2.10. miRNA data generation for HIV (iPrEx study), MRSA, MSSA & influenza

Total RNA was extracted from PBMC cells using Norgen total RNA (including small/micro RNAs) extraction kit (Norgen) according to manufacturer's instructions. Libraries were prepared using QIAseq miRNAseq Library Kit (QIAGEN) with 100 ng total RNA input. Sequencing was performed using Illumina NextSeq 500 High-Output 75 bp single-end read (SR), with 5–10 M read clusters per sample.

2.11. Bulk ATAC data generation for HIV (iPrEx study), MRSA, MSSA, & influenza

Frozen PBMC cryovials were thawed following standard procedures and counted. 50,000 PBMCs were washed with 50 μl cold 1X phosphate buffered saline (PBS). After centrifugation, the pellet was resuspended in 50 μl of cold Lysis Buffer (48.5 μl Resuspension Buffer, 0.5 μl 10% NP-40 (final 0.1% v/v), 0.5 μl 10% Tween-20 (final 0.1% v/v), 0.5 μl 1% Digitonin (final 0.01% v/v)), pipetted up and down 3× gently to resuspend the cells and incubated on ice for 3 min. 1 ml of wash buffer (990 μl Resuspension Buffer and 10 μl 10% Tween-20 (final 0.1% v/v)) was added and centrifuged at 500 × g for 10 min at 4°C. 50 μl of transposition mix (25 μl 2X TD Buffer (Tagment DNA Buffer), 16.5 μl 1X PBS, 0.5 μl 10% Tween-20 (final 0.1% v/v), 0.5 μl 1% Digitonin (final 0.01% v/v), 2.5 μl Tn5 Transposase (Tagment DNA Enzyme 1) and 5 μl nuclease-free H2O) was added to the pellet, pipetted up and down 6× gently to resuspend the nuclei and incubated at 37°C for 30 min in a thermomixer at 1,000 RPM. DNA was isolated using Qiagen MinElute Reaction Cleanup Kit and eluted in 10 μl of Elution Buffer. To prepare ATAC sequencing libraries, the following were combined: 10 μl purified transposed DNA, 10 μl nuclease-free H2O, 2.5 μl Ad1 indexing primer (25 μM), 2.5 μl Ad2 indexing primer (25 μM), 25 μl NEBNext High-Fidelity 2X PCR Master Mix. The samples were amplified in a thermal cycler with the following program: 72°C 5 min, 98°C 30 s, [98°C 10 s, 63°C 30 s and 72°C 1 min for 5 cycles]. At the end of the amplification, 5 μl of each partially-amplified library was quantified by SYBR green quantitative PCR (qPCR) to determine the number of additional PCR cycles needed. Each sample was then amplified by the number of additional thermal cycler cycles needed to reach 1/3 of the maximum. A double-sided bead purification was performed and the quality and quantity assessed by Qubit and Bioanalyzer HS DNA kit for each library. Since individual libraries were dual indexed during the PCR amplification step with unique 8 bp-index sets (using a combination of Ad1 and Ad2 indexing primers), the libraries were pooled for sequencing and sequenced 2 × 100 bp in a NovaSeq to a minimal depth of 30 million paired end reads per library at New York Genome Center.

2.12. EPIC methylation data generation for SARS-CoV-2, influenza, MRSA & MSSA

All samples were stored and processed for methylation analysis by Illumina Infinium Human Methylation EPIC Bead Chip array (Illumina Inc., CA, USA) as described in [6]. To avoid method induced effects, samples from the same individual were analyzed in the same batch.

2.13. MNase digestion & ChIPmentation Data Generation for HIV (iPrEx study), MRSA & MSSA

MNase digestion of chromatin and ChIPmentation was performed as per Schmidl et al. [8] with following modifications. Frozen PBMCs were centrifuged twice at 600 × g for 10 min at 20°C and washed with RPMI+ 20% FBS. The cell pellet was resuspended in PBS to achieve a final volume of 7.5 × 106 cells/ml. An equal volume of the MNase digestion buffer (0.02 U/μl MNase (Thermo Scientific, Cat. No: 88216) in 2× lysis buffer) was added to the cell suspension, then incubated on ice for 20 min, followed by incubation at 55°C water bath for 10 min. The digestion was terminated using 30 mM EGTA and the DNA was quantified using Qubit™ 4 fluorometer using Qubit DNA quantification kit (Q22231, ThermoFisher Scientific, USA). 500 ng of MNase treated chromatin was incubated with 50 μl of Pierce™ protein A/G magnetic beads (Cat# 88802, ThermoFisher Scientific, USA) with equal volume of RIPA-LS buffer (0.1% sodium deoxycholate (DOC), 0.1% SDS, 1% Triton-X-100, 10 mM Tris-HCL, 1 mM EDTA, 140 mM NaCl) and appropriate histone modification specific antibodies (H3K27me3, Cat# 31618020; H3K4me1, Cat# 2119002; H3K4me3, Cat# 22118006; H3K9me3, Cat# 9919003; H3K27ac, Cat# 9919002; H3K36me3, Cat# 11518001, all from Active Motif, Carlsbad, USA) at 4°C overnight at a rotating shaker. Chromatin absorbed bead was subjected to successive washes, two-times each in RIPA-LS, RIPA-HS (0.1% sodium deoxycholate (DOC), 0.1% SDS, 1% Triton-X-100, 10 mM Tris-HCL, 1 mM EDTA and 360 mM NaCl), RIPA-LiCl (250 mM LiCl, 0.5% NP40, 0.5% SDS, 10 mM Tris-HCL and 1 mM EDTA) and resuspended in 24 μl of 1X tagmentation buffer and 1 μl of tagmentation enzyme at 55°C for 10 min. Tagmented chromatin was washed two-times in RIPA-LS buffer and eluted with elution buffer (10 mM Tris, 0.1% SDS and 300 mM NaCl) with 1 mg/ml proteinase K at 55°C for 60 min. DNA was size selected using KAPA Pure Beads and Illumina sequencing index was added to the DNA using KAPA HiFi HotStart Ready Mix (Cat #KK2601, Roche Diagnostics, IN, USA). Sequencing was performed at the Arizona State University NGS core facility.

2.14. Multiplexed indexed T7 ChIP-seq (Mint-ChIP) data generation for influenza

Histone modification sites were identified using the Multiplexed Indexed T7 Chromatin Immunoprecipitation (Mint-ChIP) assay described by van Galen et al. [9] and an updated version 3 of the Mint-ChIP protocol [10]. Cryopreserved PBMCs were thawed on ice and cell number and viability were assessed using microscope and hemocytometer. Approximately 50 k to 100 k viable cells were used from each sample, or ∼7 k to 15 k cells per antibody. Briefly, cells were lysed for 60 min at 4°C and nuclei were digested with 300 units of micrococcal nuclease (MNase, New England Biolabs) for 10 min at 37°C. Annealed, double-stranded T7 adapters with unique sample barcodes were then ligated to chromatin DNA fragments for 2 h at 25°C. Samples were pooled and split for incubation with histone modification antibodies (Histone H3 (H3) antibody (#39763, Active Motif), Histone H3 lysine 27 acetylation (H3K27ac) antibody (#39133, Active Motif), Histone H3 lysine 36 trimethylation (H3K36me3) antibody (#61101, Active Motif), Histone H3 lysine 4 trimethylation (H3K4me3) antibody (#ab8580, Abcam), Histone H3 lysine 4 monomethylation (H3K4me1) antibody (#ab8895, Abcam), Histone H3 lysine 27 trimethylation (H3K27me3) antibody (#07-449, Millipore), Histone H3 lysine 9 trimethylation (H3K9me3) antibody (#ab8898, Abcam)) at 4°C overnight. Chromatin fragments bound to each antibody were isolated using Protein G Dynabeads (Thermo Fisher) and digested with proteinase K for 1 h at 63°C. DNA was purified using 1X SPRIselect bead clean-up (Beckman Coulter, retains fragments >100 bp). In vitro transcription was performed using HiScribe T7 RNA polymerase (New England BioLabs). RNA was purified using MyOne Silane Dynabeads (Thermo Fisher) followed by reverse transcription using SBS3-random hexamer primer with Superscript III enzyme (Thermo Fisher). cDNA was purified using SPRIselect beads and used for library preparation PCR using Terra DNA polymerase (Takara) with P7_bc primer (forward with antibody-specific barcode) and P5_SBS3 (reverse) primer. Double size selection of libraries was performed using SPRIselect beads (selects >200 bp and <700 bp fragments). Libraries were sequenced using NovaSeq 6000 with S4 flow cell in 150-bp PE configuration.

2.15. Single cell RNA sequencing (scRNAseq) – B. anthracis vaccination & SARS-CoV-2

ScRNAseq was performed (10× Genomics, Pleasanton, CA) and quality controlled using a Bioanalyzer prior to sequencing as described in [6].

2.16. Single Nucleus (sn) ATAC sequencing (snATACseq, separate) data generation for HIV, SARS-CoV-2, & influenza

Frozen PBMC vials were thawed and nuclei were isolated. Nuclei concentration was assessed using a Nexcelom Cellaca fluorescence cell counter. SnATACseq following the Chromium Single Cell ATAC Reagent Kits V1.1 User Guide (10× Genomics, Pleasanton, CA)was performed on 10,000 targeted nuclei for good quality samples and fewer when the sample quality was lower with 12 PCR cycles in the barcoding emulsion. Libraries were quantified by Qubit 3 fluorometer (Invitrogen) and quality was assessed by Bioanalyzer (Agilent). Equivalent molar concentrations of the libraries were pooled and the reads were adjusted after sequencing the pools on a Miseq (Illumina). Libraries were then sequenced on a NovaSeq 6000 (Illumina) following 10X Genomics recommendations. The sequencing parameters were the following:Sn ATACseq: Read 1 = 100; Read 2 = 100; I7 index = 8; I5 index = 24

Sn RNAseq: Read 1 = 100; Read 2 = 100; I7 index = 10; I5 index = 10

2.17. Epigenetic signature discovery algorithm

The purpose of the epigenetic signature discovery algorithm (ESDA) is to allow for the automated identification of multi-modal epigenetic biomarkers associated with exposure to biological and chemical sources of concern. The methodology used within the algorithm is also intentionally designed to avoid reliance on annotated regions of the epigenome in public databases, as these annotations are typically skewed by findings from well-funded and often-studied domains, creating a potential bias toward discovering such markers if one is not careful. Instead, the approach is data-driven in the identification of the markers, after which a secondary analysis is needed to link those markers back to annotated (or new) biological functions for interpretation.

The ESDA (Figure 1A–D) begins with PBMC samples drawn from populations that represent multiple exposure states. The typical use case is a two-state model for ‘exposed’ and ‘not exposed’, but the states could also represent ‘exposure Type A’ and ‘exposure Type B’ or ‘acute exposure’ and ‘chronic exposure’. The ESDA can accommodate more than two exposure states as well to find sets of biomarkers that collectively discriminate between multiple types of exposures at once. The next step in the process is to analyze the samples using one or more epigenetic assays to generate sets of features that are used to train assay-level models to distinguish between the exposure states. Each of these models captures the relationship between the exposure states and the sample features generated within a specific assay, allowing for the identification of assay-specific markers that could be used for diagnosis or serve as therapeutic targets. The next stage of the pipeline is to combine those assay-level models into a single exposure-level model. This exposure-level model relates all the features across the assays to the exposure states, and in an iterative process removes from consideration features that are uninformative or redundant with other features. By examining the remaining features, one can then ascertain which of the features are providing the discriminatory power and are good candidates for secondary investigation for biological relevance.

Figure 1. The Epigenetic Signature Discovery Algorithm (ESDA) comprises models built at the assay level and the exposure level. (A) An overview of the ESDA model. The goal of the algorithm is to find epigenetic features across the assays associated with one or more types of exposure, including non-exposure (i.e., healthy controls). (B) To build assay-level models within the ESDA, data are first normalized by genome position, either through binning (for seq-based assays) or by locus (for methylation panels) (1). Optionally, the features can be ranked and reduced through a statistical test (e.g., ANOVA p-value or fold-change criterion) (2). Finally, a recursive feature elimination procedure is used to iteratively reduce the set of features to those that are most predictive of the exposure outcomes of interest (3). (C) To combine assay-level models (1), the assay-level features are combined into an overall profile (2). A Sparse-Group LASSO model then combines those features, with the assay as the group label and selects the optimal feature set that is both parsimonious and predictive of the exposure outcomes (3). (D) When all of the assay-level features are not measured for some of the samples (1), an ensemble-based exposure model can be used, where exposure predictions are made at the assay level, then the posteriors are combined with weights proportional to the cross-validated accuracy of the constituent models of the ensemble (2).

Note the ESDA, as shown in the study below, is intended to be a screening tool to sift through large numbers of potential features and find the ones that appear to be correlated with an exposure response. However, these scenarios will typically suffer from multiple complications, including low sample sizes relative to the number of features, batch effects when multiple cohorts are involved, epigenetic changes from sources other than the exposure of interest (e.g., age, ethnicity, environment, lifestyle) and unclear connections back to underlying biological mechanisms for regions of the epigenome that are not well studied. The features selected by the ESDA would still need to be examined in a more targeted and controlled study to establish reproducibility of the biomarker and its relationship to physiology, thereby allowing the identification of therapeutic targets or development of a rapid diagnostic assay.

2.18. Description of the ESDA

2.18.1. Data pre-processing & normalization

The ESDA was designed to work with multi-omic assays for chromatin accessibility (e.g., ATAC-seq), transcription (e.g., RNAseq, miRNAseq), histone modification (e.g., MINT-ChIP-seq, ChIPmentation) and methylation (e.g., EPIC). In most cases, these assays rely on sequencing and alignment of reads to the human genome to better understand regulation of coding regions of the DNA or gene expression. EPIC arrays, while not sequencing and alignment based, rely on specific loci in the human genome. The first step of the ESDA (Figure 1B, Step 1) is to represent these data sources in a common format to allow them to be modeled in similar ways. For assays like RNAseq or MINT-ChIP-seq that produce sequence alignments at any position in the genome, this common format is a series of binned sequence counts. Non-overlapping bins of a 500 bp are defined across the positions in the human genome, then each aligned read is mapped to one and only one bin by finding the bin that most overlaps the positions occupied by the read. The features then become the count of reads mapped to each bin, ranging from zero to the maximum reads in the sample, with one feature per bin. For assays like the EPIC Infinium array, which focuses on specific loci within the genome, the typical representation is as a series of proportions (e.g., beta values) at those loci. In methylation data, for example, these proportions represent the fraction of reads that aligned at the locus and were methylated. Here, the features range from 0 to 1 and there is one feature per locus.

Since all of the samples are to be combined into a single assay-level model, a crucial step is to standardize the features (Figure 1B, Step 2) to account for variation in the number of reads generated and/or batch effects in the analysis (e.g., different reagent stocks, technicians, protocol variations across labs, etc.). For FASTQ data we apply a quality filter, align the data to the hg38 reference, convert the sorted .bam file to sorted bed format and bin the bed file using BEDOPS [11]. For count-based binned data, we use DESeq2 for normalization. For proportional locus-level data, such as the EPIC methylation array data, the beta values can be processed by the ESDA model.

2.18.2. Assay-level modeling

Once the samples have been standardized, the ESDA uses a recursive feature elimination [12] approach to triage the large number of features simultaneously (typically, the number of features p is much larger than the sample size n) and identify the features that carry the most information to help identify the exposure states for the samples (Figure 1B, Step 3). Recursive feature elimination (RFE) is a procedure that iteratively builds a predictive model, assesses the importance of the predictors in the model, then removes the least important predictors from consideration in future iterations. RFE is model-agnostic, so any kind of classification model can be used in this step. During ESDA development, we found that a ridge classifier worked well in this role, as the shrinkage of the coefficients helped to accentuate the gradient between feature importance and the features could be more easily ranked at each step than in a LASSO model [13] since the coefficients are non-zero and ties do not occur. The RFE procedure produces a succession of models, each of which contains fewer features than the last and one can explore the tradeoff of model prediction accuracy vs. the number of features to settle on a model that is parsimonious for the application and still has reasonable predictive power. The features for this model are candidate assay-level biomarkers associated with the exposure states of interest and can be prioritized for review based on the magnitudes of the coefficients, which provide a natural ranking of their importance. We used a step size of 0.01 for RFE, meaning that the least important 1% of features was removed from consideration at each iteration.

2.18.3. Exposure-level modeling

The assay-level models developed using RFE produce many candidate features (i.e., bins or loci in the genome) that could have potential value in diagnostics or therapeutics. However, since all of the assays are measuring different aspects of the same biological mechanisms, it is likely that many of these features will carry redundant information. To streamline the feature set and arrive at a more parsimonious model, the ESDA builds an exposure-level model using the Sparse-Group LASSO (SGL, [14]) model (1C).

The SGL is similar to a LASSO model [13] in that it is a linear model that can assign coefficients with a value of zero to the features, effectively selecting for features through the non-zero coefficients. The difference is that the SGL contains an additional penalty term in the loss function that enforces a group-level selection as well. In the ESDA, the groups are the assays, so the SGL will simultaneously select the most important assays as well as the most important features within those assays.

The primary requirement for using the SGL is that there is a full set of features for each observation, meaning all of the assays need to be run on the same cohort of individuals. In many cases, due to budget or technology access limitations, this may not be true. In cases where the cohorts vary across one or more assays, the ESDA can also employ an ensemble modeling approach for the exposure-level modeling (1D). The idea here is that each observation in the dataset will have a subset of features that allow one or more assay-level models to be run; each of these models will produce a posterior probability for each exposure state. The overall posterior probability for that exposure state is then calculated as a weighted average of those assay-level posteriors, where the weights are larger for the assay-level models that had the best cross-validated error rates. The overall posteriors for the different exposure states can then be used in downstream algorithms for biomarker discovery or predictive modeling.

In the results shown below, the assay level models used the ensemble approach (Figure 1D) because assays in general were not measured on the same samples in the cohort. Biomarkers for the top three assay-level models were pooled and used for the ensemble model for the exposure, with weights proportional to the root mean squared error (RMSE) of the cross-validated residuals for those models. The heatmaps depicting model performance show the ensemble model under the label “ESDA” and the assay-level models below it.

2.19. Mapping selected bins to biological function

Selected bins were mapped to biological function based on genome annotation using the Homer bioinformatic tool [15]. Homer was implemented using the annotatePeaks.pl script to generate a tab-delimited output file containing all the selected bins from the ESDA output plus additional gene information for each line in the ESDA document including Gene Ontology enrichment category data. The Homer output is processed through an R script to gather and extract relevant information to annotate the original features with biological process ontologies, export the results in readable format and generate a visualization of the results.

3. Results

Our analysis, funded under the DARPA Epigenetic CHaracterization and Observation (ECHO) program, used pre-existing cohorts from past exposure studies to test the ESDA's ability to find epigenetic signatures associated with those exposures. A summary of the findings appears in Table 1, and includes a list of the exposures we investigated, the sample sizes of the cohorts, the final assays used in the ESDA model and several performance metrics for the resulting models. Details of each of our studies under this program are captured in the sub-sections below describing the features selected from the top three assays by the ESDA model and the ESDA performance metrics. While the utility of the current ESs identified in this work are limited due to the small sample sizes of the cohorts, these signatures can be validated in a future, more thorough study of additional cohorts using more targeted sequencing approaches.

Table 1. The number of samples per group.

Exposure model (group 1 vs. group 2)	Samples in group 1	Samples in group 2	Assays in ESDA model	PPV	NPV	AUC	
MRSA vs. control	13	12	RNA, miRNA, H3K4me3	1.00	1.00	1.00	
MSSA vs. control	7	12	RNA, miRNA, H3K4me3	0.88	1.00	1.00	
MRSA vs MSSA	13	7	RNA, miRNA, H3K4me3	1.00	0.88	1.00	
HIV acute vs. control	10	10	RNA, miRNA, H3K27ac	0.60	0.60	0.76	
HIV chronic vs. control	10	10	RNA, miRNA, H3K27ac	0.75	0.88	0.91	
BA vaccine vs. control	19	35	RNA	0.46	0.70	0.67	
Influenza vs. control	13	13	RNA, miRNA, H3K4me1	0.78	0.65	0.81	
COVID positive vs. prediagnosis	72	101	RNA, EPIC	0.81	0.86	0.89	
COVID positive vs. convalescent	72	104	RNA, EPIC	0.57	0.76	0.73	
COVID convalescent vs. prediagnosis	93	60	RNA, EPIC	0.84	0.83	0.85	
Assays utilized, the positive and negative predictive value and area under the curve for each ESDA exposure model.

Bold rows had an AUC equal to or greater than 0.85 which is a threshold for a diagnostic.

3.1. Negative case study

When there are far more predictors than samples, one concern with biomarker discovery algorithms is that they may be finding noisy markers that happen to separate the samples in the dataset on which they are trained, but which do not generalize to new cohorts within the population. To assess the tendency of the ESDA to “hallucinate” signatures for exposures that do not exist, we conducted a negative case study where we told the algorithm one group of samples was exposed when in fact it was not. The expected behavior would be for the model to fail to find a signature that allowed one to predict the exposure labels with better than coin-flip accuracy. Specifically, 35 control samples were utilized from the B. anthracis vaccination study. RNAseq data from 17 of the control samples were randomly assigned labels of “exposed” and the other 18 controls were assigned labels of ‘normal’. RNAseq assay-level models were built and evaluated using fivefold cross-validation, where model building was repeated five-times with different samples left out and cross-validation accuracies were calculated. The negative case study model produced cross-validation accuracies of 0.40, 0.43, 0.57, 0.51 and 0.54 across the five folds, demonstrating that without a true exposure present in the data, the model cannot accurately differentiate unexposed from exposed and can only achieve coin-flip accuracy. Figure 2 displays the top 10 bins from the negative case study model. As expected, there is no differentiation between the two groups of samples, since they both contain control samples from the same study.

Figure 2. The top 10 bins identified using the pipeline on this sample show no differentiation between exposed (blue) and normal (orange), as expected, since the ‘exposed’ labels were arbitrarily assigned and all of the samples were controls from the same study.

Chr: Chromosome; n: Negative strand; p: Positive strand. Number indicates the first base coordinates in the 500 bp bin.

3.2. Positive case study

While the negative case study shows that the ESDA will not identify a signature when no exposure is present, there is still a question of whether it will find a signature when there are truly exposed samples in the dataset, and whether that signature generalizes to other cohorts within the same population. In the positive case study we conducted, an RNAseq assay level model was built and evaluated using leave one out cross-validation (LOOCV) on 10 acute infected and 10 control samples for HIV data from the iPrEx study (patients in North and South America). The top 10 features from the iPrEx model (Figure 3A) were then applied to an HIV cohort (11 infected and 6 control samples) from the SCOPE study (patients from California, Figure 3B). These features showed similar differentiation between HIV exposed and controls. This is evidence that the ESDA is finding epigenetic biomarkers that not only differentiate HIV samples from controls in the training dataset (the iPrEx samples), but also for an independent cohort of samples from a different HIV study (the SCOPE samples).

Figure 3. Our analytical pipelines identified features that differentiate exposures and are maintained even in cohorts from different continents.

3.3. Exposure signatures

ESDA ensemble models of exposure were built for HIV (acute and chronic), Staphylococcus aureus (MRSA and MSSA), SARS-CoV-2 (naive, acute and convalescent), influenza A and anthrax vaccinations for a total of ten models. Based on our analyses, six of the models identify epigenetic signatures that can differentiate exposure to a chemical or pathogen with an area under the ROC curve (AUC) over 0.85 (Table 1), demonstrating that some signatures are stronger than others and further data may be necessary to identify more subtle differentiating signatures.

3.4. HIV

The iPrEx study cohort (ten paired samples from donors who were negative, at the time of acute diagnosis of HIV, and chronic HIV (on treatment)) was extensively studied through multiple analyses including ChIPmentation on six histone modifications, RNAseq and miRNAseq. Assay level models and two ESDA models were developed comparing acute and chronic HIV to negative, matched samples, respectively. Figure 4A & B depict cross-validated model predictions as heatmaps, where red to yellow shades indicate high probabilities of HIV exposure and blue shades indicate low probability of exposure for each sample in the control or exposed groups. While most individual assay level models showed high variability with false positive or negative calls, three models with the least variation and most consistently accurate calls were used for development of the ESDA models, which generally called more samples correctly than individual assay models. These three models were used in the ensemble model, as described in the methods section above. In the heatmap, the ‘ESDA’ row at the top shows the cross-validated posteriors from the ensemble model. The next section of three rows shows posteriors for the three constituent assay-level models that went into the ensemble. The final rows show additional assay-level models that did not perform well enough to be included in the ensemble.

Figure 4. The predictive posterior probabilities for each assay and ESDA model for (A) acute HIV and (B) chronic HIV. Each box represents the posterior probability for a sample. Rows are assay level models or the ESDA exposure level model. Each column represents a sample. The samples are grouped with the control group under a gray bar on the left side of the figure and the exposure group under a black bar on the right side of figure. The top row is the exposure level ESDA which was developed from the best assay level models (next three rows). Not every assay level model was able to distinguish between groups. Assay level models that were not utilized in the ESDA are shown below the best assay level models. The ESDA model has more blue boxes in the control samples and more red in the exposure samples, showing high posterior probability that the exposure samples belong to the exposure group.

Both ESDA models (HIV Acute and HIV Chronic) utilized the same three data types in the ensemble: RNAseq, miRNAseq and H3K27ac ChIPmentation. The model for acute HIV, which utilized 2518 features, was not strong enough for an accurate diagnosis of HIV with an AUC of 0.76 (see Table 1). However, the chronic model, utilizing 139 features, produced an AUC of 0.91, demonstrating a clear differentiation which could potentially be utilized as part of a diagnostic algorithm where differentiating HIV (on treatment) from other conditions is desired (Table 1). This can be seen visually in Figure 4 as well, where the ‘ESDA’ row for the chronic HIV model shows a clearer pattern of blue (low posteriors) for the control samples and red (high posteriors) for the exposed samples compared with the ESDA model in the acute HIV case.

Figure 5. The predictive posterior probabilities for each assay and ESDA model for SARS-CoV-2 infected and convalescent individuals compared to pre-diagnosis and each other. (A) confirmed positive SARS-CoV-2 (black bar) v pre-diagnosis (gray bar), (B) convalescent infection (gray bar) vs confirmed positive SARS-CoV-2 (black bar), (C) convalescent infection (black bar) v. pre-diagnosis (gray bar). Each box represents the posterior probability for a sample. Rows are assay level models or the ESDA exposure level model. Each column represents a sample. The top row is the exposure level ESDA which was developed from the two assay level models (RNA-seq and EPIC). The ESDA model more separation between the two groups than each assay level model alone.

3.5. SARS-CoV-2

A subset of the CHARM study cohort (101 prediagnosis, 72 acute infection and 104 convalescent participants) was analyzed by RNAseq and EPIC arrays. Prediagnosis samples were defined as being collected at least 14 days prior to the first positive SARS-CoV-2 test. Similarly, convalescent infections were defined as being collected from individuals who have tested negative for SARS-CoV-2 for at least 14 days after having had an infection. Assay level models and three ESDA ensemble models were developed comparing acute and naive (prediagnosis), acute and convalescent infection, and naive and convalescent infection samples (Figure 5A–C). Each ESDA model utilized both EPIC and RNAseq data, which were the only available assays for this dataset. The models that compared patients with acute (Figure 5A & Table 1) or convalescent (Figure 5C & Table 1) SARS-CoV-2 infections readily distinguished these samples from SARS-CoV-2 prediagnosis samples with AUC of 0.89 and 0.85, respectively, with the majority of the exposed group with probability of SARS-CoV-2 (orange-red in the ‘ESDA’ row of the heatmaps) and the majority of the pre-diagnosis group with a probability indicating no exposure (blue in the ‘ESDA’ row of the heatmaps). However, the model comparing acute and convalescent SARS-CoV-2 infection samples was not strong enough for an accurate differentiation with an AUC of 0.73 and using 2366 features (Figure 5B & Table 1) demonstrating the effect of SARS-CoV-2 on the epigenome and transcriptome of people who have been infected are similar to people who are COVID positive confounding differentiation of these two groups. The ESDA model comparing naive, prediagnosis participants and acute, confirmed positive SARS-CoV-2 infected participants selected 34 features: four RNAseq features and thirty EPIC features (Supplementary Table S1). Several of the features were involved in apoptosis, the immune response to viruses, or regulation of gene expression. Interestingly, a methylation feature was found associated with OR5B3, an olfactory receptor. As one symptom of SARS-CoV-2 is a loss of smell, the epigenetic regulation of OR5B3 may be involved and more study will be necessary to elucidate the role of this olfactory receptor in SARS-CoV-2.

3.6. Influenza A H3N2

The BARDA-Vaccitech FLU010 Study cohort (13 matched individuals pre-exposure and post exposure to influenza A H3N2) were analyzed by MINT-ChIPseq on six histone modifications, RNAseq, and miRNAseq. Assay level models and the ESDA model were developed comparing pre- and post-exposure to the influenza challenge in the control group for the vaccine study. The MINT-ChIP data on histone H3K4me1 was best at differentiating exposure from unexposed with an AUC of 0.83, positive predictive value of 0.72 (PPV) and negative predictive value (NPV) of 0.81. The ensemble ESDA model utilized 123 features from histone H3K4me1, RNAseq and miRNAseq models and although the NPV increased to 0.85, the PPV and AUC dropped to 0.53 and 0.81, respectively. The slight drop in predictive power for the ESDA model is due to the lower predictive power of the second and third most accurate assay level models combined with the H3K4me1 (Table 1). This demonstrates that combining multiple low predictive power models in an ensemble does not always improve exposure prediction, and an AUC threshold for inclusion in the ESDA will improve final exposure models.

3.7. Staphylococcus aureus

Seven Methicillin sensitive S. aureus (MSSA), 13 Methicillin resistant S. aureus (MRSA) and 12 control samples) were analyzed by ChIPmentation on six histone modifications, RNAseq and miRNAseq. Assay level models and three ESDA models were developed comparing control to MRSA or MSSA and comparing MSSA to MRSA samples (Figure 6A–C). As shown in Figure 6, several assay level models were not able to accurately predict infection status, with samples from both the control and exposed groups having predictions of the opposite group. However, although this cohort had a low power, the signals were very strong enabling the ESDA model to accurately distinguish control from MRSA or MSSA, and MSSA from MRSA with an AUC of 1.0 (Table 1). Each ESDA model utilized the same three data types, RNAseq, miRNAseq and H3K4me3 ChIPmentation data and neither group had a false prediction (Figure 6A–C). The MRSA vs control model utilized 27 features (5 H3K4me3, 15 miRNA and 7 RNAseq, Supplementary Table S1) while the MSSA model utilized 17 features (9 H3K4me3, 6 miRNA and 2 RNAseq, Supplementary Table S1). While some of the MRSA and MSSA features are on the same chromosome, none were identical loci, and when comparing MRSA to MSSA directly, the model identified 6 features: 2 H3K4me3, 2 miRNAseq and 2 RNAseq (Supplementary Table S1) that could differentiate methicillin resistant from sensitive S. aureus infections based on host response. Interestingly, in the MRSA ESDA model, one of the RNAseq features selected functions in H3-K4 histone methylation and this histone modification was the modification selected for incorporation in the ESDA model as one of the top three significant assays. Furthermore, the functions of several of the features were identified to be involved in immune system regulation, cytoskeleton dynamics, protein turnover and regulation which aligns with a response to the S. aureus infection.

Figure 6. The predictive posterior probabilities for each assay and ESDA model for MSSA and MRSA infection compared to controls and each other. (A) MRSA (black bar) v control (gray bar), (B) MSSA (black bar) vs control (gray bar) and (C) MRSA (black bar) vs MSSA (gray bar. Each box represents the posterior probability for a sample. Rows are assay level models or the ESDA exposure level model. Each column represents a sample. The samples are grouped with the control group under a gray bar on the left side of the figure and the exposure group under a black bar on the right side of figure. The top row is the exposure level ESDA which was developed from the best assay level models (next three rows). Not every assay level model was able to distinguish between groups. Assay level models that were not utilized in the ESDA are shown below the best assay level models. The ESDA model has more blue boxes in the control samples and more red/orange in the exposure samples, showing high posterior probability that the exposure samples belong to the exposure group.

3.8. Anthrax vaccinated

Epigenetic changes could indicate past exposures to not only infectious diseases, but also be utilized to determine vaccination status. Therefore, individuals who were vaccinated against anthrax and worked with Bacillus anthracis in a controlled facility while wearing personal protective equipment were sampled and compared with nonvaccinated individuals who worked in the same facility and did not handle pathogens. RNAseq data were generated from 19 vaccinated and 35 nonvaccinated individuals, and a small subset of these (5 vaccinated and 8 nonvaccinated) individuals were also assessed by ATAC-seq for effects on chromatin remodeling. While some differentiation was observed with RNAseq, the ATACseq data was insufficient to develop a predictive vaccination status model. The ESDA model utilized 2 features and had an AUC of 0.69 (Table 1), demonstrating that vaccination status could not be accurately predicted from RNAseq data alone. Since RNAseq measures transcriptional responses which tend to be more prominent in acute stages of an exposure, such as a vaccination, it was expected that RNAseq data would be the least useful for modeling purposes in this scenario. Further data generation is therefore needed on more long-term epigenetic markers such as methylation and histone modifications to determine if a multiome model can be developed to predict vaccination status.

4. Discussion

Our Exposure Signature Discovery Algorithm (ESDA) provides a solution for identifying key diagnostic biomarkers in the host epigenome that are indicative of exposure to a variety of biological and non-biological materials that alter molecular regulatory mechanisms. By focusing on the commonalities of the host response to these materials, it is possible to determine the cause of a disease or disorder whether it is infectious (virus, bacteria and etiological agent) or non-infectious (chemical exposure or autoimmune) in its origin. Using such a technique has powerful ramifications in terms of developing assays that can rapidly characterize and diagnose the cause of disease in a symptomatic patient, thereby ensuring that the proper treatment can be selected in time. Host-based diagnostics have been in development to enable better discernment between bacterial, viral and noninfectious etiologies of disease [16,17], focusing on the human transcriptional response [16,17]. Here, we are expanding the biomarkers that can be used in a host-based diagnostic to include miRNA, DNA methylation, chromatin remodeling and histone modifications, enabling more precise diagnostics. By identifying a combination of features that are unique to a specific etiological agent and common to a category of disease such as bacterial or viral infection, the diagnostic paradigm can be shifted from one diagnostic assay per target to multiple target or eventually a universal assay to reduce the time it takes to diagnose rare disorders. While the epigenome is complex, with effects from unknown or difficult-to-measure sources such as natural aging, lifestyle decisions, living environment, which could confound the ability to utilize the epigenome for diagnostics, the effect of these unknowns can be mitigated. This is done through the design of the study itself using traditional approaches to select participants that control demographic variables, geographic location and presence of known impacts to the epigenome (e.g., smoking). With large enough cohorts, many of these signals will be dampened through averaging and features associated with these factors will not be selected by the ESDA. Biomarker discovery across multi-omic datasets often generates scenarios where the number of potential markers far outweighs the number of samples with the cohort, and the discussed studies in this work reflect that reality. In these situations, the possibility of selecting false positives is a concern. Exacerbating this issue is the fact that for some exposure types, available data are sparse either because of the low occurrence of exposure or the difficulty in collecting samples from those who were exposed. To pressure test the ESDA, we assembled datasets for a variety of exposure routes, generated multi-omic assay data, standardized those data, used the ESDA to discover biomarkers and then measured the discriminatory power of those markers through a cross-validated predictive modeling process. The overview of those results (Table 1) showed good performance in detecting many of the exposure types, and in cases where the markers failed to predict exposure with an AUC higher than our threshold of 0.85, low sample sizes appear to have played a role.

To further test the ESDA, we conducted validation studies. In one case, we tested the ESDA's tendency to produce false positives by creating an artificial “exposure” group within the control cohort. The ESDA failed to find any biomarkers that could predict the groupings beyond a coin-flip accuracy, which is the expected result in this scenario. In the other validation study, we tested the ESDA's generalizability by training a predictive model on a cohort from the iPrEx HIV study, then predicting on a cohort from SCOPE, an independent study of HIV exposure. Results from this comparison showed that the markers from the iPrEx cohort were able to differentiate positive from negative HIV cases in the SCOPE cohort, demonstrating that the ESDA's biomarkers generalize across cohorts with similar infections or exposures.

Finally, to ensure that the ESDA is finding biologically relevant biomarkers, we used an open-source workflow to map the epigenetic signatures identified to the annotated human genome to determine the function. For example, analysis of the ESDA signature for SARS-CoV-2 identified functions related to response to viral infections, immune regulation and cell death. This demonstrates that the algorithm detects not only discriminatory markers, but also regulatory elements involved in the biological processes underlying host response to an exposure.

Large epigenomic assays, such as EPIC arrays, methylation sequencing, histone modification assays, RNA-seq and miRNA-seq can be thought of as biomarker screening tools. These assays produce millions of data points; however, one needs to parse the data and identify those with diagnostic potential. The ESDA is a way to sift through the millions of potential diagnostic targets to identify the best subset. Once these are known, e.g., in the form of the ESs presented in this work, one can then develop targeted assays specific to these signatures that could be used for diagnostics. Therefore, a future direction is to develop and test rapid assays such as a qPCR based on the ESs described here, thus developing infectious disease diagnostics based on targeting host biomarkers, not the pathogen. Furthermore, while the current work developed binary models, the ESDA is not limited to a binary prediction. It can also accommodate multiple exposure models to find features that discriminate, e.g., between Exposure A vs. Exposure B vs. Exposure C vs. Not Exposed. The features identified by such models would then represent epigenetic changes that are unique to subsets of the exposures in question in such a way that, in combination, the model can distinguish between each of the exposures separately, even when the signatures are overlapping in the epigenome. Therefore, future work will develop multi-way exposure ESDA models to increase specificity of signatures prior to diagnostic development.

5. Conclusion

It is imperative to continue the development of exposure signatures to improve medical diagnostics. By developing a database of multi-omic signatures of infections, chemical exposures, cancer, autoimmune disorders and other disorders of medical importance, host-based diagnostic tests can be developed, which could lower the limit of detection. Furthermore, as more host-based exposure signatures are developed, the future of diagnostics may shift toward a universal, host-based diagnostic to enable diagnosis based on exposure signature with symptomology as supporting information to guide the diagnostic decision tree. Such a universal test based on an exposure signature database could be symptom and pathogen agnostic, saving critical time on rare disorders which currently require multiple tests to identify the underlying issue, providing patients and healthcare providers with answers needed for health management. In its current form, the ESDA is a tool that can be utilized to build such an exposure signature database or identify diagnostic markers from which rapid diagnostic tests can be derived. In future work, the ESDA will be applied on datasets with larger sample size to reduce the signal to noise ratio contributed by confounding factors and identify higher accuracy and generalizable exposure signatures. One aspect of the ESDA that can be exploited for diagnostic development is the use of recursive feature elimination. This affords a unique opportunity to consider a continuum of exposure signatures ranging from a large set of markers that can be used for accurate universal exposure prediction to a small set of markers that may allow the development of limited and targeted assays that still retain reasonable levels of accuracy. We are continuing to expand the utility of the ESDA to develop multi-exposure models, and plan to work toward a pipeline to develop the exposure signature database and achieve the vision of a universal diagnostic exposure screening tool.

Supplementary Material

Supplementary Table S1

Acknowledgments

The authors would like to thank Jack Anderson and Felicia Ruffin for their compassionate field clinic work. We are grateful to the study participants and their families for their wisdom and generosity. We would like to thank Po-Hsu Allen Chen for his contributions to development of the ESDA model and the final revision of the manuscript.

Supplemental material

Supplemental data for this article can be accessed at https://doi.org/10.1080/17501911.2024.2375187

Author contributions

RR Spurbeck was the principal investigator for this work, providing intellectual content, project oversight and execution for analytics, sample collection, manuscript writing and data interpretation.

J Schuetter, A Minard-Smith and B Hill provided substantial contributions to the conception, design and execution of the ESDA, intellectual content, data interpretation and manuscript writing.

JL Beare, AK Smith, HJ Shamma, SC Kahaian and KL Fillinger provided substantial contributions to the conception, design and execution of sample collections for B. anthracis vaccination study. They contributed to the manuscript writing and data collection.

T Chandrasekaran, NC and V Murugan provided substantial contributions to the conception, design and execution of ChIPMentation data collection and manuscript writing.

T Burke, VG Fowler, CW Woods, MT McClain provided substantial contributions to the conception, design and execution of Mint-ChIP, miRNAseq, RNAseq and sample collections for MRSA, MSSA and HIV iPrEx study and manuscript writing.

MAS Amper, W-S Cheng, Y Ge, MC George, K Guevara, N Lovette-Okwara, A Mahajan, N Marjanovic, N Mendelev, CM Miller, S Mofsowitz, VD Nair, G Nudelman, I Ramos, S Rirak, F Ruf-Zamojski, SC Sealfon, N Seenarine, A Soares-Shanoski, S Vangeti, M Vasoya, A Vornholt, X Yu and E Zaslavsky provided substantial contributions to the conception, design and execution of EPIC, multi-ome, scRNAseq, ATACseq assays and sample collections for SARS-CoV-2 and influenza.

L Ndhlovu, M Corley, S Bowler and S Deeks provided substantial contributions to the conception, design and execution of sample collections for HIV from the SCOPE study.

AG Letizia was the Principal Investigator for the COVID-19 Health Action Response for Marines (CHARM) Study.

Financial disclosure

This work was funded by the Defense Advanced Research Programs Agency Biological Technologies Office contract number W911NF19C0041 (RRS) and Defense Health Agency grant 9700140 through the Naval Medical Research Center (AGL). The investigators have adhered to the policies for protection of human subjects as prescribed in AR 70-25. The views expressed in this article are those of the authors and should not be construed to reflect the official policy or represent the positions of the Department of the Navy, Department of the Army, Department of Defense, nor the United States Government. AG Letizia is a military service member. This work was prepared as part of his official duties. Title 17 U.S.C. 105 provides that `copyright protection under this title is not available for any work of the United States Government. Title 17 U.S.C. 101 defines a U.S. Government work as work prepared by a military service member or employee of the U.S. Government as part of that person's official duties. The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed.

Competing interests disclosure

Barinthus Biotherapeutics is a company with financial interest in the influenza vaccine. The work here utilizes samples from their clinical trial provided by TGE and FC. The authors have no other competing interests or relevant affiliations with any organization or entity with the subject matter or materials discussed in the manuscript apart from those disclosed.

Writing disclosure

No writing assistance was utilized in the production of this manuscript.

Ethical conduct of research

All studies involving human subjects or data derived from human subjects were completed under supervision of an institutional review board and all subsequent analyses were carried out after informed consent was obtained. The investigators have adhered to the policies for protection of human subjects as prescribed in AR 70-25. The views expressed are those of the authors and should not be construed to represent the positions of the US Army or the Department of Defense.

Data availability statement

Data will be made available to the research community through DTRA after the DARPA ECHO program has been closed in August 2024.
==== Refs
References

Papers of special note have been highlighted as: • of interest; •• of considerable interest

1. Chakraborty R, Burns B. Systemic Inflammatory Response Syndrome. Treasure Island (FL): StatPearls Publishing; 2023.
2. Herrmann IK, Bertazzo S, O'Callaghan DJP, et al. Differentiating sepsis from non-infectious systemic inflammation based on microvesicle-bacteria aggregation. Nanoscale. 2015;7 (32 ):13511–13520. doi:10.1039/C5NR01851J 26201870
3. Langelier C, Kalantar KL, Moazed F, et al. Integrating host response and unbiased microbe detection for lower respiratory tract infection diagnosis in critically ill adults. Proc Natl Acad Sci. 2018;115 (52 ):E12353–E12362. doi:10.1073/pnas.1809700115 30482864
• Demonstrates the diagnostic power of the host response to enable healthcare providers to better determine cause of respiratory infections and treat patients accordingly.

4. Grant RM, Lama JR, Anderson PL, et al. Preexposure chemoprophylaxis for HIV prevention in men who have sex with men. N Engl J Med. 2010;363 (27 ):2587–2599. doi:10.1056/NEJMoa1011205 21091279
5. Evans TG, Bussey L, Eagling-Vose E, et al. Efficacy and safety of a universal influenza A vaccine (MVA-NP+M1) in adults when given after seasonal quadrivalent influenza vaccine immunisation (FLU009): a phase 2b, randomised, double-blind trial. Lancet Infect Dis. 2022;22 (6 ):857–866. doi:10.1016/S1473-3099(21)00702-7 35305317
6. Chen X, Wang Y, Cappuccio A, et al. Mapping disease regulatory circuits at cell-type resolution from single-cell multiomics data. Nat Computat Sci. 2023;3 (7 ):644–657. doi:10.1038/s43588-023-00476-5
•• Utilizes the same SARS-CoV-2 and MRSA sample data to identify disease regulatory circuits from single cell multi-omics data.

7. Letizia AG, Ge Y, Vangeti S, et al. SARS-CoV-2 seropositivity and subsequent infection risk in healthy young adults: a prospective cohort study. Lancet Respir Med. 2021;9 (7 ):712–720. doi:10.1016/S2213-2600(21)00158-2 33865504
8. Schmidl C, Rendeiro AF, Sheffield NC, et al. ChIPmentation: fast, robust, low-input ChIP-seq for histones and transcription factors. Nat Methods. 2015;12 (10 ):963–965. doi:10.1038/nmeth.3542 26280331
• Method utilized for determination of effect of exposures on histone modifications. This method was used on the S. aureus and the HIV data analysis presented in the current manuscript.

9. van Galen P, Viny AD, Ram O, et al. A multiplexed system for quantitative comparisons of chromatin landscapes. Mol Cell. 2016;61 (1 ):170–180. doi:10.1016/j.molcel.2015.11.003 26687680
• Method utilized for determination of effect of exposures on histone modifications. This method was used on the influenza and the organophosphate/agricultural worker data analysis presented in the current manuscript.

10. Walter L, Van Galen P, Bernstein B, et al. Mint-ChIP3: a low-input ChIP-seq protocol using multiplexed chromatin and T7 amplification. Protocolsio. 2019. https://www.protocols.io/view/mint-chip3-a-low-input-chip-seq-protocol-using-mul-rm7vznk84vx1/v1
11. Neph S, Kuehn MS, Reynolds AP, et al. BEDOPS: high-performance genomic feature operations. Bioinformatics. 2012;28 (14 ):1919–1920. doi:10.1093/bioinformatics/bts277 22576172
12. Guyon I, Weston J, Barnhill S, et al. Gene selection for cancer classification using support vector machines. Mach Learn. 2002;46 (1 ):389–422. doi:10.1023/A:1012487302797
13. Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc B. 1996;58 (1 ):267–288. doi:10.1111/j.2517-6161.1996.tb02080.x
14. Simon N, Friedman J, Hastie T, et al. A sparse-group lasso. J Computat Graph Statist. 2013;22 (2 ):231–245. doi:10.1080/10618600.2012.681250
15. Heinz S, Benner C, Spann N, et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38 (4 ):576–589. doi:10.1016/j.molcel.2010.05.004 20513432
16. Ross MH, Zick BL, Tsalik EL. Host-based diagnostics for acute respiratory infections. Clin Ther. 2019;41 (10 ):1923–1938. doi:10.1016/j.clinthera.2019.06.007 31353133
• Host-based diagnostics can enable precise differentiation of respiratory illnesses of viral, bacterial and noninfectious etiologies, leading to more specific treatment regimens and better antibiotic stewardship.

17. Atallah J, Mansour MK. Implications of using host response-based molecular diagnostics on the management of bacterial and viral infections: a review. Front Med (Lausanne). 2022;9 :805107. doi:10.3389/fmed.2022.805107 35186993
• Addresses the effect on healthcare for infectious disease when implementing a host-based molecular diagnostic assay. By utilizing a host based assay, the diagnostic paradigm shifts from a confirmation of a differential diagnosis based on symptoms to identifying the causative agent based on the host response, reducing misdiagnosis and time between presentation to a physician and diagnosis.
