Clinical characterisation and classification of dengue patient cohort

In this study, 112 NS1 antigen–positive hospital-admitted primary dengue infection patients were enrolled to explore the host transcriptome dynamics and alternative splicing to elucidate their potential role in modulating disease severity. Blood samples were collected from these patients, and total RNA was isolated for a comprehensive assessment of the whole transcriptome (Fig 1A). The patients were stratified into three severity sub-phenotypes based on the 2009 WHO classification—dengue without warning signs referred to as mild (n = 45), and dengue with warning signs that was further sub-grouped into moderate (n = 46) and severe (n = 21); based on the presence of two key clinical parameters of leukopenia and thrombocytopenia (5). Using the bulk RNA-sequencing data, we performed both gene-level and transcript-level differential expression analyses between severe versus mild, severe versus moderate, and moderate versus mild patient groups, followed by differential splicing analysis and differential transcript isoform usage analysis to investigate factors modulating dengue severity sub-phenotypes.

Upon assessing the clinical parameters associated with dengue-infected patients, we observed no statistically significant differences in age and gender between the severity sub-groups (Table 1 and Fig 1B). This ensured the absence of age- and gender-related potential biases in subsequent analyses. Severe thrombocytopenia and leukopenia were observed in the severe dengue patients, whereas moderate patients showed intermediate platelet counts but severe leukopenia compared with the mild patients (Fig 1C and D). Notably, we also observed an increase in lymphocyte count in moderate and severe conditions (Fig 1E). Moderate and severe forms of dengue were also characterised by neutropenia (Fig 1F). We also observed high bilirubin and serum glutamic–oxaloacetic transaminase levels in the severe dengue (Fig 1G), which are indicators of liver damage. The evaluated clinical parameters thus re-affirm the severity-based classification of the dengue patients. Other symptoms of dengue, including vomiting, diarrhoea, fever, and headache, were, however, not significantly different across the patients of varying severity (Fig 1H). Although study samples were NS1 antigen–positive for dengue, we also tried to identify dengue serotype(s) for these patients. Towards that, we captured viral reads from the bulk RNA-seq data using the non-human mapped reads. Dengue genome analysis revealed that most of the samples showed DENV2 infection (Table S1). For some of the samples, we had inadequate viral reads for serotype identification, which were distributed across the disease severity.

Expression dynamics of genes and transcripts across severity sub-phenotypes

In response to the viral infection, the host transcriptome is rapidly altered to either promote host protective functions to limit infection-driven pathology, or augment survival and spread of viral progenies, leading to severe infection. We performed a detailed analysis of the host gene- and transcript-level differential expression, identifying events that contribute to the dynamic changes in the host transcriptome across disease severity (mild, moderate, and severe). Taking a significance threshold of P adj < 0.05 and log₂ fold change ≥ |1.5|, for both differential gene expression (DGE) and differential transcript expression (DTE), we identified 214 DTEs between severe and mild patients (Fig 2A and Table S2A). Among these, 10 were down-regulated, and 204 were up-regulated in the severe patients. In the comparison between the severe and moderate cases, 11 DTEs were up-regulated, and 15 were down-regulated in the severe patients (Fig 2B and Table S2B). Interestingly, in the moderate versus mild comparison, all 122 DTEs were down-regulated in the moderate (Fig 2C and Table S2C).

In the case of DGEs, we identified 595 genes between the severe and mild patients (21 down-regulated, 574 up-regulated) (Fig S1A and Table S3A). Between severe and moderate, 37 genes were captured (2 down-regulated, 35 up-regulated) (Fig S1B and Table S3B). In the comparison between moderate and mild, 17 genes were down-regulated, whereas four were up-regulated (Fig S1C and Table S3C). There was a minimal overlap between the DGEs and DTEs across the disease severity sub-phenotypes. Notably, only 47 genes were common in the severe versus mild comparison, with three genes shared between DGEs and DTEs in both the severe versus moderate and moderate versus mild comparison groups (Fig 2D). These findings emphasise unique patterns in the expression dynamics of both genes and transcripts associated with various dengue disease severity in our study cohort.

We delved further into specific transcript features to identify potential enrichment/depletion of severity-specific isoform patterns that were significantly expressed compared with their non-significant counterparts/isoforms, referred to as differential isoforms (DIs) (Fig 2E). In contrast to the global distribution of biotypes, a notable decrease in the conventional protein-coding biotype was evident (P = 0.03), coupled with a simultaneous increase in non-canonical transcripts, such as long non-coding RNA (lncRNA) (P = 0.01) and retained intron transcript biotypes (P < 0.001) among severe dengue patients (Fig 2F). Protein-coding transcripts were, however, significantly enriched in both the mild (P < 0.001) and moderate (P = 0.003) patients. Although the moderate group did not show any presence of significantly expressed lncRNAs, mild patients showed significant deficiency of lncRNA transcripts (P < 0.001). Furthermore, transcripts destined for nonsense-mediated decay (NMD) were marginally but significantly enriched in the moderate but depleted in mild patients (P < 0.001), whereas we observed a significant decrease in such transcript biotypes in the severe patients (P < 0.001). Overall, a significant variation was observed in the proportions of transcript biotypes across the disease severity.

We performed the gene set enrichment analysis, considering pathways with P < 0.05 and with at least two genes as significant; we noted significant positive enrichment of immune response pathways including B-cell receptor signalling, complement activation, Fc receptor–mediated signalling, innate immune response, and infectious disease response, at both the gene and transcript isoform levels (Table S4 and Fig 2G). Interestingly, we also observed differential regulation of critical pathways exclusively at the transcript isoform level. Specifically, there was enrichment of transcript isoforms that regulate chromatin modifications and cell cycle in the severe versus mild comparison. Although there was significant positive enrichment of DTEs regulating programmed cell death (PCD) pathways in the severe versus mild condition, there was negative enrichment of DTEs affecting this pathway in the moderate versus mild condition. We also observed the ablated expression of DTEs affecting interferon signalling in severe versus mild and an antiviral response in severe versus moderate. Given that increased haemolysis is a characteristic of severe dengue, we observed enrichment of DTEs affecting haemolysis in the severe versus mild conditions.

Finally, we explored the genomic features of the DTEs across the disease severity (Fig S1D–F). We found that significant DTEs have lower numbers of exons in the severe versus mild patients, when compared to the non-significant differential isoforms (DIs) (P < 0.001) (Fig 2H). In addition, we noted a decrease in transcript length, particularly in the coding sequence (CDS), in the severe patients, between significant DTEs compared with the DIs (P = 0.008). These findings suggest a reduction in the length of transcript isoforms associated with severe dengue similar to our previous study on SARS-CoV-2 differential disease severity (Fig S1G–I) (27). The decreased number of exons observed in severe dengue aligns with the diminished length of transcript isoforms that exhibit differential expression (DTEs) in these individuals. Importantly, despite the shorter transcripts, severe dengue is characterised by greater transcript diversity, as evidenced by an increase in the ratio of differentially expressed transcript isoforms in significant DTEs when compared to the DIs (P < 0.001) (Fig 2I). These findings suggest that dengue severity is linked to a distinct transcriptional/post-transcriptional regulation that results in the formation of shorter yet diverse transcripts.

Escalated differential splicing associated with dengue disease severity

In our earlier findings, we noted an increase in transcript diversity with the severity of dengue, accompanied by a reduction in the proportion of protein-coding biotypes. To understand the distribution of isoform diversity, we delved into alternative splicing, which is one of the main contributors to the heightened transcript diversity. In our assessment of splicing patterns across severity, we identified 1,040 genes with significant differential splicing between the severe and mild patients (Table S5). In contrast, there were only 141 genes with significant differential splicing events in the severe versus moderate and 235 genes in the moderate versus mild patients (Fig 3A). Despite variations in transcript expression levels being potentially influenced by splicing events, there was a limited overlap observed between DTEs, DGEs, and differentially spliced (DS) genes across the disease severity. Specifically, we noted an overlap of 13 transcript isoforms between the DTEs and DS in severe versus mild, compared with only four and six overlaps in severe versus moderate and moderate versus mild, respectively. Hence, although dengue severe states are characterised by remodelling of the host transcriptome splicing landscape, it seems that these changes do not necessarily impact the levels of gene or transcript expression.

To elucidate the functional implications of differential splicing on gene function, we investigated the pathways associated with the genes undergoing differential splicing. We specifically observed the enrichment of immune response pathways in severe compared with mild conditions (Fig 3B). These pathways included TLR, nod-like receptor, and C-type lectin receptor innate recognition; neutrophil response; interferon, TNF, and interleukin cytokine response; Fc receptor–mediated response; MHC-I and MHC-II; and T-cell and B-cell signalling pathways. Molecular signalling pathways associated with spliced transcripts included NFκB, PI3-Akt, Wnt/beta-catenin, Rho GTPase, and VEGF signalling pathways. Other homeostatic functions affected by splicing changes included the cell cycle, post-translational modifications, alternate splicing, epigenetic regulation, DNA damage response, and programmed cell death pathways. Platelet activation and erythropoietin signalling were also enriched among spliced transcripts in the severe patients. In addition, we observed significant splicing changes in the genes associated with pathogen entry into the host, endocytosis, and phagocytosis. The severe versus moderate exhibited enrichment of spliced transcripts associated with the innate immune response, along with RNA and protein modification pathways (Fig 3C). In contrast, exclusive enrichment of innate immune response pathways, such as neutrophil activation, TLR recognition, and interferon and interleukin signalling, was observed in moderate versus mild patients (Fig 3D).

We also observed severity-specific enrichment or reduction in certain splicing events in comparison with their global distribution (Fig 3E). All splicing events, except for skipped exons (SE), showed significant differences across all three groups in differentially spliced (DS) genes. The alternative 5′ splicing site (A5SS) and retained intron (RI) events exhibited a consistent reduction in the DS genes across all severity conditions in comparison with the global splicing event (P < 0.001 for all three groups). Interestingly, in severe DTE genes, RI was significantly higher compared with the global and DS genes (P = 0.001), which aligns with the higher presence of a retained intron biotype in the severe group compared with the global and moderate/mild in our previous result (Fig 2F). In DS genes, we see a slight albeit significant reduction in the prevalence of alternate 3′ splicing site (A3SS) events (P = 0.016, P < 0.001, and P = 0.002 in severe, moderate, and mild groups, respectively). A3SS was significantly higher in the DTEs of severe and moderate patients (P < 0.001 and 0.017, respectively), but lower in mild (P < 0.001), when compared to the overall global splicing events in our dataset. In addition, we observed significant enrichment of mutually exclusive exon events in the DS transcripts of all severity groups (P < 0.001 across all three groups) and DTEs of moderate and mild compared with the global events (P = 0.012 and <0.001, respectively). The observed changes in the distribution of splicing events indicate potential alterations in transcript isoform diversity, likely influencing the functional repertoire of transcripts associated with severe disease.

Our previous investigation identified preferential enrichment of repeat elements, specifically short interspersed nuclear elements (SINEs), near the splicing sites of differentially expressed transcripts in the COVID-19 patients (27). This observation prompted us to investigate the overall presence of repeat elements, specifically within 200 bases downstream of 5′ splice sites (5′ss) and upstream of 3′ splice sites (3′ss) of DTEs across dengue severity (Fig 4A). Intriguingly, we noted a substantial increase in the overall abundance of repeat elements in significant DTEs of severe patients when compared to their non-significant isoforms (DIs) (Fig 4B). This contrasted sharply with DTEs in mild and moderate patients, where we observed significant depletion of repeat elements at both the 5′ss and 3′ss. We then examined the distribution of specific repeat classes, including simple repeats, LTRs, SINEs, and long interspersed nuclear elements (LINEs), at the 5′ss and 3′ss. In the severe group, we observed abundance of LINEs (P < 0.001) and LTR elements (P < 0.001), accompanied by depletion of SINEs (P < 0.001) and simple repeats (P = 0.038) in significant DTEs at 5′ss compared with DIs (Fig 4C). Conversely, in the moderate group, there was notable enrichment of SINEs (P < 0.001), whereas both moderate and mild groups exhibited significant depletion of LINE (P = 0.014 and 0.009, respectively) and LTR elements (P < 0.001 in both) at the 5′ss of significant DTEs. At 3′ss, SINE depletion (P < 0.001) and LTR enrichment (P < 0.001) were observed in the severe patients (Fig 4D). On the contrary, in the moderate group, there were SINE enrichment (P = 0.004) and LTR reduction (P = 0.0016) at the 3′ss of DTEs. Significant DTEs in the patients with mild symptoms displayed LINE depletion (P < 0.001) and LTR enrichment (P < 0.001), along with simple repeat enrichment (P = 0.035) at 3′ss. Thus, the alteration observed in the repeat distribution patterns across 5′ss and -3′ss from mild to severe patients could potentially be associated with alternative splicing events and isoform diversity in DTEs across severity groups. Understanding how these diverse transcript biotypes collectively contribute to the ultimate gene function is crucial.

Differential isoforms impact on overall gene expression and functional outcomes

Given the substantial functional changes observed in DS transcripts across different levels of severity, we further examined the expression patterns of diverse transcript isoforms to comprehend the functional consequences of genetic heterogeneity arising from transcript diversity on dengue disease severity. To assess this, we examined the DTU across different severity levels. The DTU provides an accurate measurement of how individual transcript isoforms contribute to the overall gene expression (Fig 5A). In total, we identified 121 DTU events among the expressed transcripts in severe versus mild, 92 DTU events in moderate versus mild, and 7 DTU events in severe versus moderate patients (Fig 5B and Table S6). It is notable that there was a minimal overlap between the groups and no overlaps across all three comparisons. Interestingly, we observed varying degrees of an overlap between DGE, DTE, DS, and DTU events across paired comparisons between the three severity groups (Table S7). Specifically, we identified 40 genes that exhibited both DS and DTU events in the severe versus mild comparison, and only nine genes in moderate versus mild and no overlaps in severe versus moderate groups (Fig 5C–E). This underscores the inherent complexities/dynamics associated with gene expression, emphasising that no single regulatory mechanism at a particular stage can be deemed to exert sole dominant influence over the ultimate expression and functional identity of a gene.

Nevertheless, given the substantial diversification of functional pathways associated with DS transcripts in severe dengue (Fig 3B), we chose to explore the functional enrichment or reduction in specific splicing biotypes among genes that exhibited both DS and DTU events (Fig 5F). Remarkably, we observed a severity-dependent, significant decrease in the overall abundance of protein-coding biotypes within various functional pathways related to dengue pathogenesis and the host immune response. This suggests a potential impairment of these human host responses in severe dengue. We observed specific ablation of protein-coding DTU and DS transcripts involved in blood clotting, a process in severe dengue that could potentially influence increased haemorrhagic response in severe dengue. Despite the observed up-regulation of DTE related to the cell cycle and chromatin modifications in severe dengue (Fig 2G), transcript usage analysis revealed that DS transcripts involved in epigenetic regulation and the cell cycle exhibited higher abundance of non-functional biotypes such as protein-coding CDS not defined, and retained intron, respectively, in severe patients (Fig 5F). Our observations also indicate a general decrease in functional biotypes related to cell adhesion, cytoskeletal remodelling, cellular metabolism, and intracellular signalling functions during severe dengue infection.

Furthermore, there was suppression of protein-coding biotypes associated with innate and adaptive immune responses, which could potentially drive attenuated immune responses resulting in severe infection. Surprisingly, we also observed a decrease in the protein-coding biotypes of transcripts associated with immunosuppression in severe dengue patients. In transcripts related to RNA processing, we observed increased abundance of functional protein-coding biotypes and the elimination of non-functional NMD transcripts. This result could potentially account for the dramatic enrichment of DS genes in severe patients attributable to enhanced splicing.

Functional consequence of DS, DTE, and DTU events on gene expression underlying dengue disease severity

To elucidate the functional implications of the DS and DTU with respect to disease severity, we focused on two genes, MX1 and RBM39, which are key components of the host immune response. As both these genes also exhibit DTE between the severe and mild patients but showcasing distinct directions of regulation, MX1 shows down-regulation, whereas RBM39 shows up-regulation in the severe patients compared with the mild. The MX1 gene, recognised as a type I/type III interferon–inducible viral restriction factor, demonstrates GTPase activity known to impede viral replication by disrupting the formation of the viral polymerase–ribonucleoprotein complex (28, 29, 30). Notably, differential intron excision was observed at two specific junctions of the gene (Fig S2A), where both the junctions corresponded to the smallest protein-coding isoform (MX1-216) and the significantly expressed transcript being up-regulated in the mild patients (Fig S2B). The differential splicing generates multiple transcript biotypes in the MX1 gene, predominantly comprising protein-coding variants and one retained intron, collectively resulting in overall suppression of MX1 gene function in severe conditions (Fig S2C).

On the contrary, the RNA binding motif protein 39 (RBM39) gene, identified as a novel regulator of innate immune responses and alternative splicing during viral infection, has been demonstrated to directly influence the transcription and/or splicing of interferon-stimulated genes, notably affecting the expression of interferon regulatory factor 3 (31 Preprint). RBM39 displays distinct exon usage patterns associated with exons in the shortest protein-coding isoform (RBM39-214), and these patterns are also significantly differentially expressed (Fig 6A). RBM39 splicing activity is potentially linked to the efficient recognition of 3′ splice sites of RNA (32). It exhibits high expression in the immune compartment, including natural killer cells, T and B cells, myeloid cells, and other CD34-positive haematopoietic cells. Although an up-regulation of splicing around protein-coding transcripts is observed, differential splicing contributes to NMD biotype diversity (Fig 6B). However, a down-regulation in an NMD biotype is noted across severity, coupled with an up-regulation in protein-coding transcripts (Fig 6C). Differential splicing events in both the genes contribute not only to an increase in transcript biodiversity but also to the regulation of differential gene expression. Solely focusing on gene expression analysis might overlook the nuances associated with these transcript biotypes, which, when considered collectively, provide a more comprehensive understanding of the host response at the gene expression level. By examining results at the transcript level, we can pinpoint specific isoforms that play crucial roles in the major functions of these genes underlying differential dengue disease severity.

Sample collection and clinical classification

The present investigation, carried out at CSIR-Institute of Genomics and Integrative Biology (IGIB) in New Delhi, India, included a cohort of 112 hospital-admitted individuals diagnosed with dengue through NS1 antigen testing. After the acquisition of written consent from the participants, the paramedical team at MAX Super Specialty Healthcare Hospital, New Delhi, India, gathered 2–3 ml blood samples during the initial hospital visit using EDTA-coated vials, adhering to the principles outlined in the Declaration of Helsinki. This sample collection took place between August 2022 and November 2022, which coincides with the peak dengue infections in New Delhi, India. All pertinent clinical information was retrospectively obtained by scrutinising the electronic medical records of each participant.

Based on the complete blood count reports of the patients, the study categorised the 112 dengue NS1-positive patients into three groups: dengue without warning signs (45 patients with normal platelet and leucocyte counts), and dengue with warning signs (comprising two sub-groups—46 patients with normal platelet but decreased leucocyte counts [leukopenia] and 21 patients with both low platelet counts [thrombocytopenia] and decreased leucocyte counts [leukopenia]). None of the patients exhibited severe bleeding, organ failure, or abnormal liver parameters, in accordance with the World Health Organization (WHO) dengue classification and management scheme (5). Consequently, a separate group for severe dengue was not formed. For clarity, the groups are designated as mild (no thrombocytopenia and leukopenia), moderate (no thrombocytopenia but leukopenia), and severe (thrombocytopenia and leukopenia) dengue cases.

Nucleic acid isolation

RNA extraction from the blood samples was carried out using QIAGEN QIAamp RNA Blood Mini Kit, following the manufacturer’s instructions with specific modifications. The incubation time and centrifugation time during the erythrocyte lysis step were reduced to 5 min. In addition, for enhanced RNA purification, a 3-min incubation period was introduced during all washing steps. The purity of the isolated RNA was assessed using NanoDrop Microvolume Spectrophotometer and confirmed by running it through an agarose gel. Subsequently, the RNA was stored for short time at −80°C until thawed for total RNA-seq library preparation.

Library preparation and sequencing

Libraries were prepared with a total of 250 ng of total RNA using Illumina TruSeq Stranded Total RNA Library Prep Globin Kit (Cat. No. 20020612; Illumina). Initial steps involved the depletion of globin mRNA and ribosomal RNA (both cytoplasmic and mitochondrial), as these two forms of abundant RNA are present in high levels in the whole blood. Cleaved RNA fragments underwent first-strand cDNA synthesis through reverse transcriptase and random primers. Subsequently, double-stranded cDNA was synthesised using DNA polymerase 1 and RNase H, followed by purification using AMPure XP (A63881; Beckman Coulter). Adenylation at the 3′ blunt end of double-stranded cDNA was performed, and sequencing libraries were uniquely indexed and enriched through PCR amplification.

Library quality was confirmed through size analysis on Agilent 2100 Bioanalyzer, and library concentrations were determined using the Qubit double-stranded DNA (dsDNA) high-sensitivity (HS) assay kit (catalogue no. Q32854; Thermo Fisher Scientific). Libraries were then diluted to 4 nM and combined in an equimolar fashion, with 24 samples per library pool. Paired-end, 2 × 151 read length sequencing was performed on a NextSeq 2000 platform (Illumina) at a final loading concentration of 650 pM.

Differential gene and transcript expression analyses

The raw sequencing reads were demultiplexed using bcl2fastq (v), and subsequently, read quality was assessed using FastQC (64, 65). To filter out low-quality bases (phred score < 20) and trim any adapter content, Trimmomatic v.0.39 was used (66). The filtered reads underwent a quality check to ensure effective filtration, and quantification was carried out using the pseudo-alignment method in Salmon (v1.8.0) against the human reference genome (GRCh38.110.primary assembly from Ensembl) (Fig S3) (67).

During quantification with Salmon, the --numGibbsSamples parameter was set to 10 inferential replicates to generate bootstrapped abundance estimates for each sample using posterior Gibbs sampling to account for the variability in transcript abundance. The quantified files were imported to the R environment using the tximport R package for differential expression (DE) analysis (68).

For gene-level DE analysis, the DESeq2 package was employed and genes with >10 reads in at least 21 samples (the size of the smallest group) were retained for stringent further analysis (69). To control batch effects in the data, the removeBatchEffect function from the limma R package was used in conjunction with the model design in DESeq2 (70). For transcript-level DE analysis, the swish method from the fishpond (v 1.0.0) R package was used, using the inferential replicates (71). Similar to DESeq2 analysis, transcripts with >10 reads in at least 21 samples were kept for differential expression analysis. Unannotated significant genes were further filtered.

Differential gene and transcript expression analyses were conducted between different severity patients’ groups, namely, severe versus mild, severe versus moderate, and moderate versus mild. The Benjamini–Hochberg (BH) correction was applied to address multiple comparisons. DEGs and DETs with an adjusted P < 0.05 and | log₂ fold change (log₂FC)| ≥ 1.5 were considered statistically significant.

Gene set enrichment analysis

Gene set enrichment analysis using the filtered lists of differentially expressed (DE) genes and transcripts was performed. The human MsigDB collection, specifically the Reactome pathway database (c2.cp.reactome.v2023.1.Hs.symbols.gmt), was employed for this analysis through the fgsea R package (72, 73 Preprint). DE statistics served as the basis for ranking genes, and for transcripts, the average of only the significant isoforms belonging to the same gene was considered for ranking. Pathways with at least two genes mapping to the pathway and a P-value cut-off of < 0.05 were considered significant (Table S4). Redundant or similar pathways were grouped together for clarity. We used ggplot2 in R to visually represent the significant pathways, considering both the combined score and the number of genes involved in these pathways (74).

Differential splicing analysis

The filtered fastq files were aligned against the reference human genome (GRCh38.110) using the splice-aware aligner STAR (v2.7.11) (75). The aligned bam files were used for detecting both event-based splicing detection and differential splicing quantification between the severity groups. For event-based differential splicing detection, rMATS (v4.2.0) (replicate Multivariate Analysis of Transcript Splicing) was employed with the ‐‐variable-read-length and ‐‐allow-clipping parameters (76). The statistical model of rMATS detects differential splicing by employing a hierarchical framework to model exon inclusion levels, considering estimation uncertainty and variability among replicates.

For quantification of differential splicing between severity groups, the LeafCutter v0.2.9 tool was used (77). For this, the bam files were transformed into junction files using regtools. Subsequently, intron clustering and analysis of differential splicing were performed using LeafCutter, with allowed intron sizes set between 50 and 500 kb. The human reference was used for annotation. The groups for differential splicing analysis were organised in a manner consistent with those used for the comparison in the differential gene/transcript-level analysis. The results were visualised using LeafViz, and a report was generated using RStudio.

Viral read quantification

Using the reads unmapped to humans, alignment was performed against the serotype-specific DENV genomes for DENV1, DENV2, DENV3, and DENV4 (NC_001477, NC_001474, NC_001475, NC_002640) using BWA-MEM to generate bam files (78). The viral genome coverage and depth were determined using SAMtools (Table S1) (79).

Repeat element enrichment analysis

We analysed the repeat elements around the splicing sites for the DTEs and their non-significantly expressed isoforms (referred to as differential isoforms [DIs]). The genomic sequences 200 bases upstream of the 3′ splice site and 200 bases downstream of the 5′ splice site were extracted for each transcript from the GRCh38 reference genome using BEDTools (80). To query for repeat elements, all the sequences were uploaded to RepeatMasker using hmmer search engine and Human as DNA source (81). The repeat elements considered were LINEs (L1 and L2), SINEs (Alu and MIR), LTR, and the simple repeats (Table S8). The repeat bases were normalised for total repeat bases and represented as proportion. The significance was compared between DTEs and the DIs using the χ2 test in GraphPad Prism, and the data were visualised in Excel and modified using Inkscape.

Differential transcript usage analysis

Transcript usage was evaluated, and variations in the proportional impact of a transcript on the overall gene expression were computed by determining its relative abundance. The DTU analysis was performed using the salmon alignment-free abundance estimations, employing direct read quantification at the transcript level. The DRIMSeq R package was used for the analysis, relying on a Dirichlet-multinomial model for individual genes and estimating a gene-wise precision parameter (82). This approach directly modelled the correlation among transcripts within their parent gene.

To mitigate false discovery rates stemming from the complexity of transcript diversity and the number of transcripts per gene, transcripts with fewer than 10 read counts or with a relative contribution to the overall gene expression of less than 0.1 were excluded. In addition, genes with fewer than 10 counts in any samples were omitted. The resulting filtered set of transcript-level counts was used as input for the DRIMSeq analysis to investigate alterations in transcript usage between severity groups. The P-values reported at the gene level were further corrected within stageR using the BH FDR procedure (83). For the screened genes, transcript-level P-value is adjusted to control for a family-wise error rate. A transcript was deemed a DTU event if the P-value was less than 0.05. A gene was considered differentially used if it had at least one DTU event.

Statistical analysis and data visualisation

Wherever appropriate, differences between continuous data points were assessed using the two-tailed Mann–Whitney U test, whereas categorical data comparisons were conducted through chi-square testing. All statistical analyses were carried out using a licensed version of GraphPad Prism. For data visualisations, a combination of ggbio (86), GenomicRanges (87), biovizBase (88), ggplot2 (74), ggbreak (89), ggraph (90), dendextend (91), igraph (92), colormap (93), EnhancedVolcano, GViz (94), rawgraphs (95), SRPlots (96), and Microsoft Excel was employed. Figures were subsequently modified using Inkscape (97). A significance threshold of P < 0.05 was applied, unless otherwise specified.