Utilization of nanopore direct RNA sequencing to analyze viral RNA modifications

Abstract

Modifications on viral RNAs (vRNAs), either genomic RNAs or RNA transcripts, have complex effects on the viral life cycle and cellular responses to viral infection. The advent of Oxford Nanopore Technologies Direct RNA Sequencing provides a new strategy for studying RNA modifications. To this end, multiple computational tools have been developed, but a systemic evaluation of their performance in mapping vRNA modifications is lacking. Here, 10 computational tools were tested using the Sindbis virus (SINV) RNAs isolated from infected mammalian (BHK-21) or mosquito (C6/36) cells, with in vitro-transcribed RNAs serving as modification-free control. Three single-mode approaches were shown to be inapplicable in the viral context, and three out of seven comparative methods required cutoff adjustments to reduce false-positive predictions. Utilizing optimized cutoffs, an integrated analysis of comparative tools suggested that the intersected predictions of Tombo_com and xPore were significantly enriched compared with the background. Consequently, a pipeline integrating Tombo_com and xPore was proposed for vRNA modification detection; the performance of which was supported by N⁶-methyladenosine prediction in severe acute respiratory syndrome coronavirus 2 RNAs using publicly available data. When applied to SINV RNAs, this pipeline revealed more intensive modifications in subgenomic RNAs than in genomic RNAs. Modified uridines were frequently identified, exhibiting substantive overlapping between vRNAs generated in different cell lines. On the other hand, the interpretation of other modifications remained unclear, underlining the limitations of the current computational tools despite their notable potential.IMPORTANCEComputational approaches utilizing Oxford Nanopore Technologies Direct RNA Sequencing data were almost exclusively designed to map eukaryotic epitranscriptomes. Therefore, extra caution must be exercised when using these tools to detect vRNA modifications, as in most cases, vRNA modification profiles should be regarded as unknown epitranscriptomes without prior knowledge. Here, we comprehensively evaluated the performance of 10 computational tools in detecting vRNA modification sites. All tested single-mode methods failed to differentiate native and in vitro-transcribed samples. Using optimized cutoff values, seven tested comparative tools generated very different predictions. An integrated analysis showed significant enrichment of Tombo_com and xPore predictions against the background. A pipeline for vRNA modification detection was proposed accordingly and applied to Sindbis virus RNAs. In conclusion, our study underscores the need for the careful application of computational tools to analyze viral epitranscriptomics. It also offers insights into alphaviral RNA modifications, although further validation is required.

Keywords: RNA modification; Sindbis virus; alphavirus; comparative methods; nanopore.

Fig 1

Data sets and computational tools used for mapping SINV RNA modifications. (A) Read-level data features of full-length viral reads base called using the Dorado model. Baby hamster kidney BHK-21 and Aedes albopictus C6/36 cells were infected with SINV at a multiplicity of infection of 0.1, followed by polyadenylated RNA isolation and ONT DRS. IVT RNA was synthesized from a plasmid encoding the SINV cDNA. Full-length viral reads were extracted from the alignment files. Read-level accuracy, identity, mismatch, insertion, and deletion were calculated using a custom Python script. (B) Base-level data features of full-length viral reads base called using the Dorado model. The confusion matrixes show the frequencies of each base being correctly base called, miscalled, or deleted in individual samples. (C) Classification of modification analysis tools used in the present study.

Fig 2

Evaluation of 10 computational tools in detecting vRNA modifications with default or recommended parameters. (A) Correlations between probability values computed for the IVT and BHK-21 samples at each testable site using single-mode methods. (B) Correlations between probability values computed for the IVT and C6/36 samples. For panels A and B, equal numbers of reads were subsampled from individual data sets, followed by analyses using m6Anet, Nanom6A, and Tombo_de novo; all predictable A sites are included. (C) Workflow of identifying false-positive predictions resulting from comparative approaches. Briefly, full-length IVT reads were divided into two subsets and used as input data sets for comparison method analyses. The default or recommended cutoffs were applied, and the final predictions were regarded as false positives. (D) Volcano-like plot for output from error rate-based comparative method Differr. (E) Volcano-like plots for outputs from current signal-based comparative methods Tombo_com and Xpore. For panels D and E, dashed lines indicate recommended cutoff values; false-positive predictions are labeled in red; all types of modifications are included. FPR, false-positive rate. (F) Comparative tools with a FPR lower than 0.05%.

Fig 3

Cutoff optimization for comparative methods based on unmodified IVT reads. (A) Effect size cutoff values determined by computing 99% confidence intervals (CIs). One-sided 99% CIs were calculated for error-based methods, and two-sided 99% CIs were computed for signal-based approaches. (B) Significance cutoff values decided by computing percentile. The final percentiles were rounded to the nearest integers and used as new significance thresholds.

Fig 4

Integrated analysis of comparative methods. An equal number of reads, 772 in amount, were randomly subsampled from the BHK-21, C6/36, and IVT samples, followed by analyses using comparative methods. (A and B) UpSet plots for modification sites reported by different comparative methods in the BHK-21 or C6/36 samples. Red asterisks indicate significantly enriched intersections from the background. (C and D) Veen diagrams showing the intersections of Tombo_com and xPore predictions. The P-values were calculated by Fisher’s exact test. (E and F) Boxplots showing normalized signal differences at predicted modification sites and adjacent regions. The plots were generated using Tombo. Predicted sites are indicated by black arrows.

Fig 5

Proposal of a vRNA modification detection pipeline and its performance evaluation. (A) Proposed pipeline for detecting vRNA modifications utilizing computational methods. Briefly, native and IVT reads are subsampled to an equal coverage (>500), followed by Tombo_com and xPore analyses and intersection identification. Cutoff values at the coverage depth of 500 are recommended for use. (B) Venn diagram showing the relationship between modification sites detected at different coverage depths while using the same cutoff values. Different numbers of reads (500, 1,000, and 1,500) were randomly subsampled from the IVT and BHK-21 full-length viral reads and analyzed the Tombo_com and xPore pipeline. (C) Detection of modification sites in the SARS-CoV-2 vRNAs using the Tombo_com and xPore pipeline. DRS fats5 data were obtained from reference . (D) Comparison between the computational pipeline outputs and MeRIP data in detecting SARS-CoV-2 vRNA m6A modifications. Modified A sites were extracted from the total prediction result of the Tombo_com and xPore pipeline. MeRIP data were obtained from reference . Caco-2 and Vero indicate MeRIP data measured in Caco-2 and Vero cells, respectively.

Fig 6

Predictions of modification sites in the gRNAs and sgRNAs. As many reads as possible were subsampled from the original sequencing data while maintaining balanced coverage depths for the input dataset pairs, followed by analysis using the Tombo_com and xPore pipeline. (A) Intersections of modification sites in the SINV gRNAs generated in BHK-21 and C6/36 cells. (B) Positions of predicted modification sites in the SINV gRNAs. Read coverage depths of individual analyses are indicated in parentheses. (C) Intersections of modification sites in the SINV sgRNAs, generated in BHK-21 and C6/36 cells. (D) Positions of modification sites in the SINV sgRNAs.