Meta-analytic approach for transcriptome profiling of herpes simplex virus type 1

In this meta-analysis, we re-analysed and compared herpes simplex virus type 1 transcriptomic data generated by eight studies using various short- and long-read sequencing techniques and different library preparation methods. We identified a large number of novel mRNAs, non-coding RNAs and transcript isoforms, and validated many previously published transcripts. Here, we present the most complete HSV-1 transcriptome to date. Furthermore, we also demonstrate that various sequencing techniques, including both cDNA and direct RNA sequencing approaches, are error-prone, which can be circumvented by using integrated approaches. This work draws attention to the need for using multiple sequencing approaches and meta-analyses in transcriptome profiling studies to obtain reliable results.


Introduction
Second-generation short-read sequencing (SRS) technology -launched in the mid-2000s-, has revolutionized both genomic and transcriptomic research because of its ability to sequence millions of nucleic acid fragments simultaneously at a relatively low expenditure per base. Third generation long read sequencing (LRS) approaches have emerged in recent years. Currently, two LRS methods are in use: single-molecule real-time technology and nanopore-based sequencing developed by Pacific Biosciences (PacBio) and by Oxford Nanopore Technologies (ONT), respectively.
LRS can overcome several shortcomings of SRS in transcriptome analysis mainly based on the ability of LRS techniques to read full-length RNA molecules. However, similarly to SRS, LRS techniques often produce spurious transcripts owing to issues such as template switching and mispriming in reverse transcription (RT) and PCR. The major problem is that no efficient bioinformatic tools are currently available to detect these errors. Direct RNA (dRNA) sequencing is considered superior to cDNA sequencing because dRNA sequencing does not involve RT, second strand synthesis and amplification by PCR, which are prone to generate artefacts (however, notably, direct cDNA sequencing without PCR amplification is also possible using both LRS platforms). Nonetheless, dRNA-Seq has also limitations, such as low throughput, 15-30 bases and missing from the transcription start site. Moreover, errors produced by, for example the ligation used for the attachment of adapters, single-strand cDNA formation, or the potential slippage of RNA molecules during their passage across the nanopore as a result of temporary improper functioning of the ratcheting enzyme are also noteworthy drawbacks of dRNA-seq. The low throughput of dRNA-Seq makes both transcript identification and the annotation of nucleic acid sequences at base-pair resolution difficult, which is especially critical in species with large genomes. LRS has already been applied for the transcriptome analysis of various organisms 1,2 , including viruses [3][4][5][6][7][8][9] . This approach has revealed extremely complex transcriptome profiles in every examined species. LRS techniques can be used in analyses that are challenging for SRS approaches, such as the detection of multi-spliced transcripts, parallel transcriptional overlaps, low-abundance transcripts, very long RNA molecules and embedded transcripts, including 5′-truncated ORF-containing mRNAs and non-coding transcripts. A single technique may fail to detect certain transcripts or transcript isoforms, and to precisely map transcript ends or intron boundaries. Additionally, platform-and library preparation-dependent sequencing errors may produce false isoforms. A meta-analysis including multiplatform approaches, such as various LRS and SRS techniques, as well as different auxiliary methods, such as cap selection, and 5′-and 3′-ends mapping can circumvent the aforementioned problems, especially if different library preparation protocols are used. Furthermore, the comparison of various datasets provides a tool for identifying novel transcripts, validating already-described RNA molecules or removing putative transcripts if not confirmed by other techniques.
Herpes simplex virus type 1 (HSV-1) is a member of Alphaherpesvirinae subfamily of the Herpesviridae family. According to estimates by the WHO, more than 3.7 billion people are infected with this virus worldwide 10 . HSV-1 has a 152-kbp linear double-stranded DNA genome, which is transcribed by the host RNA polymerase in a cascade-like manner producing three kinetic classes of transcripts and proteins: immediate-early (IE), early (E), and late (L) 11 . IE genes code for transcription activators required for the expression of E and L genes. The viral E genes primarily specify proteins playing a role in DNA synthesis, whereas L genes encode structural proteins. The identification of HSV-1 transcripts faces an important challenge due to the polycistronic and overlapping nature of viral transcripts. However, polycistronic units of herpesviruses are different from those of bacterial operons, in that only the most upstream genes are translated due to the use of cap-dependent translation initiation 12 . The herpesvirus genes are organized into tandem clusters generating transcripts with co-terminal transcription end sites (TESs). Previous studies have revealed several novel mRNAs, long non-coding RNAs (lncRNAs) 13-17 and microRNAs 18 .

Results
In this study, we employed an integrated approach based on the meta-analysis of the HSV-1 transcriptome data published by Depledge and colleagues (using ONT dRNA-Seq and Illumina RNA-Seq) 19 24 and our laboratory (Tombácz and colleagues using PacBio RSII 25 , as well as Boldogkői et al. 26 , and Tombácz et al. 27 using PacBio Sequel, ONT dRNA-Seq and cDNA sequencing with multiple library preparation methods; Fig. 1, Supplementary Table 1). Our investigations led to the discovery of several novel transcripts, especially of novel multigenic RNA molecules (Fig. 2), and novel splice sites (Figs. 3-5; Tables 1 and  2, and Supplementary Tables 2 and 3). As Figs. 3 and 4 show, a relatively high percentage of introns identified in a study were not detected in other studies, probably due to the varying strictness of criteria used for the annotations. Another possible reason for the large number of unique introns may derive from the variance between the methodologies (e. g. viral titre of infection, virus strain, etc.) used for dataset generation. We note that a large number of unique introns share the splice donor and/or acceptor sites with other introns, which suggests the existence of these splice sites.   25,27 are already multiplatform-based (various combinations of library preparation techniques of Pacific Biosciences RSII and Sequel, and Oxford Nanopore Technologies MinION sequencing). These datasets were compared with the intron datasets generated by Tang et al. 20 and Whisnant et al. 22 . We also used raw sequencing reads from Depledge's direct RNA-Seq study 19 . The obtained data were aligned to the HSV-1 genome and then analysed using LoRTIA. This analysis detected 214 introns. Four large raw Illumina datasets [19][20][21][22]24 were also mapped and reanalysed. Only the introns that were present in at least two independent datasets were accepted and plotted. We obtained 3,848 additional potential introns from this part of the work (see in Supplementary Table 2). (a) Introns identified by Tombácz and colleagues. Altogether, 44.7% of these introns have been validated by the other studies. (b) Introns identified in Depledge and co-workers' dataset using the LoRTIA tool. Our analysis of the raw dRNA-Seq reads detected 309 potential introns, from which 104 were also found in the other studies. The LoRTIA tool did not identify the previously www.nature.com/scientificdata www.nature.com/scientificdata/ published intron within the RNA encoding the fusion protein RL2-UL1 19 ; however, it was verified by the dataset from Tang and colleagues' publication 20 . (c) Introns published by Tang and colleagues. These authors published a large number of potential introns (2352), but only 5% of them were validated in the other datasets. (d) High-coverage introns from Tang and co-workers' publication. 59 out of 2,352 detected introns were identified as highly abundant by the authors. From these 59, only 16 (27%) were detected in at least one of the other three published intron datasets. (e) Introns from Whisnant and colleagues' publication. They have published 79 introns, 84% of which were also found in other datasets. The authors have analysed our previous dataset 22 and found that seven of the eleven published introns are low-abundance isoforms. Therefore, they considered them as unconfirmed. We found and validated five out of these seven introns in our novel dataset, which were also present in Tang  www.nature.com/scientificdata www.nature.com/scientificdata/ Additionally, we confirmed putative RNA molecules and transcript isoforms which were previously unpublished because of inadequate evidence supporting their existence (Supplementary Table 3). This meta-analysis also revealed that practically all HSV-1 genes contain at least one shorter transcript variant with 5′-truncated in-frame ORFs (Fig. 6). Loosening the annotation criteria probably would lead to the identification of truncated genes in every canonical gene. We also identified several fusion genes with relatively long introns spanning across gene boundaries (Supplementary Table 3b). We confirmed the RL2-UL1 and UL52-UL54 fusion transcripts described by Depledge and colleagues 19 but the longer intron of the RL2-UL1 transcript was only present in a very low abundance in the remapped Illumina dataset published by Tang end co-workers 20 . However, it was undetected in the re-mapped, LoRTIA-filtered Depledge-dataset. In most fusion genes, only introns were identified but not precise transcript termini. Additionally, a large number of low-abundance transcript isoforms -including splice and length variants-were detected in this and other studies 28 also. Whether these molecules have functional significance, or are merely the result of transcriptional noise remains to be ascertained. The general functions of embedded and fusion genes are also unknown. This work also revealed longer transcription start site (TSS) isoforms of several RNA molecules (Supplementary Table 3c). For example, we discovered longer TSS variants for the replication-associated RNAs (raRNAs) 29 that overlap OriL or OriS (Fig. 2), and for latency-associated transcripts (LATs) (Fig. 2). The meta-analytic approach is also suitable for the elimination or addition of unconfirmed transcripts into the "putative" category. For example, a minor fraction of 5′-and 3′-truncated RNA molecules sequenced by the PacBio RS II platform 25 , were undetectable by other techniques, therefore they were removed from the latest list of HSV-1 transcripts.

SRS
Direct RNA sequencing is considered to be the golden standard of transcriptome analysis due its apparent lack of errors. However, we demonstrated here that this technique produces a relatively high level of errors. The fact that we could not detect a large number of dRNA introns in either cDNA database (205 introns in Depledge's  www.nature.com/scientificdata www.nature.com/scientificdata/ dataset and a single intron in our dRNA-Seq dataset) indicates lower than expected fidelity rate of dRNA-seq. In this study, the dRNA dataset produced the shortest average intron length. This dissimilarity in the two datasets is explained by the differences in the depths of coverage. However, the most abundant introns were present in databases of both approaches. Our meta-analysis confirmed the existence of an extremely complex meshwork of transcription overlaps (described by Tombacz and co-workers 27 ), which is produced by transcriptional readthroughs between tandem and convergent genes and by the head-to-head overlap between divergent genes. Here, we identified several very long readthrough RNAs, including complex transcripts (containing at least two genes www.nature.com/scientificdata www.nature.com/scientificdata/ in opposite orientations), and transcript isoforms with long 5′-untranslated regions (5′-UTRs) ( Supplementary  Table 3c). Except for most parallel and some convergent overlaps, the majority of transcription readthroughs generate low-abundance transcripts the function of which, if any, is currently unknown. Transcriptional readthroughs might be the by-products of a genome-wide interference mechanism operating via the collision and competition of various elements of the transcription machinery 30 . This hypothesis does not exclude the possibility that the generated RNA strands also have functions of their own. Our comparative study clearly demonstrates the need for multiplatform and meta-analytic approaches for transcriptome profiling to obtain reliable results.
We assembled the sequence of HSV-1 transcripts using ReadConsensus script SeqTools (https://github.com/ moldovannorbert/seqtools) and our previously published LRS data 27 . We found 157 transcripts after removing those with a read depth less than 30 × (Supplementary Table 4) Their sequence consensus can be found in CITE. We note here that de novo or reference-guided transcriptome annotation is more challenging than genome annotation due to the fact that the same DNA region generally codes for multiple RNA isoforms, including splice, TSS and TES variants. Higher read depths and multiple biological replicates are needed for such analyses.

Discussion
In this study, we re-analysed and compared the datasets on HSV-1 transcriptome generated by eight studies [19][20][21][22][23][24][25][26][27] . Here, we provide the most complete transcriptome of HSV-1 to date. We identified a number of novel RNA molecules and transcript isoforms, including intron and length variants. We also confirmed the existence of previously published transcripts. This multiplatform study also identified and confirmed several low-abundance transcripts, such as mono-and multi-spliced transcripts, 5′-truncated mRNAs with short in-frame ORFs, and very long TSS variants, polycistronic and complex transcripts. The functions of these RNA molecules (if any)  Table 3. Data records that were used in this study.
www.nature.com/scientificdata www.nature.com/scientificdata/ have to be demonstrated experimentally. Furthermore, we also demonstrated that various sequencing techniques, including dRNA-Seq, are error-prone, which can be circumvented by using integrated approaches. This study showed that using different reference genomes for mapping, the same transcripts can lead to somewhat different results with respect to the splice sites, especially in SRS. Taken together, employing multiplatform approaches with distinct library preparation methods is especially important in transcriptome research, because of the high error-rate and variance in the results obtained using various library preparation, sequencing and annotation methods. Furthermore, meta-analyses can account for the potential errors derived from using different kits and protocols, as well as from dissimilar work styles and conditions in different laboratories.

Methods
Datasets In this study, several datasets (Depledge et al. 19 , Tombácz et al. 25 24 ) were reanalysed to define the complete HSV-1 transcriptome. The datasets from our laboratory are filtered from data derived from PacBio and ONT cDNA sequencing, and various ONT library preparation approaches including cDNA-, Cap-selected cDNA and dRNA sequencing 26,27 . The wetlab and in silico protocols are detailed in the above mentioned studies. The SRS datasets were used only for the  Table 2). (b) An example for an embedded ORF-containing transcript detected by various techniques. Visualization of the UL2 transcript and one of its truncated transcripts (ul2.5) using Integrative Genomics Viewer. The sequencing reads are from long-read (LRS) sequencing and short-read sequencing (SRS) datasets including direct RNA (dRNA) and cDNA sequencing. It can be seen that the dRNA-seq and the two LRS cDNA techniques detected the same TSS (note that dRNA sequencing produces shorter 5'-UTRs [on average, 23 bp are missing]). The figure also shows that SRS without a specialized library preparation method (e.g., CAGE) is not sufficient to identify 5'-ends of transcripts. (2020) 7:223 | https://doi.org/10.1038/s41597-020-0558-8 www.nature.com/scientificdata www.nature.com/scientificdata/ identification of intron donor and acceptor sites, whereas the LRS data were used to detect novel splice variants and TES-and TSS-isoforms.
Data analysis Detection of introns was carried out by a two-step analysis (Fig. 7). First, we compared published introns 19,20,22,26 with each other, then we remapped the raw data and used them to identify potentially novel introns. We also analysed the effect of the selected reference genome and the aligner on the obtained results [19][20][21][22][23][24][25][26][27] ( Table 3). The adapter sequences from raw reads of each SRS run were removed using the Cutadapt v2.6 software. The fastp tool was used for validation. Next, we aligned the sequencing reads to the HSV-1 reference genome (GenBank: X14112.1) using minimap2 or STAR mapper for the LRS or the SRS data, respectively. The LoRTIA tool (https://github.com/zsolt-balazs/LoRTIA) was used to annotate introns, TSSs, and TESs from the LRS data (Fig. 8); whereas we used the STAR software to detect introns from the SRS samples. The previously published introns (Tang et al. 20 , Wishnant et al. 22 , and Tombácz et al. 25,27 ) were compared with each other, reanalysed, and validated using the datasets from all of the aforementioned publications.
In this work, we assembled the sequence of HSV-1 transcripts with SeqTools/ReadConsensus scripts using our previously published LRS data 27 . The alignment of a transcript's sequencing reads annotated by LoRTA to the section of the reference genome overlapped by the annotation was performed using minimap2. This was followed by variant calling using bcftools' mpileup and call functions, and consensus sequence generation using bcftools' consensus function. Read depth for each transcript was calculated by LoRTIA. To avoid sequencing errors, transcripts with a coverage of less than 30x were eliminated. This read depth is standard for MinION genome assembly [31][32][33] .

Data availability
The datasets used in this work were publicly available and were obtained from the original publications (  26 , and from Tombácz et al. 25,27 . All data generated in this study are included in Supplementary Tables 2 and 3. The data of introns plotted in this study were obtained from Tang et al. 20 , Rutkowski et al. 21 , and from Tombácz et al. 25