Compartment-Aware Benchmarking of Respiratory-Virus Transcriptomes for Nasal and Blood Host-Response Modules: A Computational Study

Summary

This computational study presents a compartment-aware workflow for benchmarking public respiratory-virus transcriptomes, demonstrating that nasal and blood host responses show limited gene-level concordance yet yield reproducible, biologically interpretable compartment-specific modules across independent datasets, clinical comparisons, longitudinal recovery, and robustness analyses.

Abstract

Public respiratory-virus transcriptomes are valuable for studying host responses, but differences in tissue source, control definition, and study design can confound pooled analyses. We developed a compartment-aware computational workflow to determine whether reproducible host-response activity can be identified while preserving nasal and blood biological context. The paired GSE117827 paediatric cohort served as the anchor dataset, comprising nasal-swab and whole-blood transcriptomes from symptomatic picornavirus infection, symptomatic respiratory syncytial virus infection, asymptomatic picornavirus detection, and virus-negative controls. After HUGO Gene Nomenclature Committee (HGNC) filtering, 27,685 genes were analyzed. Separate 50-gene protein-coding nasal and blood modules were defined from the top positive responses and locked before external evaluation. Gene-level nasal and blood effects were nearly independent (Pearson r = 0.015), and the modules shared six exploratory rank-overlap genes (Jaccard index = 0.064). Nevertheless, the nasal module separated infection from controls in independent upper-airway cohorts, with areas under the receiver operating characteristic curve (AUROCs) of 0.749, 0.693, and 0.609, whereas the blood module achieved AUROCs of 0.832, 0.924, and 0.870 in external blood cohorts. In longitudinal natural-infection data, matched scores decreased from acute illness to discharge, with paired deltas of 0.436 for nasal samples and 0.330 for blood. Three additional benchmark datasets comprising 666 external samples, together with random-gene nulls, module-size sweeps, bootstrap stability, marker-program correlations, and variance partitioning, defined the robustness and limitations of the workflow. The Pandya 33-messenger ribonucleic acid (mRNA) set remained stronger for viral-versus-bacterial discrimination, while the blood module also increased in bacterial pneumonia. These findings support compartment-specific modules as reusable host-response activity scores for cohort comparison and recovery tracking, rather than universal pan-tissue biomarkers or stand-alone pathogen classifiers.

Introduction

Host transcriptomic signatures are widely used to classify infectious syndromes and compare immune responses across pathogens^1,2,3,4,5. Respiratory viruses often activate overlapping interferon-stimulated genes (ISGs), inflammatory mediators, and antigen-presentation pathways^6,7,8,9,10. Recent paired and longitudinal studies also show that local airway and systemic responses can differ in timing, cellular composition, and magnitude¹¹. This issue extends beyond influenza, respiratory syncytial virus (RSV), rhinovirus, and SARS-CoV-2. Human metapneumovirus remains an important cause of respiratory illness, yet its host-response and intervention literature is still developing^12,13. These observations across respiratory-virus infections further highlight differences between local airway and systemic host responses¹⁴. For computational host-virus studies, the practical question is therefore not only whether a response exists. It is whether public data can be reused through a transparent workflow that preserves tissue context and yields reproducible host-response scores.

Most public respiratory-virus transcriptome datasets were not designed for clean cross-virus or cross-tissue inference. Tissue source, timing, severity, age, control definition, and platform often vary together. High-throughput expression datasets are also vulnerable to unwanted technical and study-level variation^15,16. A pooled analysis can therefore recover a strong interferon or inflammatory signature while obscuring its origin. Nasal samples capture epithelial and mucosal immune biology, whereas whole blood reflects systemic leukocyte responses. Treating these compartments as exchangeable makes a score easy to calculate but difficult to interpret.

We built the analysis around GSE117827, a same-study paired nasal-blood cohort, rather than the largest available discovery collection¹⁷. The cohort is small, but it reduces between-study confounding when nasal and blood effect sizes are compared. It includes symptomatic picornavirus infection, symptomatic RSV infection, asymptomatic picornavirus detection, and virus-negative controls. We therefore used it only as an anchor for separate compartment-specific modules. Module membership was locked before external testing. The stability of the small-cohort selection was examined by stratified bootstrap resampling, random-gene nulls, and module-size sensitivity analyses.

We added a benchmarking layer with bacterial and non-infectious comparators. GSE63990 contains whole-blood acute respiratory illness samples labelled viral, bacterial, or non-infectious¹⁸. GSE40012 contains severe community-acquired pneumonia, systemic inflammatory response syndrome (SIRS), and healthy controls¹⁹. These cohorts test whether a module reflects general host-response activity or pathogen-class specificity. GSE53543 was analysed separately as an ex vivo peripheral blood mononuclear cell (PBMC) rhinovirus perturbation benchmark. It measures leukocyte responsiveness under a controlled stimulus but is not equivalent to natural infection.

Our premise is deliberately conservative. We do not try to prove a universal antiviral signature from heterogeneous public data. Instead, we ask whether a compartment-aware workflow can produce module scores that remain useful after external testing. The intended use is complementary to diagnostic classifiers such as Pandya 33-mRNA. These modules score nasal and blood host-response activity across cohorts, recovery, and symptom gradients. They are not designed to assign pathogen class.

Protocol

This study reanalyzed de-identified, publicly available data and involved no new recruitment, intervention, or specimen collection. Institutional ethics approval and new informed consent were therefore not required for the present secondary analysis. Ethics approval and informed consent for the original studies were reported by the respective data generators.

Study design and workflow logic
We conducted a retrospective public-data bioinformatics study using datasets deposited in the Gene Expression Omnibus^20,21. No new patient samples, cell-culture experiments, animal models, or wet-lab validation were generated. The workflow used an anchor-first design. GSE117827 was used for effect-size estimation, module construction, cross-compartment concordance, and symptom-gradient analysis. Independent datasets were introduced only after module membership had been fixed. The workflow comprised paired-compartment anchoring, HGNC mapping and protein-coding filtering, separate nasal and blood module construction, tissue-matched and longitudinal validation, clinical benchmarking, and robustness analyses. Confirmatory analyses comprised the locked tissue-matched and longitudinal tests. The overlap-gene, symptom-gradient, marker-program, and variance analyses were exploratory or descriptive.

Anchor cohort and primary contrast
GSE117827 contains nasal-swab and whole-blood expression profiles from 26 children: 9 symptomatic picornavirus cases, 5 asymptomatic picornavirus detections, 6 symptomatic RSV cases, and 6 asymptomatic virus-negative controls¹⁷. Two RSV cases lacked blood specimens. The complete anchor design therefore contained 50 samples, and the primary infected-versus-control contrast contained 40 samples after excluding asymptomatic picornavirus detections from module construction. For the primary contrast, symptomatic picornavirus and symptomatic RSV samples were labelled infected, and asymptomatic virus-negative samples were labelled controls. Asymptomatic picornavirus detections were not treated as controls because viral detection without symptoms is biologically distinct from virus-negative health. These samples were held out from module construction and later used for a symptom-gradient analysis.

Probe mapping and preprocessing
GSE117827 was profiled on GPL23126. Transcript-cluster identifiers were mapped to gene symbols using the Clariom D Human na36, hg38 annotation file distributed through GEO platform GPL24539. Only HGNC-approved symbols were retained²². Control probes, ERCC probes, unmapped transcript clusters, and non-approved symbols were removed. Multiple transcript clusters mapping to the same approved symbol were collapsed by median expression. After filtering and collapsing, 27,685 genes were available for the anchor analysis. A log₂(x + 1) transformation was applied only when the 95th percentile of the expression matrix exceeded 50, indicating an unlogged intensity scale. Genes with more than 20% missing values were removed; remaining missing values were replaced by the within-gene median. Each gene was then standardized across all samples in that dataset as z = (x − mean)/sample standard deviation (ddof = 1). Zero-variance genes were assigned a standardized value of 0. Module construction was restricted to entries classified as protein-coding genes in the HGNC complete set.

Effect-size analysis and module construction
Nasal and blood compartments were analyzed separately. Symptomatic viral samples were compared with virus-negative controls using Hedges g standardized mean differences and two-sided Welch tests²³. Hedges g was calculated as the infected-minus-control mean difference divided by the pooled standard deviation and multiplied by the small-sample correction 1 − 3/(4df − 1), where df = n_infected + n_control − 2. Welch tests were implemented with unequal variances. Benjamini–Hochberg false discovery rates were calculated across all tested genes within each compartment²⁴. Because no protein-coding gene met the prespecified strict combination of FDR < 0.05 and Hedges g > 0.8 in the small anchor cohort, positive-effect genes were ranked first by ascending two-sided Welch P value, then by descending Hedges g, with gene symbol as the deterministic final tie-breaker. The top 50 genes formed each primary module. Fifty genes were chosen a priori as a moderate module size that retained biological breadth while limiting missing-gene loss across platforms. This module size was not tuned to external AUROC. Sensitivity analyses using 10, 25, 50, 100, and 200 genes yielded the same qualitative anchor-separation result. The aim was module-level reproducibility, not single-gene discovery.

Cross-compartment concordance
Nasal and blood Hedges g estimates were aligned by gene and compared using Pearson and Spearman correlations. The overlap between the top 50 nasal and top 50 blood modules was summarized by overlap count and Jaccard index. Overlap genes were interpreted as exploratory cross-compartment rank-overlap candidates rather than validated conserved biomarkers.

External tissue-matched validation
Module-score validation used independent public datasets with interpretable source groups. Upper-airway validation included GSE41374 RSV nasal-wash samples and the SARS-CoV-2 datasets GSE152075 and GSE156063. Blood validation included GSE171110 and two separately normalized GSE38900 units: GPL10558 (36 samples; 28 RSV and 8 controls) and GPL6884 (138 samples; 107 RSV and 31 controls). For every external dataset, probes were mapped to HGNC symbols, and duplicate symbols were collapsed by median expression. The same transformation, missing-value filtering, median imputation, and per-gene z-standardization rules used for the anchor dataset were applied within each dataset. A module score was calculated as the unweighted mean of standardized expression values for the represented module genes. AUROC, average precision, and Welch-test FDR were reported as portability metrics, not as estimates of clinical diagnostic performance.

Large clinical benchmark and ex vivo perturbation datasets
Two whole-blood datasets were used as clinical benchmarks rather than discovery cohorts. In GSE63990, samples were assigned from the deposited infection-status metadata to viral respiratory infection (n = 117), bacterial respiratory infection (n = 73), or non-infectious illness (n = 90); samples without one of these unambiguous labels were excluded¹⁸. In GSE40012, deposited diagnosis fields were mapped to influenza A pneumonia (n = 39), bacterial pneumonia without influenza (n = 61), SIRS without pneumonia (n = 40), healthy controls (n = 36), or mixed bacterial/influenza pneumonia (n = 14)¹⁹. Mixed bacterial/influenza samples were described but excluded from binary benchmark contrasts.

GSE53543 contains 196 ex vivo PBMC profiles from 98 individuals. Each individual contributed one media-only sample and one sample exposed to rhinovirus 16 for 24 h; subject identifiers confirmed 98 complete pairs. The principal benchmark reported AUROC across the 98 stimulated and 98 unstimulated samples because AUROC is a rank-separation metric. Subject pairing was retained in the metadata and used in a paired sensitivity analysis of score differences. This experiment was analyzed separately from natural clinical infection and was not used to support clinical diagnostic claims.

Raw GEO series matrices were downloaded for GSE63990, GSE40012, and GSE53543. GPL571 was used for GSE63990, GPL6947 for GSE40012, and GPL10558 for GSE53543. Probes were mapped to HGNC-approved symbols, and duplicate mappings were collapsed by median expression. Deposited normalized matrices were used. The general 95th-percentile rule triggered log₂(x + 1) transformation only for matrices on an unlogged scale. GSE53543 was already log₂-transformed, rank-invariant normalized, and adjusted for processing day by the original investigators, so no additional log transformation was applied. After removing genes with more than 20% missing values and median-imputing the remaining missing values, each gene was z-standardized across all samples within its dataset. The benchmark layer contained 470 clinical whole-blood samples and 196 ex vivo PBMC samples.

Reference signature benchmarking
The anchor nasal and blood modules were benchmarked against three unweighted reference sets. The official Pandya 33-mRNA set was transcribed from Supplementary Table 1 of the original classifier report⁵. The Andres-Terre comparator comprised 33 interferon-oriented genes curated from the reported multi-virus signature³. The Hallmark comparator comprised a 33-gene interferon-alpha core subset linked to the MSigDB HALLMARK_INTERFERON_ALPHA_RESPONSE set^25,26. Complete gene lists for all three reference sets are provided in Supplementary Data 1. These reference scores do not recreate the original weighted classifiers. For each dataset, a score was calculated as the unweighted mean of the available per-gene z scores; at least three represented genes were required, and the number of represented genes was reported. Benchmark contrasts were viral versus bacterial, viral versus non-infectious, viral versus bacterial-or-non-infectious, viral versus healthy/control, and bacterial versus healthy/control where the required groups were available. AUROC and average precision were interpreted as benchmarking metrics. Welch-test P values were adjusted using the Benjamini–Hochberg procedure across all valid dataset-contrast-module combinations in the benchmark table.

Independent longitudinal validation
The paired acute-versus-discharge datasets GSE97741 and GSE97742 were used only for longitudinal validation²⁷. Samples labelled RSV single infection (RSVsi) or rhinovirus (hRV) in the deposited metadata formed the primary analysis; RSV co-infections (RSVco) were retained only in the supplementary output. Acute and discharge labels were parsed from sample titles. Samples were matched by dataset, compartment, virus group, and subject identifier, and only subjects with both time points were retained. The primary combined group comprised 38 RSVsi and 30 hRV pairs (68 pairs per compartment).

No gene selection, module refinement, or threshold tuning used GSE97741 or GSE97742. These datasets were reserved for external validation. Probe identifiers were mapped to approved HGNC symbols, and duplicate symbols were collapsed by median expression. The same 95th-percentile log-transformation rule, missing-value filter, median imputation, and within-dataset per-gene z-standardization procedures were applied before the fixed GSE117827 module scores were calculated.

For each tissue and module, acute and discharge scores were paired by subject and virus group. The acute-minus-discharge difference was calculated for each pair. Cohen dz was calculated as the mean paired difference divided by its sample standard deviation. Two-sided paired t tests and two-sided Wilcoxon signed-rank tests were reported, together with AUROC for acute-versus-discharge score separation. Benjamini–Hochberg FDR values were calculated across all 20 valid paired t tests in the full longitudinal output (two datasets, two module sources, and five prespecified virus-group summaries).

Symptom-gradient and post-hoc analysis
After module membership had been fixed, held-out asymptomatic picornavirus detections in GSE117827 were used only for symptom-gradient analysis. The ordinal clinical score was 0 for virus-negative controls, 1 for asymptomatic picornavirus detections, and 2 for symptomatic infection. Spearman correlation assessed the ordinal trend, whereas the Kruskal–Wallis test assessed overall differences among the three groups. Two-sided Mann–Whitney U post-hoc tests compared the three group pairs. Benjamini–Hochberg correction was applied separately within each compartment to its family of three pairwise comparisons.

Functional enrichment
Nasal and blood module genes were analyzed with Enrichr through gseapy using MSigDB Hallmark 2020, Reactome 2022, and GO Biological Process 2023 libraries^{25,26,28,29,30,31}. Enrichr used its standard over-representation framework based on Fisher's exact test. Library-reported Benjamini–Hochberg-adjusted P values < 0.05 were considered significant. The eight most significant terms per module across the three queried libraries were ranked by adjusted P value. Enrichment results were used only for interpretation.

Robustness, marker-program, and variance analyses
Robustness analyses tested whether the results depended on module size or chance selection. Module sizes of 10, 25, 50, 100, and 200 genes were evaluated. For the random-gene null analysis, the eligible universe comprised HGNC protein-coding genes represented in the processed GSE117827 matrix. Five hundred 50-gene sets were sampled without replacement using a NumPy random seed of 20260622. Bootstrap reselection was restricted to the 1,000 positive-effect protein-coding genes ranked highest in the original anchor analysis for each compartment. In each resample, positive protein-coding genes were ranked by ascending two-sided Welch P value and then by descending mean difference. The top 50 genes were selected. Gene stability was calculated as the number of selections divided by 100. The 20 genes with the highest selection frequency in each compartment were displayed.

Marker-program scores were used as descriptive aids rather than as cell-fraction estimates. Six curated programs were evaluated: epithelial (EPCAM, KRT8, KRT18, KRT19, MUC1, KRT5, KRT14, SCGB1A1, FOXJ1), monocyte/macrophage (LYZ, LST1, S100A8, S100A9, FCGR3A, FCGR1A, CD14, MS4A7, CTSS), neutrophil (S100A8, S100A9, MPO, ELANE, CEACAM8, FCGR3B, OLFM4, MMP8), T/NK (CD3D, CD3E, TRAC, NKG7, GNLY, PRF1, GZMB, KLRD1, IL7R), B/plasma (MS4A1, CD79A, CD79B, MZB1, JCHAIN, IGHG1, SDC1), and myeloid interferon (SIGLEC1, IFI27, IFI44L, ISG15, MX1, OAS1, RSAD2, IFIT3). Each program score was calculated as the mean of the available per-gene z scores, with at least three represented genes required. Spearman P values were adjusted using the Benjamini–Hochberg procedure across the complete family of dataset-module-program correlations. One-way eta-squared was calculated as the between-group sum of squares divided by the total sum of squares for each module-factor pair. Rows missing the relevant factor were excluded from that calculation. This analysis was descriptive and did not adjust for mutually correlated factors.

Reproducibility
All analyses were conducted using scripted workflows in Python 3.12.13, with the software packages and versions reported in the Table of Materials. Randomized procedures used a fixed seed of 20260622. Public GEO expression matrices were processed using a consistent project structure that separated raw inputs, processed data, statistical outputs, tables, and figures. Module definitions, sample-curation records, effect-size estimates, validation statistics, enrichment results, robustness analyses, and figure-source data were retained to support independent verification. The analysis code, dependency specifications, and supporting derived data are available as described in the Data Availability statement.

Results

Anchor-step benchmarking separated tissue-specific from pan-tissue gene effects
The curated GSE117827 anchor matrix contained 27,685 HGNC-approved genes. The primary nasal contrast included 15 symptomatic infected and 6 virus-negative control samples; the blood contrast included 13 symptomatic infected and 6 controls (Table 1). Because this small cohort cannot support stable single-gene discovery, it was used as a paired-compartment anchor. Stability was evaluated at the module level...

Discussion

The main lesson is operational: start with the sampled compartment, not with the largest pooled matrix. GSE117827 is small, but its paired design permits a direct nasal-blood comparison within one study. The near-zero gene-level effect correlation shows that a pooled pan-tissue rank list would hide substantial compartmental structure. Separate modules were therefore the more defensible design. Their external and longitudinal performance supports reproducible host-response activity at the module level, not a universal gen...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Andres-Terre multi-virus interferon-stimulated gene set	Andres-Terre et al.	33-gene comparator; Supplementary Data 1	Unweighted interferon-oriented comparator curated from the reported multi-virus signature; the complete gene list is provided in Supplementary Data 1.
Enrichr	Ma'ayan Laboratory	Accessed 18 August 2026; RRID:SCR_001575	Gene-set over-representation resource accessed through gseapy.
Gene Expression Omnibus	National Center for Biotechnology Information	GEO; RRID:SCR_005012	Public transcriptome repository used to access the analysed datasets.
Gene Ontology	Gene Ontology Consortium	GO Biological Process 2023; RRID:SCR_002811	Functional annotation library used through Enrichr.
GPL10558	NCBI GEO	GPL10558	Platform for GSE53543 and the 36-sample GSE38900 validation unit.
GPL23126	NCBI GEO	GPL23126	Expression platform for GSE117827.
GPL24539	NCBI GEO	Clariom_D_Human.na36.hg38. probeset.csv	Clariom D Human na36, hg38 annotation used for GSE117827 transcript-cluster mapping.
GPL571	NCBI GEO	GPL571	Platform annotation applied to GSE63990.
GPL6884	NCBI GEO	GPL6884	Platform for the 138-sample GSE38900 RSV whole-blood validation unit.
GPL6947	NCBI GEO	GPL6947	Platform annotation applied to GSE40012.
GSE117827	NCBI GEO	GSE117827	Primary paired nasal-swab and whole-blood anchor dataset.
GSE152075	NCBI GEO	GSE152075	External SARS-CoV-2 upper-airway validation dataset.
GSE156063	NCBI GEO	GSE156063	External SARS-CoV-2 upper-airway validation dataset.
GSE171110	NCBI GEO	GSE171110	External SARS-CoV-2 whole-blood validation dataset.
GSE38900	NCBI GEO	GSE38900; GPL10558 and GPL6884	External RSV whole-blood validation analysed as separately normalized 36- and 138-sample platform units.
GSE40012	NCBI GEO	GSE40012	Clinical influenza A pneumonia, bacterial pneumonia, SIRS, and healthy-control benchmark.
GSE41374	NCBI GEO	GSE41374	External RSV nasal-wash validation dataset.
GSE53543	NCBI GEO	GSE53543	Paired ex vivo PBMC rhinovirus perturbation benchmark involving 98 subjects and 196 samples.
GSE63990	NCBI GEO	GSE63990	Clinical viral, bacterial, and non-infectious whole-blood benchmark.
GSE97741	NCBI GEO	GSE97741	Longitudinal whole-blood acute-versus-discharge validation dataset.
GSE97742	NCBI GEO	GSE97742	Longitudinal nasopharyngeal acute-versus-discharge validation dataset.
gseapy	gseapy developers	Version 1.3.1	Python interface used to query Enrichr.
HGNC gene symbol resource	HUGO Gene Nomenclature Committee	Accessed 18 August 2026; RRID:SCR_002827	Approved gene-symbol normalization and protein-coding classification resource.
Matplotlib	Matplotlib development team	Version 3.11.1	Figure generation.
MSigDB Hallmark gene sets	Broad Institute	Hallmark 2020; RRID:SCR_016863	Hallmark enrichment library accessed through Enrichr.
MSigDB Hallmark interferon-alpha response gene set	Broad Institute	HALLMARK_INTERFERON_ALPHA_ RESPONSE; 33-gene core subset in Supplementary Data 1	Unweighted interferon-alpha comparator; the exact subset is provided in Supplementary Data 1.
NumPy	NumPy developers	Version 2.5.2	Numerical computing and seeded random sampling.
pandas	pandas development team	Version 3.0.5	Tabular data processing.
Pandya 33-mRNA host-response gene set	Pandya et al.	Supplementary Table 1; Supplementary Data 1	Official 33-gene set scored as an unweighted comparator.
Python	Python Software Foundation	Version 3.12.13	Computational analysis environment.
Reactome	Reactome	Reactome 2022; RRID:SCR_003485	Pathway enrichment library accessed through Enrichr.
scikit-learn	scikit-learn developers	Version 1.9.0	AUROC and average-precision calculations.
SciPy	SciPy developers	Version 1.18.0	Welch, paired t, Wilcoxon, Mann–Whitney, and correlation tests.
Seaborn	Seaborn development team	Version 0.13.2	Statistical graphics.
statsmodels	statsmodels developers	Version 0.14.6	Statistical utilities.