Accurate Quantification of Overlapping Herpesvirus Transcripts from RNA Sequencing Data

Abstract

Herpesviruses employ extensive bidirectional transcription of overlapping genes to overcome length constraints on their gene product repertoire. As a consequence, many lytic transcripts cannot be measured individually by RT-qPCR or conventional RNA-seq analysis. Bruce et al. (Pathogens 2017, 6, 11; doi:10.3390/pathogens6010011) proposed an approximation method using Unique CoDing Sequences (UCDS) to estimate lytic gene abundance from KSHV RNA-seq data. Although UCDS has been widely employed, its accuracy, to our knowledge, has never been rigorously validated for any herpesvirus. In this study, we use CAGE-seq as a gold-standard to determine the accuracy of UCDS for estimating EBV lytic gene expression levels from RNA-seq data. We also introduce the Unique TranScript (UTS) method that, like UCDS, estimates transcript abundance from changes in mean RNA-seq read-depth. UTS is distinguished by its use of empirically determined 5’ and 3’ transcript ends, rather than coding sequence annotations. Compared to conventional read assignment, both UCDS and UTS improved quantitation accuracy of overlapping genes, with UTS giving the most accurate results. The UTS method discards fewer reads and may be advantageous for experiments with less sequencing depth. UTS is compatible with any aligner and, unlike isoform-aware alignment methods, can be implemented on a laptop computer. Our findings demonstrate that accuracy achieved by complex and expensive techniques such as CAGE-seq can be approximated using conventional short-read RNA-seq data when read assignment methods address transcript overlap. Although our study focuses on EBV transcription, the UTS method should be applicable across all herpesviruses and other genomes with extensively overlapping transcriptomes. IMPORTANCE Many viruses employ extensively overlapping transcript structures. This complexity makes it difficult to quantify gene expression using conventional methods including RNA-seq. Although high-throughput techniques that overcome these limitations exist, they are complex, expensive, and scarce in herpesvirus literature relative to short-read RNA-seq. Here, using Epstein-Barr virus (EBV) as a model, we demonstrate that conventional RNA-seq analysis methods fail to accurately quantify abundance of many overlapping transcripts. We further show that the previously described Unique CoDing Sequence (UCDS) and our Unique TranScript (UTS) methods greatly improve the accuracy of EBV lytic gene measurements obtained from RNA-seq data. The UTS method has the advantages of discarding fewer reads and being implementable on a laptop computer. Although this study focuses on EBV, the UCDS and UTS methods should be applicable across herpesviruses and for other viruses that make extensive use of overlapping transcription.

Reduced IRF4 expression promotes lytic phenotype in Type 2 EBV-infected B cells

“…The heatmap in Fig 1 was generated as follows. EBV gene expression levels were quantified in RPKM using the UTS method [ 86 ]. Hierarchical clustering was calculated and visualized using the Euclidean distance method with the ComplexHeatmap (v2.8.0) R package [ 87 ].…”

Rta is the principal activator of Epstein-Barr virus epithelial lytic transcription

“…In order to classify EBV lytic genes according to our schema of responsiveness (Rta responsive , Rta synergy , and Rta+Zta), we used our recently described UTS method [ 41 ] to quantify expression of each transcript in our different trans-complementation conditions. For each early or leaky transcript with sufficient depth to interpret (see materials and methods for details) we applied the following definitions: Genes that were only expressed in response to trans-complementation with both Rta and Zta were classified as Rta+Zta genes ( Fig 5 ).…”

“…Resulting reads were mapped to the Akata_BAC_GFP genome [37] using STAR (v2.7.6a) with default settings [78]. Reads were assigned to EBV genes using with mmquant and UTS using default settings [41,79] and a previously described EBV1 annotation file [41] except that the LMP2A and LMP2B genes were considered as a single gene whose expression level was PLOS PATHOGENS calculated based on sequencing depth of the common exons (2)(3)(4)(5)(6)(7)(8)(9). The heatmap in Fig 5 was generated using the ComplexHeatmap (v2.8.0) R package [80].…”

A comprehensive single cell data analysis of in lymphoblastoid cells reveals the role of Super-enhancers in maintaining EBV latency

“…Specifically, we identified three prominent clusters that were marked by the expression of the EBV genes GP350 and LMP1/BNLF2 . Of note, the EBV genome has extensive numbers of overlapping genes such as LMP1 and BNLF2a/b , making the quantification of their mRNA expression more challenging (61). This challenge could be further exacerbated by the 3’ mRNA capture bias in some of the current single cell technologies.…”

“…10x Genomics Cell Ranger 6.0.2 count (77) was used to align the raw sequencing reads to a customized human (GRCh38) and EBV (NC_007605, obtained from (61)) hybrid reference genome to generate barcode and UMI counts. Seurat (v4) (78) was applied for the downstream analysis and visualization of the data as following: Genes that were expressed in less than 3 cells were discarded.…”