Clair3-RNA: a deep learning-based small variant caller for long-read RNA sequencing data

Domain: physics
Published: 2025-12-22
Editorial outcome: major_revision

Abstract
Variant calling with long-read RNA sequencing (lrRNA-seq) helps to analyze full-length isoforms and gene expression but is complicated by high error rates, transcript diversity, RNA editing events, etc. Here, we propose Clair3-RNA, the first deep learning-based variant caller tailored for lrRNA-seq data. Building upon Clair series’ pipelines, Clair3-RNA enhances lrRNA-seq variant calling using optimized techniques, such as uneven coverage normalization, refined training data, editing site discovery, and haplotype phasing to enhance performance. Clair3-RNA supports various platforms, including PacBio, ONT complementary DNA sequencing (cDNA), and ONT direct RNA sequencing (dRNA). Clair3-RNA achieved a ~91% SNP F1-score on the ONT platform using the latest ONT SQK-RNA004 kit (dRNA004) and a ~92% SNP F1-score in PacBio Iso-Seq and MAS-Seq for variants with at least 4x coverage. With least 10x coverage and disregarding zygosity, the performance reached a ~95% and ~96% F1-score for ONT and PacBio, respectively. After phasing, the performance reached ~97% for ONT and ~98% for PacBio. Across GIAB samples, Clair3-RNA consistently outperformed existing callers and accurately identified RNA editing sites. Clair3-RNA is open-source at (https://github.com/HKU-BAL/Clair3-RNA).

Introduction
RNA sequencing (RNA-seq) is extensively utilized for gene expression quantification, isoform identification and quantification, and gene fusion analysis1,2. Recent advancements highlight the potential of variant calling3–6 in RNA-seq, showcasing its numerous advantages for this application. First, it is more cost-effective than DNA sequencing (DNA-seq). Second, the presence of multiple isoforms of some genes increases sequencing depth, enhancing the efficacy of variant detection. Third, RNA-seq allows the identification of post-transcriptional modifications7, which are not detectable using DNA-seq.

Long-read RNA sequencing (lrRNA-seq) has become popular for its capability to capture full-length isoforms8 and novel transcript discovery9. Pacific Biosciences (PacBio) and Oxford Nanopore Technology (ONT) are two leading platforms in this field. PacBio introduced Iso-Seq for full-length cDNA sequencing using PacBio SMRT sequencing technology10. PacBio also provides the advanced Multiplexed Arrays Sequencing (MAS-Seq) method, which increases throughput by combining cDNA molecules into longer concatenated fragments. ONT also introduced direct RNA (dRNA) sequencing, which is a nanopore-based technique that directly sequences RNA without enzymatic conversions11. ONT dRNA sequencing enables the observation of full-length transcripts, RNA modifications, and RNA epigenetics analysis12, which are challenging to detect using cDNA sequencing. Recently, deep-learning methods have improved the capabilities of germline and somatic DNA variant calling13–15 and short-read RNA variant calling3. Although previous works have tried repurposing DNA-seq variant callers for RNA-seq germline variant calling with transformed alignments16, or readapting a short-read RNA-seq variant calling pipeline17 and DNA-seq variant calling workflow3,14 to work with long-read, there still exists no variant caller natively designed for lrRNA-seq that considers RNA-seq specific characteristics including splice junctions and RNA editing in model training.

RNA-seq offers advantages for variant calling and interpretation, such as that the identified allele is guaranteed to be expressed in contrast to DNA variant calling. Nevertheless, several disadvantages need to be considered such as a higher error rate than DNA-seq18,19, which averages a 1–5% error rate. This necessitates robust variant-calling systems for distinguishing true variants from sequencing artifacts. Second, unlike the uniform read coverage typically observed in DNA-seq data4, the coverage is uneven across genomic regions in RNA-seq, and the variable coverage poses challenges for accurate variant calling, particularly in regions characterized by inadequate read support or excessive read coverage that is not accounted for in neural network training. Additionally, RNA-seq faces RNA editing events, especially A-to-I editing catalyzed by adenosine deaminases acting on RNA (ADAR)3,20, leading to A-to-G or T-to-C substitutions that can resemble genuine variants, thus introducing false positives in RNA-seq. Moreover, while phasing has been demonstrated to be advantageous for lrDNA-seq analyses13,14, no studies have yet investigated the impact of phasing on variant calling performance in lrRNA-seq, particularly when reads are shorter, and zygosity flipping can happen.

In this study, we introduce Clair3-RNA, the first deep learning-based small-variant caller designed for long-read RNA-seq data. Clair3-RNA utilizes a Bi-LSTM pileup model with multi-task output for variant prediction. The model was trained on lrRNA-seq datasets using various techniques, such as coverage normalization, variant zygosity flipping, training sample augmentation, and RNA editing-site tagging to optimize variant-calling accuracy. We further developed a workflow that incorporates additional phasing information into the neural network, thereby enhancing variant calling performance and demonstrating the efficacy of phasing for lrRNA-seq. We also designed a benchmarking pipeline that defines callable regions and benchmarking metrics, enabling a fair evaluation of variant-calling performance. Our benchmarking results demonstrate that Clair3-RNA outperformed other variant callers and can achieve ~93%, ~93%, and ~91% F1-score in PacBio Iso-Seq, PacBio MAS-Seq, and the ONT dRNA004 platform, respectively, for variants with at least four supporting reads. With the incorporation of additional pileup phasing features, Clair3-RNA reached 98.59% SNP F1-score for PacBio MAS-Seq, 98.24% for Iso-Seq, and 97.16% for ONT dRNA004 when excluding zygosity and requiring a minimum of 10 supporting reads. Furthermore, in our analysis of GENCODE-annotated genes, we identified 1085 and 507 genes that could be completely phased, as well as 1165 and 655 genes without SNP, Indel or switch errors in the PacBio MAS-Seq HG004 and ONT dRNA004 HG004 datasets, respectively. Furthermore, we established a high-quality benchmark set of RNA editing sites by cross-evaluating RNA samples with corresponding DNA samples. Our results show that Clair3-RNA can identify 5677 and 2549 of editing sites in PacBio and ONT HG004 datasets, respectively. Clair3-RNA also demonstrated outperforming runtime and memory usage compared to other variant callers, requiring ~20 min for PacBio Iso-Seq and ~30 min for ONT dRNA004 variant calling using 16 CPU cores, with a memory footprint of about 7 GB.

Results
Overview of Clair3-RNA
Clair3-RNA is a multi-task bidirectional long short-term memory (Bi-LSTM) neural network designed for lrRNA-seq variant calling and Fig. 1 shows an overview of the Clair3-RNA workflow. Clair3-RNA takes a BAM/CRAM file as input to identify candidate variants with any minor allelic fraction exceeding a specific threshold (typically 0.08, configurable by option) and then computes a summarized pileup tensor of the candidate and employs a neural network for prediction. Clair3-RNA includes two probabilistic output tasks - 21-genotype and zygosity - facilitating the prediction of variants and artifacts. The “Methods” section details techniques for training data augmentation and optimizations specific to lrRNA-seq, addressing challenges like uneven coverage, variant zygosity inconsistency, and RNA editing events. Clair3-RNA was trained with the Genome in a Bottle (GIAB) HG002 sample and supports variant calling on multiple PacBio and ONT instruments and chemistry configurations. Supplementary Fig. 1 illustrates the statistics of the training samples, including the distribution of variant types, output task categories, and allelic fraction.

In this section, we assess the performance of Clair3-RNA on PacBio cDNA, ONT cDNA, and ONT dRNA datasets. Our evaluation involved comparing results with four other variant callers, including LongcallR21 (a long-read SNP caller for RNA data), Clair313 (our deep learning-based germline variant caller for DNA data), DeepVariant14 (the first deep learning-based variant caller, developed by Google), and CTAT-Mutations17 (a variant calling pipeline focused on detecting variants from RNA-seq). For Clair3-RNA and longcallR, we used the corresponding specific RNA model and parameter settings for RNA variant calling. For longcallR, we also evaluated its neural network-based extension, longcallR-nn21 into benchmarking. For DeepVariant, which has extended its functionality to support RNA variant calling in Illumina short-read and PacBio long-read3, we selected the corresponding PacBio model for benchmarking. Clair3 and the other tools in the series15,22–25 before Clair3-RNA were used for germline and somatic small variant calling, but none of them were designed for lrRNA-seq. Therefore, for performance comparison on long-read data, the best available DNA models of Clair3 were employed for benchmarking. For benchmarking, we selected two GIAB samples: HG004 (the mother of HG002) and HG005 (from an unrelated Chinese trio family). This dual-sample approach enables comprehensive evaluation while controlling for potential biases from biological relationship, as HG004 shares genetic background with HG002. Since LongcallR is limited to SNP variant detection, the comparisons with Clair3-RNA were conducted on SNP performance only. For other methods, we benchmarked both SNP and Indel performance.

Generating callable high-confidence benchmarking regions for evaluation
RNA transcripts represent only subsets of the corresponding DNA sequence. Certain regions may have insufficient coverage, resulting in undetectable variants in those regions. To address this, we identified callable high-confidence benchmarking regions by intersecting callable regions with adequate read coverage, with those high-confidence regions defined by GIAB. We also considered the zygosity imbalance in the RNA data. In certain cases, a heterozygous variant in DNA could be either a homozygous variant or a homozygous reference in RNA. This discrepancy arises from the dominant expression of one of the haplotypes. By disregarding zygosity, we aim to understand the impact of zygosity during evaluation better. More details are listed in the “Methods–Benchmarking metrics” section.

We further introduced two key benchmarking factors for performance assessment: read coverage (DP) and allele depth (AD). DP signifies the read coverage of the candidate site meeting minimum mapping quality support and excluding all supplementary and secondary alignments. AD represents the read coverage that supports the true alternate allele at the site. We also considered RNA editing events and created a reliable benchmark set of editing sites for comprehensive evaluation. More details are provided in the “Methods — Orthogonal high-quality RNA editing benchmark set generation” section.

Performance on ONT
The performance on the ONT platform was carried out in six datasets. Details of the benchmarking datasets are provided in Supplementary Table 1. These ONT datasets include ONT cDNA sequencing, dRNA sequencing using the SQK-RNA002 kit (hereafter referred to as dRNA002), and dRNA sequencing employing the ONT latest SQK-RNA004 kit (hereafter dRNA004). The ONT cDNA and dRNA002 data were made available by GIAB. The ONT cDNA data was sequenced using R9.4.1 chemistry and base called using Guppy26 version 6.4.6. The ONT dRNA002 data was sequenced using the dRNA SQK-RNA002 sequencing kit and base called with Guppy 6.4.6. In 2023, ONT released the dRNA SQK-RNA004 kit27 with an enhanced throughput and a reduced error rates compared to the SQK-RNA002 kit. We sequenced GIAB HG002, HG004, and HG005 using the dRNA SQK-RNA004 kit at the University of Hong Kong. The total RNA of three samples was acquired from the Coriell Institute, sequenced using the SQK-RNA004 kit, and base called using the Dorado basecaller28 version 0.7.2. All base called alignments were mapped to the GRCh38 reference genome using the minimap2 aligner in splicing alignment mode. More details can be found in the “Methods—ONT SQK-RNA004 kit library preparation and sequencing” section. As shown in Fig. 2a, the overall throughput increased from an average of 3 GB to 18 GB per flow cell using dRNA004 compared with dRNA002 and the average error rate decreased from 9.7 to 1.8%.

As shown in Fig. 3b and Supplementary Table 2, Clair3-RNA consistently outperformed other variant callers in the ONT dRNA004 and cDNA datasets for both the HG004 and HG005 samples. Specifically, Clair3-RNA achieved 91.00% and 91.73% F1-score in HG004 and HG005, respectively, using dRNA004 sequencing when DP≥4 and AD≥2. In comparison, LongcallR had HG004/61.07% and HG005/65.20% F1-score, respectively; Clair3 had HG004/74.85% and HG005/77.21% F1-score. LongcallR, as explained by its authors, expects higher sequencing coverage for a better accuracy21. When we increased the minimal requirement of DP to 10 and AD to 4, we observed that LongcallR performance improved; the F1-score reached HG004/72.42% and HG005/74.42%, resulting in an average ~10% improvement against the default setting. Using the same setting, while the performance gap narrowed, Clair3-RNA still performed better than longcallR, reaching HG004/94.75% and HG005/95.17%, respectively.

In the cDNA dataset, Clair3-RNA yielded HG004/68.60% and HG005/79.83% F1-score and achieved better performance than longcallR (HG004/+23.72%, HG005/+23.72%) and Clair3 (HG004/+17.85%, HG005/+19.19%). The lower performance across all callers might be because of a higher base error rate of the cDNA dataset generated with R9.4.1 chemistry. We expect the performance to improve notably with cDNA sequencing using R10 chemistry.

Boosting the performance across ONT RNA-specific flow cells. In our evaluation, we compared the performance of different callers in the GIAB HG004 and HG005 dRNA002 and dRNA004 datasets. As illustrated in Fig. 2b and c, the data throughput increased ~six times per flow cell, from ~3 GB in dRNA002 to ~18 GB in dRNA004. The read N50 also increased from ~1100 bp in dRNA002 to ~1500 bp in dRNA004. Across all callers, we observed improvements in precision, recall, and F1-score, signifying improved read accuracy in the dRNA004 data. Notably, Clair3-RNA exhibited an average F1-score improvement of 22.07% (rising from 68.61 to 91.00% in HG004 and from 69.98 to 91.73% in HG005), where longcallR improved by 16.35% and Clair3 improved by 48.69%.

Among the three ONT benchmarking datasets, we noted that Clair3 had lower precision and longcallR had lower recall. In the dRNA004 dataset, Clair3 had only HG004/66.93% and HG005/70.33% for precision, but HG004/84.89% and HG005/85.58% for recall. LongcallR had better precision than its recall, with HG004/94.49% and HG005/94.89% precision, but HG004/45.11% and HG005/49.67% recall. Compared with other callers, Clair3-RNA demonstrated balanced precision and recall, underscoring the necessity of a variant caller specific for lrRNA-seq.

Performance on PacBio
We benchmarked the performance on the PacBio platform using two sequencing technologies, Iso-Seq and MAS-Seq, both developed for PacBio cDNA sequencing. We collected the HG004 and HG005 Iso-Seq RNA sequencing datasets from two sources, Baylor College of Medicine (BCM) and Google, and the HG004 MAS-Seq RNA sequencing dataset from BCM. More details of sequencing data are provided in Supplementary Table 1. Compared with ONT sequencing, PacBio MAS-Seq and Iso-Seq have a higher N50 read length, with ~3000 bp for PacBio Iso-Seq and ~1700 bp for MAS-Seq, versus ~1200 bp for ONT sequencing.

As indicated in Fig. 4 and Supplementary Table 3, Clair3-RNA consistently achieved a SNP F1-score exceeding 91% across all Iso-Seq and MAS-Seq datasets for HG004 and HG005, showcasing better performance than the second-best DeepVariant RNA model, with an averaged 2.56% F1-score improvement, and outperforming other variant callers including longcallR, longcallR-nn, and CTAT-Mutations. Clair3-RNA also demonstrated consistent performance across all PacBio sequencing sources, achieving HG004/94.93% and HG005/95.50% F1-score for BCM Iso-Seq and an HG004/93.38% F1-score for BCM MAS-Seq. In terms of precision and recall balance, Clair3-RNA exhibited a balanced performance with 95.96% precision and 93.29% recall. Notably, DeepVariant achieved an average 91.98% precision and 92.19% recall, and CTAT-Mutations has high recall (92.58%) while inferior precision (76.67%). Clair3 maintained an average recall rate of 95.38%, but had a low precision (70.30%). Conversely, LongcallR had 93.57% precision but 68.18% recall. In addition to whole-genome sequencing benchmarking, we also evaluated performance on chromosome 20—which was held out from training for both Clair3-RNA and DeepVariant RNA model. As shown in Supplementary Table 5, Clair3-RNA consistently outperformed the DeepVariant RNA model by an average 3.26% F1-score in SNP and 5.51% F1-score in Indel across all read coverages on PacBio MAS-Seq chr20 datasets.

The PacBio datasets were originally aligned with the GRCh38 reference genome using pbmm2, a minimap2 aligner wrapper specifically designed for the PacBio platform. During our experimental analysis focusing on alignment methods, we observed that the pbmm2 aligner tended to introduce more false positives and false negatives compared to the minimap2 aligner when using Clair3-RNA for small variant calling on PacBio. We aligned five datasets using pbmm2 and minimap2 splice mode to assess the impact of different alignment strategies. Further details regarding the alignment procedure are included in the Supplementary Notes - Command line used section.

As shown in Fig. 5 and Supplementary Table 6, when utilizing minimap2 in place of pbmm2, Clair3-RNA’s performance improved across all datasets. The average number of false positives and false negatives decreased by 25.42% and 0.79%, respectively. The averaged F1-score of Clair3-RNA improved from 94.57 to 95.10% (+0.53%) for BCM Iso-Seq, from 91.37 to 91.98% (+0.61%) for Google Iso-Seq, and from 92.76 to 93.38% (+0.62%) for BCM MAS-Seq. Similar trends were observed in other callers; longcallR and Clair3 demonstrated average F1-score improvements of 0.46% and 14.40%, respectively. These improvements suggest that using minimap2 had better performance compared to using pbmm2 when using the default splice mode setting in PacBio.

A recent work introduced a method that repurposes long-read DNA variant callers to use with transformed RNA-seq alignments for long-read RNA-seq based variant calling16. We have also benchmarked against the method. The workflow and command lines used for alignment transformation and repurposing Clair3 (Clair3-mix) and DeepVariant (DeepVariant-mix) are detailed in the Supplementary Notes - Command line used section. The results, presented in Supplementary Table 7, revealed that both the repurposed Clair3-mix and DeepVariant-mix exhibited improved performance on RNA-seq when used with the transformed alignments. However, Clair3-RNA still outperformed both Clair3-mix and DeepVariant-mix across all benchmarked datasets, highlighting the advantage of variant callers naively designed for lrRNA-seq.

Sequencing coverage and allele depth analysis
Insufficient coverage and inadequate read support for alternate bases may lead to false negatives in evaluations. Supplementary Table 8 shows the size and percentage of genomic regions at certain read coverages (DP) or above. Figure 5a shows the benchmarking performance across different DP and AD thresholds. In the ONT dRNA004 HG004 dataset, increasing the depth from 4 to 8 resulted in a 1.54% increase in precision and a 2.67% improvement in recall. Raising the depth to 10 showed marginal gains of +0.33% in precision and +0.58% in recall, highlighting the difficulty in capturing low-coverage variants. Similar trends were observed in the PacBio datasets. In the BCM Iso-Seq HG004 dataset, raising the depth benchmark from 4 to 8 led to precision and recall improvements of 2.15% and 1.90%, respectively. Increasing to a depth of 10 showed modest gains of 0.38% in precision and 0.18% in recall. AD is another crucial factor for variant evaluation. When the AD was increased from 2 to 4, the ONT dRNA004 HG004 dataset showed a 0.28% improvement in precision and a 1.91% increase in recall. Similarly, in the BCM Iso-Seq HG004 dataset, there was an improvement of 0.51% in precision and 1.51% in recall. This highlights the importance of allele depth in accurately calling missing variants.

Zygosity flipping analysis in RNA dataset
Zygosity flipping (an observation in variant calling due to the expression of only one haplotype) is also common in RNA data due to factors like differential gene expression and RNA splicing3. To evaluate this, we used the Clair3-RNA’s “calculate_overall_metrics” function to consider the zygosity flipping by disregarding zygosity in benchmarking. As shown in Fig. 5a, Supplementary Table 2, and Supplementary Table 3, disregarding zygosity led to an average 2.43% and 2.88% improvement in the F1-score on the ONT and PacBio, respectively. Benchmarking on the ONT dRNA004 HG004 dataset, disregarding zygosity resulted in a reduction of 2699 FPs and 2698 FNs, enhancing the F1-score by 2.76%. Similarly, in the BCM Iso-Seq HG004 dataset, the number of FPs decreased by 4334 and FNs by 4331, leading to an improvement of 2.69% in the F1-score. Moreover, the impact of disregarding zygosity lessened as DP increased. For instance, in the ONT dRNA HG004 dataset, setting DP≥8 while considering zygosity boosted the F1-score by 4.19% over DP≥4, whereas disregarding genotyping resulted in a 2.45% F1-score improvement. These results emphasize the importance of read coverage in variant detection, especially when incorrect genotypes are introduced without adequate reads to determine variant zygosities.

Performance analysis by genomic context
We also carried out benchmarking in challenging genomic regions. Specifically, we compared the SNP performance of Clair3-RNA and LongcallR in areas defined as difficult by GIAB Stratification V329. We categorized the analysis into five groups: “Low complexity”, “Segmental duplications”, “Low mappability”, “Functional regions”, and “Other difficult regions” - and conducted it on ONT dRNA004 HG004, PacBio BCM Iso-Seq HG004, and PacBio MAS-Seq HG004 datasets. As shown in Fig. 5b and Supplementary Table 9, Clair3-RNA consistently outperformed LongcallR in SNP F1-score.

Benchmarking on the ONT dRNA004 dataset, Clair3-RNA demonstrated varying performance, achieving an average 80% F1-score. Notably, Clair3-RNA excelled in coding regions, attaining a 94.58% F1-score in the Coding Sequence (CDS) regions. In contrast, its performance dropped to 61.80% in low complexity regions characterized by homopolymer regions greater than 11 bp. In comparison, LongcallR exhibited more significant variability on the same dataset, peaking at 82.51% in the CDS regions, while dropping to 29.85% in homopolymer regions greater than 11 bp. Similar trends were observed in the PacBio Iso-Seq and MAS-Seq datasets, Clair3-RNA’s F1-score averaged ~86%, with the highest at 95.22% in the CDS regions and the lowest at 74.76% in the homopolymer regions. In this dataset, LongcallR averaged about 66%, with a peak of 91.74% in the CDS regions and a low of 37.07% in the homopolymer regions. Among all categories, “Low complexity” and “Segmental duplications” were the most challenging for variant detection, with the lowest average F1-score of all categories.

Indel performance
In our evaluation, we compared the performance of various callers, excluding LongcallR, which identifies only SNP variants. For the ONT platform, as shown in Supplementary Table 2, using DP≥10 and AD≥4 coverage settings, Clair3-RNA achieved HG004/66.50% and HG005/72.37% F1-score for the dRNA004 dataset, and HG004/50.63% and HG005/60.88% for the cDNA dataset. It consistently outperformed Clair3 (dRNA004 HG004/52.65% and HG005/60.64%, cDNA HG004/45.63% and HG005/55.04%) in both the cDNA and dRNA004 datasets. In a comparison between the dRNA004 and dRNA002 datasets, Clair3-RNA demonstrated a significant performance improvement of 43.11%, with the average F1-score increasing from 26.33 to 69.44%. Even in the most precise dRNA004 dataset, the averaged indel recall was lower than 60%, attributed primarily to the presence of indel signals with varying lengths clustered at the same site, confusing the model from calling Indel precisely.

On the PacBio platform, as shown in Supplementary Table 3, using the DP≥10 and AD≥4 coverage setting, Clair3-RNA achieved an Indel F1-score of HG004/89.42% and HG005/92.90% for BCM Iso-Seq, and HG004/80.05% for BCM MAS-Seq. With an average Indel performance of 87.46%, Clair3-RNA outperformed Clair3 and DeepVariant RNA model, where the average F1-score was 69.13% (+18.33%) and 81.07% (+6.39%), respectively.

RNA editing tagging using the REDIportal database. The REDIportal has compiled a total of 36,234,771 RNA editing sites, categorized as follows: (1) 13,101,296 editing sites collected only in ATLAS (A); (2) 22,836,582 editing sites collected by ATLAS and RADAR (A&R); (3) 4,708 editing sites collected by ATLAS and DARNED (A&D), and (4) 292,185 editing sites reported by all three sources (A&R&D). We tagged the variant as an RNA editing site that was validated by at least two independent sources within the REDIportal, specifically in ATLAS, DARNED, and RADAR3,30,31. As shown in Supplementary Table 4, excluding the variants tagged by the REDIportal database, there was enhanced SNP performance across all platforms. Specifically, Clair3-RNA achieved a 0.82% improvement in the ONT dRNA004 HG004 dataset, with a total of 1898 editing sites tagged. For the PacBio datasets, the improvement was less significant, with 1033 and 912 editing sites tagged in the BCM Iso-Seq HG004 and MAS-Seq HG004 datasets, respectively, resulting in a 0.28% and 0.13% F1-score improvement.

A-to-G and T-to-C false positives analysis
A-to-I RNA editing poses a significant challenge for RNA variant calling algorithms because they must distinguish between A-to-T and T-to-C SNPs induced by editing events and genuine A-T and T-to-C SNPs. To address this issue, we analyzed A-to-G and T-C SNPs identified by Clair3-RNA and Clair3 that were not present in the GIAB ground truth. The high-quality benchmark set of the editing sites were established and analyzed using the PacBio HG004 BCM Iso-Seq dataset and the ONT HG004 dRNA004 dataset. The dRNA004 dataset was aligned with the HG004 DNA dataset, sequenced by ONT R10.4.1 Q20+ chemistries, and base called using Dorado, provided by Nanopore EPI2ME Labs32. For the BCM HG004 dataset, we compared the RNA dataset with the PacBio HiFi HG004 DNA dataset generated by Google33. We then identified potential A-to-I editing sites by looking for evident A-to-T or T-to-C substitutions in the RNA data that were absent in the DNA data. More details of the benchmarking editing site generation are described in the “Methods—Orthogonal high-quality RNA editing benchmark set generation” section. Finally, we identified 3800 and 6086 editing sites on the ONT dRNA004 and PacBio Iso-Seq datasets, respectively. 46.21% (1756 over 3800) and 55.71% (3391 over 6086) of the benchmark set of editing sites were found in the REDIportal database.

In the Venn diagrams presented in Fig. 6b, the A-to-G and T-to-C false positives identified by Clair3-RNA and Clair3 were extracted for comparison. Clair3-RNA reported 3940 false positives in ONT and 3979 in PacBio, significantly fewer than Clair3, which identified 30,411 false positives in ONT and 94,476 in PacBio. Among these false positives, 1251 (539+712) and 409 (81+328) variants intersected with the benchmark set of editing sites for ONT and PacBio, respectively, indicating that Clair3-RNA can classify 2549 (1505+1044 over 3800) and 5677 (2614+3063) high-quality RNA editing sites as non-variants through network prediction. In comparison, Clair3 detected only 1592 and 1011 RNA editing sites for ONT and PacBio, respectively. The results indicate that Clair3-RNA demonstrates a capability in rejecting RNA editing sites during variant detection. As shown in Fig. 6c, for the remaining false positives in Clair3-RNA, 29.87% (1177 over 3940) in ONT and 26.94% (1072 over 3979) in PacBio were annotated by the REDIportal database. After network prediction and REDIportal database tagging, 62.77% of the false positives in the PacBio dataset and 38.37% of the false positives in ONT were identified. After network prediction and editing site tagging by the REDIportal database, the remaining A-to-G and T-to-C false positives can be reduced from 27,175 to 2051 on ONT and from 54,942 to 2579 on the PacBio platform. The results suggest the feasibility of utilizing the REDIportal database to annotate challenging editing sites that are missed by network prediction.

Editing sites comparison and analysis
We identified the high-quality editing sites across additional datasets and analyzed their intersections, as depicted in Fig. 6d. In the Iso-Seq (Google) dataset, we identified 8049 sites, while the MAS-Seq dataset yielded 15,659 sites. It is not noteworthy that the ONT dRNA004 dataset contained fewer sites compared to the PacBio datasets, primarily due to the absence of a reverse transcription step in the ONT protocol. This is in contrast to PacBio’s cDNA preparation, which includes the reverse transcription step as shown in Supplementary Fig. 3. Additionally, the MAS-Seq dataset exhibited a significantly higher number of sites than the other datasets, which can be attributed to its broader genomic region coverage. This resulted in a greater number of detected variants on the MAS-Seq dataset. Among the four datasets, the PacBio datasets demonstrated a higher degree of overlap with each other than with the ONT dataset, highlighting the distinct methodological differences between ONT dRNA and PacBio cDNA sequencing approaches.

We used the vcf2maf tool34 to annotate the variants identified by Clair3-RNA and editing sites in the PacBio Iso-Seq HG004 and ONT dRNA004 HG004 datasets. As shown in Fig. 6e and Supplementary Table 10, both the variants and editing sites were found to be predominantly located in intronic regions. About half of the PacBio editing sites, ONT variants, and ONT editing sites were in introns, while over three-quarters of PacBio variants were in intron. Editing sites, however, occur more often in splice sites and splice regions. Specifically, 27.43% of the editing sites were in these regions, compared to just 0.19% of the variants. In the ONT dRNA004 dataset, 5.50% of the editing sites but 0.53% of the variants, were found in splice sites and splice regions. These findings align with research on the A-to-I editing survey20,35.

Enhancing variant calling performance with phasing information
A heterozygous variant identified in a single haplotype is considered more reliable than sequencing errors that affect both haplotypes13. Based on this, we attempted to incorporate phasing haplotype information into the neural network, with more details provided in the “Methods—Integrating phasing information into neural network” section. In the variant phasing phase, as illustrated in Fig. 7d, we analyzed phased heterozygous variants and found that at least 63.34% of heterozygous variants can be phased. The percentage rises to 78.05% and 79.07% when using the latest PacBio MAS-Seq HG004 and ONT dRNA004 HG004 datasets, respectively. These datasets have average read lengths of 1731 and 1404 for variant phasing.

For read haplotagging, we computed the percentage of haplotagged reads against all sequencing reads, while excluding supplementary reads, secondary reads, and reads with a mapping quality lower than 5, as recommended by WhatsHap36. Figure 7b and c display the haplotagged percentages of all datasets and the number of haplotagged reads - for PacBio MAS-Seq HG004 and ONT dRNA004 HG004 datasets. The ONT dRNA004 achieved an average haplotagging rate of 26.17%, with PacBio MAS-Seq 35.83%, and Iso-Seq 38.22%. The haplotagging rates of lrRNA-seq are significantly higher than those of the short-read datasets, which is typically less than 10%. This suggests the potential for enhancing variant detection by incorporating phasing information.

We trained the Clair3-RNA pileup model with phasing information from both the PacBio and ONT platforms using the GIAB HG002 dataset. Subsequently, we evaluated the trained models on HG004 and HG005 datasets. As shown in Fig. 7e and Supplementary Table 11, both SNP and Indel precision and recall improved across all datasets. Benchmarking results for the ONT and PacBio datasets indicate an overall improvement of 1.52% in the F1-score for SNP. Specifically, using the PacBio MAS-Seq HG004 dataset, the F1-score increased by 1.23% (94.61% versus 93.38%). The improvement was more significant in the ONT dRNA004 dataset, with an average 2.01% F1-score increase. For HG004, the F1-score increased from 91.73 to 93.56%, and for HG005, it increased from 91.00 to 93.19%. With phasing information, the F1-score reached 98.59% for PacBio MAS-Seq, 98.24% for Iso-Seq, and 97.16% for ONT dRNA004, with at least 10 supporting reads and disregarding zygosity.

The improvement in Indel detection is more significant compared to SNP, achieving an average 5.86% F1-score improvement across all benchmarking datasets. PacBio Iso-Seq demonstrates greater improvement than MAS-Seq (6.54% versus 5.01%) and ONT dRNA004 (6.54% versus 4.92%). The Indel performance reached 85.61% F1-score in the PacBio Iso-Seq HG005 dataset, representing a 4.21% improvement. ONT dRNA004 HG005 achieved a 70.63% F1-score with a 4.58% improvement. By increasing the minimum thresholds to DP 10, AD 4, and disregarding zygosity, the average Indel F1-score can reach 89.42% for PacBio and 72.56% for ONT. Notably, the highest Indel F1-score achieved was 95.17% for PacBio Iso-Seq HG005 and 75.51% for ONT dRNA004 HG005.

GENCODE annotated gene analysis
We evaluated the Clair3-RNA variant calls across gene regions to assess its applicability for functional studies in lrRNA-seq. We used the annotated genes from the GENCODE (version 46) database37. For gene region analysis, we selected those genes that overlapped the defined callable regions specified in the “Methods—Benchmarking metrics” section.

We first categorized all callable genes into three distinct groups based on their overlap with callable high-confidence regions. This overlap ratio H was computed how many percent of a gene is in callable high-confidence regions, defined as the intersected regions of the callable regions and GIAB high-confidence regions. More details are given in the “Methods—Analysis of overlap and error statistics for GENCODE annotated gene” section. Genes were classified as wholly (H=100%), partially (75%<H<100%), or not (H=0%). We then analyzed how many genes are with variants covered and has no calling mistake in each dataset. Furthermore, we compared these genes to the phasesets computed by WhatsHap to determine how many of them are wholly phased.

Figure 8a presents the number of genes covered by callable high-confidence regions for both the PacBio MAS-Seq HG004 and ONT dRNA004 HG004 datasets. For PacBio, a total of 3859 genes were fully covered, 3589 were partially covered, and 1135 were not included in defined callable high-confidence regions. In contrast, the ONT platform showed 2199 fully covered genes, 1870 partially covered, and 876 not included in callable high-confidence regions. Further analysis revealed the number of genes without SNP error. In the PacBio dataset, 1279 genes (33.1% of fully covered) and 2404 genes (67.0% of partially covered) were without SNP error. For the ONT platform, 735 genes (33.4% of fully covered) and 1167 genes (62.4% of partially covered) were without SNP error. Note that the SNP errors in the “not in callable high-confidence region” category could not be determined, because they fall outside the GIAB high-confidence region.

We also summarized the numbers by looking into only those genes annotated as “protein_coding” in the GENCODE database, which are more important for studying protein synthesis and biological functions. As shown in Fig. 8b, we identified a total of 449 and 375 protein-coding genes in PacBio and ONT, respectively. Among these, 242 and 186 genes were detected without SNP error, while 85 and 49 genes were wholly phased in PacBio and ONT. When examining only the coding sequence (CDS) regions, we found that 263 genes (58.6% of the protein-coding genes) in PacBio and 227 genes (58.7%) in ONT were free of SNP errors. Considering only the exon regions, 242 genes (53.9% of the protein-coding genes) in PacBio and 220 genes (58.7%) in ONT had no SNP errors. We further restricted our analysis to genes without SNPs, Indels, or switch errors. In this context, we observed that 261 genes (58.1%) in CDS and 251 genes (55.9%) in exons for PacBio, as well as 227 genes (60.0%) in CDS and 214 genes (57.1%) in exons for ONT, met this criterion. These results indicate a promising variant calling accuracy in the protein-coding regions.

Figure 8c presents the number of wholly phased genes across all benchmarking datasets. The MAS-Seq HG004 identified 1085 wholly phased genes, exceeding 386 and 304 genes reported by the Iso-Seq HG004 and HG005, respectively. Similarly, ONT dRNA004 identified 504 and 479 wholly phased genes in HG004 and HG005, respectively, compared to 211 and 317 in the cDNA datasets and 58 and 52 in the dRNA002 datasets. Figure 8d illustrates the number of genes free from SNPs, Indels, or switch errors, with MAS-Seq leading at 1165 genes, followed by dRNA004 with 655 in HG004 and 653 in HG005. The substantial numbers of wholly phased and error-free genes highlight the potential for reliable downstream gene analysis, particularly using the PacBio MAS-Seq and ONT dRNA004 technologies.

Performance analysis at splice junctions
Variant calling using RNA-Seq is prone to errors at splice junctions16,38,39, where read coverage is often shallow and may be confounded by sequencing artifacts. To systematically assess Clair3-RNA’s performance in these regions, we used the splice junctions annotated in the GENCODE database37. We calculated the variant calling performance at the positions that are 2 to 16 base-pairs from the nearest splice junction boundary (either the start or end of an exon), with 2 base-pairs as a step (resulting in eight groups of results). The results are shown in Supplementary Fig. 5. Across all benchmarking datasets, the F1-score gradually dropped when approaching the splice junction boundary. Specifically, the average F1-score decreased from 90.31% to 83.98% for PacBio Iso-Seq, from 92.52 to 85.45% for MAS-Seq, and from 91.45 to 85.30% for ONT dRNA004, using the best-performing Clair3-RNA model with phasing enabled. These findings underscore the inherent challenges of variant calling in splice junctions. The decline is likely because of relatively fewer variants at or near the splice sites for model training and the accumulation of sequencing error in the truncated exon regions that lead to less confident model predictions.

Performance analysis after removing shared variants and artifacts between samples
We conducted additional benchmarking to evaluate variants and artifacts in HG004 and HG005 while explicitly excluding those present in HG002, which were used during model training, to prevent any potential evaluation bias. As shown in Supplementary Table 12, this exclusion removed ~350k–~651k candidates from HG004 or HG005 benchmarking, depending on the read coverage of the sequencing datasets. After this adjustment, Clair3-RNA achieved with average SNP F1-scores of 94.50%, 90.80%, and 90.43% in PacBio Iso-Seq, MAS-Seq, and ONT dRNA004 datasets, respectively. The model showed further improvement when phasing was enabled, achieving 95.74%, 93.10%, and 92.89% on the same datasets. After removing HG002 variants, the precision remains comparable to the original HG004/HG005 benchmarks in all callable regions (average 95.31% versus 94.77% before phasing, and 96.24% versus 95.61% after phasing). The recall decreased by an average of 4.14% and 2.72% before and after phasing. This reduction is largely attributable to the exclusion of variants overlapping with the HG002 benchmark, which resulted in an average ~43% decrease in the number of benchmarked true positives. Since these overlapping variants are likely enriched for well-supported, common variants that are easier to detect, their removal naturally leads to a modest reduction in recall. Importantly, the precision and recall for the non-overlapping variants remain high, indicating that Clair3-RNA did not overfit to specific samples.

Benchmarking runtime and memory usage
As shown in Fig. 4e and Supplementary Table 13, we evaluated the runtime and memory usage of the Clair3-RNA algorithm and other variant callers. Clair3-RNA demonstrated an outperforming runtime against other variant callers. Specifically, using 16 CPU cores, Clair3-RNA required ~18 min and ~50 min, on HG004 PacBio Iso-Seq and MAS-Seq, and ~34 min and ~29 min on HG004 ONT cDNA and dRNA004 dataset, respectively. This is generally faster than Clair3, LongcallR, and DeepVariant, in which the average runtime ranges from 42 to 80 min. Clair3-RNA also showcased uniform memory consumption across all datasets, with a footprint of ~7 GB. In contrast, DeepVariant required about 26 GB and 24 GB of memory for the PacBio and ONT, respectively. Clair3 displayed a slightly higher memory footprint, on average 49 GB, due to its variant phasing process. In comparison, LongcallR recorded the highest memory consumption at PacBio MAS-Seq, with 83 GB, and averaged 47 GB.

Methods
Overview of Clair3-RNA
Figure 8 shows an overview of the Clair3-RNA workflow. Starting from the read alignments in BAM/CRAM format of an input sample, Clair3-RNA identifies potential variants with adequate read support and base alteration. Leveraging the pileup methodology inherited from Clair3, the candidates and their corresponding flanking window sites are aggregated to create pileup tensors, which are then processed through a Bi-LSTM-based pileup neural network for inference. The neural network produces probabilistic multi-task output that differentiates variants from editing sites and sequencing artifacts. Identified variants are further refined by the editing tagging module to tag known editing sites. Clair3-RNA is compatible with various platforms, including PacBio Sequel with the Iso-Seq kit, PacBio Revio with the MAS-Seq kit, ONT cDNA, ONT dRNA002, and ONT dRNA004.

Clair3-RNA’s input and output
Input. Clair3-RNA takes the input, which includes two components: (1) a reference genome of the sequencing sample, and (2) a sequencing alignment with the same reference material as (1). Clair3-RNA automates the whole variant calling process to get the VCF output, as illustrated below. The pileup tensor consists of 594 integers, covering 33 genome positions, which include the candidate position and 16 flanking positions on both the left and right sides, with 18 features at each position. These 18 pileup features are designated as “A+“, “C+“, “G+“, “T+“, “I+S“, “I+R“, “D+S“, “D+1S“, “D+R“, “A-“, “C-“, “G-“, “T-“, “I-S“, “I-R“, “D-S“, “D-1S“, and “D-R“. A, C, G, T, I, and D represent the counts of read support for the four nucleotides, insertion and deletion, while the symbols “+“ and “-“ represent the read from the forward and reverse strands, respectively. The superscripts “S” and “R” represent the start and remaining positions of Indel. A detailed description of the 18 pileup features is provided in the Supplementary Notes - Description of RNA pileup input features section. All the features above worked with sites with the reads covered, and we assigned a background intensity of 0 for all padding sites. This technique of zero-padding is also utilized in RNA-seq to account for surrounding RNA splicing regions.

Clair3-RNA model architecture. Clair3-RNA retains the bidirectional long short-term memory (Bi-LSTM) model architecture and output tasks from Clair3’s pileup model. The pileup network incorporates two Bi-LSTM layers and three fully connected layers. The output of the pileup model encompasses two probabilistic tasks: 21-genotype and zygosity. Further details about the output tasks are provided in the Supplementary Notes - Description of pileup network outputs section. The identified variants, encompassing both nucleotide alterations and genotypes, were determined by evaluating the combined probabilities from two output tasks. In the case of an Indel variant, the network specifies the variant type, and the length of the Indel is inferred from the supporting read evidence. Additionally, Clair3-RNA incorporates several successful deep-learning techniques from Clair3 for model training, including the (1) model initialization; (2) activation function; (3) optimizer; (4) dropout; (5) L2 regularization; and (6) combination of multiple coverages for model training, and (7) focal loss function for addressing the unbalanced variant type problem.

VCF output. Three categories of variants are shown in Clair3-RNA’s VCF output: (1) variants labeled as “PASS” that passed quality score threshold, (2) labeled as “LowQual” that has low-quality scores (i.e., QUAL<2 for PacBio and QUAL<8 for ONT, configurable), and (3) tagged as “RNA Editing” corresponded to known editing sites identified by the REDIportal database. Specifically, variants were tagged as “RNA Editing” if they match both the reference and alternate alleles given in REDIportal. For each variant entry, the reference, alternate variant, genotype, coverage of the reference allele, and coverage of all alternate alleles are shown.

Training techniques tailored for RNA data
Training candidate generation. In DNA model training, the ground truths were established using the golden GIAB truth variants. However, in RNA model training, owing to limited sequencing coverage in specific regions, not all truth variants are suitable for training. Initially, we scanned the training dataset to acquire the read coverage (DP) and the alternate allele (AD) supporting truth variants, focusing on evident variants with adequate read support. For the artifacts, we scanned the alignments site-by-site to identify candidates with DP≥4, AD≥2, and exceeding a minor allelic fraction (typically 0.08, which was manually decided). Additionally, we accounted for zygosity flipping and editing sites in the training data, as elaborated in the “Methods—Variant zygosity flipping” section.

Variant zygosity flipping. Zygosity flipping (an observation in variant calling due to the expression of only one haplotype) in RNA data can lead to misinterpretations of genetic variants. This occurs when a heterozygous variant appears as a homozygous reference or a homozygous variant due to allelic imbalance. Zygosity flipping in RNA data might be due to differential gene expression when one allele dominates in specific regions. Additionally, RNA splicing might create isoforms, altering the variant zygosity across different isoforms5. Ignoring zygosity flipping in RNA data can have a significant impact on model training. Therefore, before training the model, we computed the alternate allelic fraction of each candidate and adjusted the zygosity obtained from GIAB based on DNA-seq samples. By defining each observed candidate allelic fraction as VAFO and the truth zygosity as GTTruth, where GT∈{0/1, 1/1, 1/2}, we can determine the inferred zygosity GTInferred, as follows:

GTInferred = 
{0/1 if VAFO > VAFHET - ϵ and VAFO < VAFHET + ϵ and GTTruth ≠ 0/1
1/1 if VAFO > VAFHOM - ϵ and GTTruth ≠ 1/1
0/0 if VAFO < VAFArtifact else GTTruth}

where GTTruth is the observed ground truth genotype in GIAB truths. VAFHET, VAFHOM, and VAFArtifact represent the expected allelic fraction for heterozygous variants, homozygous variants, and sequencing artifacts were set as 0.5, 1.0, and 0.08, respectively. ϵ represents the tolerance epsilon parameter and was manually set as 0.2. For multi-allelic 1/2 genotypes, we considered each allele and inferred the genotype when any of the alternate alleles was missing. Finally, we used GTInferred instead of GTTruth for constructing zygosity and 21-genotype labels during model training.

Coverage normalization for extremely high-coverage data
RNA data exhibits an uneven distribution of sequencing read coverage, stemming primarily from variations in transcript abundance, non-uniform RNA expression, and RNA splicing events. These factors result in uneven coverage compared to homogeneous DNA-seq data. The presence of high-coverage sites poses challenges for model inference, as the training data may not adequately represent the full spectrum of coverage distribution. Although Clair3-RNA has implemented techniques like Focal Loss40 to address data imbalances during model training, excessively high coverage levels can significantly impede the accuracy of model inference.

One potential strategy involves subsampling reads during inference to mitigate the gap; this approach, however, may impact performance due to less alignment included in inferences. Building upon this, we proposed a coverage normalization strategy. Initially, we computed the coverage distribution of the training sample and excluded the top 10% of high-coverage sites, identifying the maximum coverage confidence interval as CT. In the inference phase, we determined the coverage of each candidate tensor CI by computing the reads supporting the nucleotide A, C, G, T, and deletions. Subsequently, for each candidate coverage CI exceeding CT, we computed the normalized pileup features Fpileup_norm using the formula:

Fpileup_norm = Fpileup * (CT / CI)

where Fpileup is the total 594 pileup integers used for prediction. In our experiment, using the normalized feature generally improved recall by ~1%, effectively mitigating the influence of extremely high coverage on the model’s inference.

Orthogonal high-quality RNA editing benchmark set generation
RNA A-to-I editing is a post-transcriptional modification, where adenosine (A) residues are deaminated to inosine (I) by ADAR enzymes. The editing process can alter the RNA structure and function, impacting gene expression and protein diversity20. In sequence alignments, A-to-I during editing is interpreted as a G by reverse transcription. RNA A-to-I editing is often represented as an A-to-G change on the edited RNA strand and as a T-to-C change on the corresponding DNA template strand. dRNA sequencing detects A-to-I editing directly from RNA data, as it omits the reverse transcription step, offering a precise view of RNA modifications without potential artifacts.

While no definitive editing site reference is available from GIAB, we generated a high-quality benchmark set of editing sites by cross-validating potential editing sites with DNA-seq data, as shown in Fig. 6a. The rationale is that an editing site is determined if present in RNA data but absent in DNA data. Initially, we identified A-to-G and T-to-C candidates from RNA alignments and computed read support in the corresponding DNA-seq data. If the sequencing coverage (DP) and variant allelic fraction (VAF) supporting either editing allele (G or C) in RNA (R) and DNA (D) is calculated as DPR, VAFR, DPD, and VAFD, respectively, then the editing site is decided as follows:

isEditing = True if DPR ≥ MDP and VAFR ≥ MAFR and DPD ≥ MDP and VAFD < MAFD
False else

where MDP represents minimal coverage depth, and VAFD and VAFR represent the VAF observed in DNA and RNA, respectively. MAFD represents maximal allelic fraction supporting A-to-C or G-to-C in DNA, and MAFR represents minimal allelic fraction supporting A-to-C or G-to-C in RNA. We set MDP, MAFD, and MAFR as 8, 0.08, and 0.75, respectively, for both the PacBio and ONT experiments. Supplementary Fig. 2 presents the detected allele frequency distribution of high-quality RNA editing sites across PacBio and ONT datasets.

We incorporated candidates from the GIAB high-confidence regions for model training. These candidates were categorized in RNA BAM files as RNA variants, RNA editing sites, or artifacts. Throughout our experiments, we explored the inclusion of a separate RNA editing site task to identify RNA edits. Additionally, we tried incorporating the editing sites as an additional category within the zygosity task. We found that both approaches failed to accurately identify true RNA variants compared to directly integrating editing sites into the homozygous reference category. This indicates that the model is able to classify the editing sites as references in model training.

RNA editing site tagging using the REDIportal database
REDIportal is a comprehensive database that consolidates a vast array of RNA editing sites sourced from three main repositories: ATLAS, DARNED, and RADAR3,30,31. Within REDIportal v2, approximately 40 million RNA editing sites are documented. In Clair3-RNA, variants are tagged as RNA editing sites if they correspond to A-to-G or T-to-C alterations found in the REDIportal database, ensuring a matching alternate representation. REDIportal prioritizes RNA editing sites validated by agreement of at least two of its constituent sources (ATLAS, DARNED, RADAR), with adjustable configuration setting for users’ customs. Clair3-RNA supports REDIportal versions customized for both the GRCh37 and GRCh38 reference genomes, providing a reliable resource for exploring established RNA editing events.

Incorporating phasing information into neural network
Haplotype assembly, or phasing, is a fundamental technique for reconstructing an individual’s haplotypes36,41. Long-read sequencing generates reads spanning tens of kilobases or more, enabling the capture of haplotype information often lost in short-read sequencing data. This method has been shown to significantly enhance variant calling performance in long-read sequencing data. In comparison to lrDNA-seq, which typically has read lengths of 10–100k for phasing, lrRNA-seq can still achieve read lengths of 1–3k, allowing for phasing in specific regions. However, to date, no methods have investigated the impact of phasing on the performance of lrRNA-seq. In this study, we attempted to incorporate additional phasing information into Clair3-RNA architecture to gain a deeper understanding of its contribution to lrRNA-seq performance.

Variant phasing and read haplotagging. As illustrated in Fig. 7a, during the training phase, we selected only the heterozygous germline SNP variants identified in both Clair3-RNA and the GIAB truths for variant phasing. This approach was chosen to mitigate the potential introduction of false positives, which could adversely affect phasing accuracy and lead to cumulative phasing errors during model training. Incorporating only GIAB truths might introduce variants that have insufficient read coverage and discrepancies in zygosity. Alignments are then haplotagged using the phased variants using the WhatsHap “haplotag” submodule, which assigns each phased read the “HP:x” tag, where x indicates the haplotype assigned.

During model training, read subsampling is adopted. We first applied variant phasing and read haplotagging to the BAM file with a full-coverage spectrum. The phased full-coverage reads were then subsampled into varying coverage levels for data augmentation. For candidates with phased reads, we randomly switched the phased reads to unphased reads by discarding the “HP” tag. This switching process simulates phasing scenarios for low-coverage BAM files, where the read haplotagging percentage is lower, while avoiding introducing incorrect phasing that could occur when applying phasing to each subsampled coverage. All subsampled candidates and switched candidates were generated with additional phasing features for model training. In the inference phase, the heterozygous variants were acquired from Clair3-RNA’s first round called variants. We enabled the “--distrust-genotypes” option in WhatsHap phasing to allow for switching variant zygosities in the event of incorrectly called zygosities.

Pileup phasing feature. We incorporated 12 pileup phasing features for each candidate position, represented as XH1 and XH2, where X∈{A, C, G, T, I, D}. The phasing features denote the total counts of nucleotides A, C, G, and T, as well as insertions (I) and deletions (D) for each candidate site in the phased paternal and maternal haplotypes, derived from the “HP” tag of the phased reads. The 12 additional pileup phasing features were concatenated with the original 18 pileup features, resulting in a total of 30 features that were input into the neural network for model training and inference.

Analysis of overlap and error statistics for GENCODE annotated gene
The gene list was obtained from the GENCODE database (version 46), focusing specifically on the “gene” category for analysis. We categorized the genes based on their overlap with callable high-confidence regions, defined as the intersection of the callable region with read support and the GIAB high-confidence regions. Giving each gene coordinate and all its overlapping high-confidence regions (hc), we calculate the intersection size using the formula:

LO(g) = sum_{i=1}^n max(0, min(gend, hcend_i) - max(gstart, hcstart_i))

Where n represents the number of overlapping regions. g and hc refer to gene and callable high-confidence region, respectively. The overlap ratio score H(g) for each gene is then computed as follows:

H(g) = LO(g) / L(g)

Where L(g) denotes the gene size, a gene is categorized as “wholly covered” if H(g) equals 1. Conversely, it is classified as “partially covered” if the total overlap length H(g) is greater than or equal to α, where α is a threshold percentage for partial overlapping and is manually set to 0.75. If H(g) equals 0, the gene is classified as “not covered”.

A gene is considered to have no SNP or Indel errors if true positives are present within the specified gene coordinates, with no false positives or false negatives in the same region, based on benchmarking results obtained from hap.py. A gene is considered to have no switch errors if all heterozygous variants in the given phasesets are phased correctly. Switch errors36 were assessed using the WhatsHap “compare” submodule within the given overlapping phasesets of lrDNA-seq and lrRNA-seq. The editing distance of the phased heterozygous variants in these two groups was computed as switch errors. The command used for this analysis is given in the Supplementary Notes – Command line used section.

Data preparation
PacBio mapping with pbmm2 and minimap2. PacBio Iso-Seq and MAS-Seq reads were aligned with pbmm2 or minimap2 with the configurations provided in the Supplementary Notes. ONT reads were aligned with minimap2 splice mode. As discussed in the “Results—Performance on PacBio” section, we found miniamp2 alignment output performed in RNA data alignment, even though pbmm2 was optimized for PacBio data.

ONT dRNA004 base-calling and mapping. The ONT dRNA002 and cDNA reads were base called using Guppy (version 6.4.6), and the ONT dRNA004 reads were base called using Dorado (version 3.2.4). Base called reads with low read quality scores (“qs”<10) were labeled as “FAILED” and discarded, as recommended by the ONT protocol. Subsequently, the base called reads were mapped to the GRCh38 human reference genome using minimap2 for evaluation.

ONT SQK-RNA004 kit library preparation and sequencing. The RNA-seq libraries were constructed following the ONT’s SQK-RNA004 protocol. RNA sourced from GIAB HG002, HG004, and HG005 samples was obtained from the Coriell Institute. The total RNA underwent assessment for concentration, purity, and integrity using Nanodrop and Qubit analyses. Subsequently, the RNA fragments were ligated with a sequencing adapter utilizing ONT’s latest SQK-RNA004 direct RNA-seq kit. Finally, the RNAs were sequenced on the PromethION flow cell using our in-house PromethION 2 Solo sequencing device, running MinKNOW software version 1.18.02 for 96 h.

Model training and evaluation. Model training was conducted using the known variants in the GIAB HG002 dataset, encompassing the full spectrum of read coverage, along with four subsampled coverage levels. For all pre-trained models described in the manuscript, we consistently used variants from chromosomes 1 to 19, 21, and 22 for training, with a random selection of 90% of the data for training and the remaining 10% for evaluation, while always holding out chromosome 20 in training. Performance evaluation was conducted in callable regions defined in the “Methods—Benchmarking metrics” section. Clair3-RNA trains the model for either 30 epochs or until the validation loss does not decrease for five consecutive epochs, employing an exponentially decaying learning rate until convergence. The epoch that performs best on the validation set will be selected for benchmarking and as the production model. The training process utilizes the Rectified Adam optimizer42 and includes an initial weight warm-up.

The variant and artifact identification in our study was significantly affected by the coverage and error rate distributions in the training datasets. As shown in Supplementary Fig. 1, artifact candidates substantially outnumbered SNPs and Indels during model training, and multi-allelic variants were far less frequent than homozygous and heterozygous variants. To address data and category imbalance issues, we used Focal loss43 during model training - an approach first introduced in our previous work23. Focal loss improves upon standard cross-entropy by incorporating a modulating factor that automatically down-weights well-classified examples while maintaining focus on challenging, underrepresented cases. This approach proved particularly effective for our application, as it ensured proper attention to rare but critical variant categories despite the overwhelming prevalence of artifacts in the training data, ultimately leading to more stable model convergence and improved performance.

Benchmarking metrics
RNA sequences represent transcribed segments of the genome, encompassing mainly the exon regions. Notably, the number of total variants in RNA was considerably lower than that in DNA, where there were ~27,000 to ~318,000 variants, depending on the sequencing platform, RNA throughput, and specific gene expression. In contrast, ~4 million true variants were available for DNA data. Figure 4a shows the benchmarking workflow for RNA variant evaluation. We utilized Clair3-RNA’s “get_rna_bed” submodule to identify callable regions with minimum read support (default set at 4, adjustable through options). These callable regions were then cross-referenced with the GIAB high-confidence regions to establish reliable benchmarking regions. Within the benchmark VCF from GIAB, variants with less than two reads supporting the alternate allele were excluded. We used hap.py44 to compare the variant-called variants against the truth variants in callable high-confidence regions. We used Clair3-RNA’s “calculate_overall_metrics” submodule to generate three metrics, “Precision”, “Recall”, and “F1-score”, for five categories: “Overall”, “SNP”, “Indel”, “Insertion”, and “Deletion”.

Computational performance
Clair3-RNA was implemented with Python and TensorFlow and leveraged PyPy for speedup. RNA variant calling using Clair3-RNA requires only CPU resources. Using one replicate sequencing data, Clair3-RNA takes ~30 min for ONT dRNA004 or PacBio MAS-Seq, using two 12-core Intel Xeon Silver 4116 processors. The memory consumption was capped at 500 MB per process. Model training requires a high-end graphics processing unit (GPU). We used TensorFlow to train models on an Nvidia GeForce RTX™ 3090 Ti and 4090 Ti GPU. We used the GNU time command to benchmark the runtime and memory usage.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Clair3-RNA is a deep learning based small variant caller for long read RNA sequence data that builds on the Clair3 long read caller by the authors and incorporates modifications specific to long-read RNA sequence data. The key positives of the method are that it achieves better accuracy that a recently released tool (longcallR) and DNA-seq variant callers on PacBio and ONT RNA-seq datasets. The authors also evaluate the impact of using haplotype phase information and RNA-editing databases on variant calling accuracy.\nMajor comments:\n1. The method itself builds on the Clair3 variant caller. The authors should clearly specify what components of the new method are the same as Clair3 and those that are novel. Variant calling using rna-seq is likely to have more errors near splice junctions, this point has been noted by previous papers. However, I did not see any mention or discussion on this point.\n2. Benchmarking\n(a) The authors compare their method to LongcallR (cited via a github link). This method is not published and there is no paper/preprint describing this method. The authors refer to the github repository (\"LongcallR employs a Bayesian statistical model for variant calling, demanding higher sequencing coverage for accuracy, as suggested by the authors\") but I am not sure how readers can find this information. Also, the results reported on the github page for this tool for the accuracy on PacBio data are much better than the authors report. Perhaps, some explanation would be useful.\n(b) The authors compare their method to DeepVariant and Clair3 (DNA-based methods). However, there is a RNA-based version of the DeepVariant caller published recently in Bioinformatics advances (ref #3). A recent paper showed that \"transformation of alignment files improves performance of variant callers for long-read RNA sequencing data ([https://genomebiology.biomedcentral.com/articles/10.1186/s13059-023-02923-y)\". This paper is not cited in the manuscript. The implication of this paper is that the DeepVariant-RNA model can easily be applied to long-read RNA-seq data by transforming the input BAM files and using the DeepVariant-RNA model. I think that this would be a fair comparison.] (https://genomebiology.biomedcentral.com/articles/10.1186/s13059-023-02923-y)\n3. The authors state \"However, there exists no deep learning-based caller designed specifically for variant calling in lrRNA-seq yet.\" Following up on the previous point, I don't see why the DeepVariant-RNA caller cannot be applied to long reads by applying the transformation reported in the Genome Biology paper. In addition, apart from the phasing information, it is not clear what difference is between the Clair3 model for variant calling on long-read and short-read RNA-seq model? Can the same model approach be applied to short read data?\n4. Numbers providing the fraction of the genome/transcriptome that is callable using RNA-seq reads (before and applying the coverage thresholds) should be provided in the results.\n5. The DeepVariant-RNA paper reported \"We find that DV RNA-seq [GTEx] GQ ≥ 18 is able to achieve a precision of 0.998 on SNPs and 0.989 on INDELs with similar levels of sensitivity we observed when deriving our cutoff (Table 3).\nMinor points:\nDiscussion first line: Variant calling in RNA-seq datasets offers more comprehensive insights into the transcriptome. Please provide references.\nLine 61 in introduction: References for the higher error rate in lr-RNA-seq compared to DNA-seq would be helpful. Also, what is the reason for the higher error rate?\nThe authors use the term \"zygosity switching\" (first in the introduction and later as well). This is not a standard term unlike RNA editing.\nPlease refer to the DeepVariant-RNA paper when describing the RNA-editing database analysis since a similar analysis was done.", "reviewer_2": "In this manuscript, the authors describe a computational method Clair3-RNA to call small variants, like SNV and short indels, from long-read RNA-seq data, which is also the first deep-learning-based method for this task. The authors have conducted very comprehensive evaluations and demonstrate that Clair3-RNA achieved substantially higher accuracy than other DNA-based approaches. In addition to the conventional analysis and techniques, the authors also evaluated the RNA editing events, and utilized phasing information to improve the method. The manuscript is well-written, where methods are clearly described and the evaluations are fair and comprehensive. Nevertheless, there are a few concerns the author may address to strengthen the manuscript.\n1. For each position, the deep learning model in Clair3-RNA takes a flanking region to form the pileup tensor. It seems that this flanking region is on the consecutive region on the DNA, so for the position near the exon boundary, the pileup may become empty and shallower for the intronic region. Similar patterns may occur in alternative splicing events. If this is the case, the authors may consider adding the evaluation of variant calling accuracy based on how far the position is away from the closest exon boundary.\n2. Since Clair3 for variant calling from DNA data is developed by the same lab, is it possible to retrain or fine-tune Clair3 with the RNA-seq data’s pileup, and then compare with Clair3-RNA with the adjusted Clair3? This could demonstrate the improvement from Clair3-RNA is due to suitable structures and other feature engineering. Simiarly, a recent paper: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-023-02923-y claims that by processing the alignment file properly, the DNA-based variant caller approach may achieve good accuracy for PacBio long-read IsoSeq data.\n3. When Clair3-RNA utilizes the REDIportal database for RNA editing site tagging, does the strand information take into account? For example, whether A-to-G or T-to-C alterations are both allowed, or one of them will be eliminated if the read alignments clearly indicate the pileup region corresponds to a specific strand RNA?\n4. Two minor concerns:\nFigure 8 is referred to between Figure 3 and Figure 4 in the text. The authors may consider reordering the figures or moving Figure 8 to the supplementary file, so the figure reference can be in order.\nIn the speed comparison, why would Clair3-RNA be so much faster than Clair3?", "reviewer_3": "The manuscript introduces Clair3-RNA, a new deep learning method for small variant detection from long-read RNA-seq data (lrRNA-seq). Clair3-RNA achieves strong performance in several ONT and PacBio benchmarks, substantially outperforming other benchmarked tools. The manuscript is very well written and includes numerous well-designed figures and tables. Notably, the authors performed extensive and informative experiments to examine the impact of key data attributes on performance (e.g., genomic context, coverage, zygosity switching), which were also thoughtfully incorporated into method design.\nMajor concerns:\n* The HG004 and HG005 benchmark results appear to be only reported for the entire genome in each sample. While the authors mention that they left out chromosome 20 during model training on the HG002 genome, I could not find the results specific only to chromosome 20 for HG004 and HG005. Since the HG002 sample (son) is related to the HG004 sample (mother), the model has been trained on variants which also appear in the test set; it is hence important to test performance specifically on variants unseen during training.\n* Additional metrics describing the generated training data are needed to evaluate the approach (i.e. how balanced/diverse was the training data), including: how many true variants, artifacts, and editing sites were included; what types of artifacts are expected in this data and were any artifact sub-categories specifically identified/selected; what was the distribution of allele frequencies in the training data.\n* In the paragraph starting with line 638, the authors mention that they explored including RNA edit site prediction as a separate task or category within the zygosity task but these approaches did not produce accurate results (why?). However, when integrating editing sites into the homozygous reference category, the model performed well (why?). This discrepancy in performance could point to a limitation of the input representation or the training dataset composition and should be further examined. How many edit sites versus reference sites versus true SNP sites were included during training and evaluation?\nGiven the features provided to the model as input, what distinguishing information is currently available to the model to differentiate between RNA edit sites and true variants?\n* The paper includes comparisons to only a single method designed specifically for lrRNA-seq variant calling (namely, LongcallR). A few other methods should be cited and benchmarked to enable a better performance assessment — e.g., CTAT-Mutations, isoLASER, LongcallR-nn, or the recent lrRNA-seq model of DeepVariant.\nMinor:\n* In addition to the reported parallel runtime performance, runtime should also be reported for each method in single-threaded mode (since each method may employ different parallel processing schemes)."}, "author_response": "Responses to\nComments of reviewers\nClair3-RNA: A deep learning-based small variant caller for long-read RNA sequencing data\nNCOMMS-24-75685\nReviewer: #1\n----------------------------------------------------------\nClair3-RNA is a deep learning based small variant caller for long read RNA sequence data\nthat builds on the Clair3 long read caller by the authors and incorporates modifications\nspecific to long-read RNA sequence data. The key positives of the method are that it\nachieves better accuracy that a recently released tool (longcallR) and DNA-seq variant\ncallers on PacBio and ONT RNA-seq datasets. The authors also evaluate the impact of\nusing haplotype phase information and RNA-editing databases on variant calling accuracy.\nMajor comments:\n1. The method itself builds on the Clair3 variant caller. The authors should clearly specify\nwhat components of the new method are the same as Clair3 and those that are novel.\nRE1: We thank the reviewer for their valuable comments for us to further improve our manuscript.\nWe have added the details in the Methods - Clair3-RNA model architecture section (line 573-584) to\nclarify the relationship between Clair3-RNA and Clair3.\n“Clair3-RNA retains the bidirectional long short-term memory (Bi-LSTM) model architecture and\noutput tasks from Clair3’s pileup model. The pileup network incorporates two Bi-LSTM layers and\nthree fully connected layers.” and “Clair3-RNA incorporates several successful deep-learning\ntechniques from Clair3 for model training, including the 1) model initialization; 2) activation\nfunction; 3) optimizer; 4) dropout; 5) L2 regularization; and 6) combination of multiple coverages for\nmodel training, and 7) focal loss function for addressing the unbalanced variant type problem.”\nEven using the same model architecture for model training, Clair3-RNA introduces several\nenhancements tailored to lrRNA-seq data, including new training sample generation strategies,\nadapted training labels for RNA-seq dataset, sample coverage normalization, incorporation of\nphasing information, etc, to improve the performance of lrRNA-seq data.\nVariant calling using rna-seq is likely to have more errors near splice junctions, this point has\nbeen noted by previous papers. However, I did not see any mention or discussion on this\npoint.\nRE2: We sincerely thank the reviewer for the valuable suggestion. We have added a dedicated section\nin the Results titled “Performance analysis at splice junctions” to address this point (line 477-489).\nTo analyze performance in splice junction regions, we focused on well-defined splice junctions from\nthe GENCODE dataset. Specifically, we extracted the boundaries of all exon regions and restricted\nour benchmark analysis to regions flanking splice junctions. The benchmarking was performed based\non the distance to the exon boundary, considering only those boundaries with read support on either\nthe left or right side.\nWe evaluated performance at distances of 2, 4, 6, 8, 10, 12, 14, and 16 base pairs from the splice\njunction boundary, to understand how the splice site would affect model performance. Notably, 16 bp\nrepresents the maximum scanning window size for neural network prediction.\nA decline in F1-score was observed across all datasets as variants approached the splice junction\nboundary (from ≤16 bp to ≤2 bp). Specifically, the average F1-score dropped from 90.31% to 83.98%\nfor PacBio Iso-Seq, from 92.52% to 85.45% for MAS-Seq, and from 91.45% to 85.30% for ONT\ndRNA004 using the best performing Clair3-RNA with phasing model. This highlights the inherent\nchallenges of variant calling in regions close to splice junctions.\nFigure R1. Performance of different callers based on the minimum distance to the nearest splice\njunction boundary in various datasets.\nWe observed that the quality score and performance tend to drop more significantly for variants\nflanking splice junction regions. We attribute this decline to two primary reasons: 1) The lack of\nsufficient variants near splice sites for a comprehensive model training. 2) The accumulation of\nsequencing errors in truncated exon regions leads to less confident model predictions.\nDP≥4 and AD≥2\nPlatform Dataset Sample Caller\n≤2 ≤4 ≤6 ≤8 ≤10 ≤12 ≤14 ≤16\nHG004 Clair3-RNA with phasing 82.31% 85.52% 87.39% 88.19% 89.05% 89.47% 90.13% 90.76%\nHG004 Clair3-RNA 81.51% 84.49% 85.81% 86.82% 87.81% 88.10% 88.87% 89.49%\nHG004 DeepVariant 63.32% 67.55% 68.38% 70.34% 72.36% 73.59% 74.89% 75.68%\nHG004 LongcallR 69.78% 74.28% 76.60% 78.40% 78.72% 79.02% 79.72% 80.03%\nIso-Seq HG004 Clair3 52.17% 60.87% 58.67% 61.05% 63.64% 67.52% 67.72% 69.41%\n(BCM) HG005 Clair3-RNA with phasing 85.45% 88.26% 90.14% 90.59% 91.17% 91.55% 91.87% 92.52%\nHG005 Clair3-RNA 83.31% 86.31% 87.56% 88.28% 89.28% 89.78% 90.40% 90.93%\nHG005 DeepVariant 64.31% 67.42% 69.07% 71.20% 73.37% 74.49% 75.51% 76.25%\nHG005 LongcallR 72.29% 74.09% 77.19% 78.94% 79.85% 80.22% 80.63% 81.10%\nHG005 Clair3 66.67% 65.45% 66.67% 63.92% 67.19% 69.94% 68.13% 67.92%\nHG004 Clair3-RNA with phasing 81.08% 82.29% 84.03% 85.13% 86.63% 87.33% 88.08% 88.83%\nHG004 Clair3-RNA 78.96% 80.55% 82.17% 83.23% 84.84% 85.64% 86.35% 87.13%\nPacBio HG004 DeepVariant 18.15% 16.11% 17.20% 17.45% 18.03% 17.79% 18.35% 18.65%\nHG004 LongcallR 66.51% 69.88% 72.75% 74.66% 75.37% 75.65% 76.44% 76.68%\nIso-Seq HG004 Clair3 22.22% 33.33% 40.00% 41.18% 47.83% 49.77% 52.92% 53.38%\n(Google) HG005 Clair3-RNA with phasing 80.76% 83.91% 86.32% 87.32% 88.10% 88.44% 89.12% 89.83%\nHG005 Clair3-RNA 79.09% 81.43% 83.29% 84.29% 85.57% 85.95% 86.69% 87.65%\nHG005 DeepVariant 19.35% 18.50% 18.27% 18.69% 18.74% 18.31% 18.84% 19.08%\nHG005 LongcallR 70.75% 72.78% 75.73% 77.66% 78.75% 79.06% 79.64% 80.04%\nHG005 Clair3 28.57% 36.62% 41.12% 39.42% 43.33% 49.33% 49.60% 50.00%\nHG004 Clair3-RNA with phasing 84.90% 86.28% 87.48% 88.07% 88.71% 89.00% 89.22% 89.74%\nHG004 Clair3-RNA 81.36% 83.02% 84.12% 84.81% 85.70% 86.24% 86.52% 87.14%\nMAS-Seq HG004 DeepVariant 61.98% 64.16% 65.46% 66.51% 68.34% 69.52% 70.49% 71.07%\nHG004 LongcallR 80.43% 82.19% 84.02% 84.57% 84.77% 84.84% 84.87% 85.07%\nHG004 Clair3 75.86% 74.19% 81.32% 79.28% 82.12% 82.54% 83.78% 85.02%\nHG004 Clair3-RNA with phasing 87.64% 90.20% 90.59% 90.80% 91.33% 91.60% 91.76% 91.89%\nHG004 Clair3-RNA 86.30% 88.68% 89.21% 89.04% 89.57% 89.83% 90.07% 90.26%\nHG004 DeepVariant 67.46% 69.43% 70.68% 71.13% 72.56% 73.12% 73.68% 74.32%\nHG004 LongcallR 59.69% 63.87% 66.12% 68.13% 69.01% 69.13% 69.82% 70.32%\nHG004 Clair3 59.26% 54.24% 58.43% 58.93% 63.09% 65.92% 67.62% 68.97%\nONT dRNA004\nHG005 Clair3-RNA with phasing 89.84% 90.84% 91.63% 92.25% 92.53% 92.66% 92.91% 93.07%\nHG005 Clair3-RNA 89.60% 90.66% 91.18% 91.26% 91.43% 91.31% 91.65% 91.94%\nHG005 DeepVariant 70.99% 71.15% 71.61% 72.59% 74.45% 75.08% 75.83% 76.38%\nHG005 LongcallR 62.54% 65.22% 68.18% 70.11% 71.34% 71.93% 72.56% 72.97%\nHG005 Clair3 48.00% 59.02% 58.43% 60.87% 64.43% 66.67% 66.35% 67.51%\nTable R1. Performance of different callers based on the minimum distance to the flanking exon\nboundary in various datasets.\nInspired by these observations, starting from v0.2.2, we introduced the following options to improve\nthe variant detection near the splice junctions:\n1) providing an “--enable_padding_in_splice_junction_regions” option. This option completes the\nmissing features in the input tensor caused by hitting the exon boundaries. The missing features\nreduced the capability of the neural network of inferencing the correct variant. We filled the missing\nfeatures assuming the reference allele at each position with missing feature and a read coverage\ncalculated as the average of the flanking covered positions. Below are two randomly picked cases\nshowing how the padding was done.\nPacBio MAS-Seq HG004 chr3:138570305-138570336 PacBio MAS-Seq HG004 chr1:236012023-236012054\nIntronicregion Candidateposition Candidateposition Intronicregion\nPileupfeaturesbeforeenablingpadding Pileupfeaturesbeforeenablingpadding\n(reversestrand) (reversestrand)\nPileupfeaturesafterenablingpadding Pileupfeaturesafterenablingpadding\n(reversestrand) (reversestrand)\nFigure R2. Visualization of pileup features in Clair3-RNA with and without enabling the\n‘enable_padding_in_splice_junction_regions’ option for two cases.\n2) providing an “--enable_variant_calling_at_sequence_head_and_tail” option. This option forces\nvariant calling at sequence heads or tails where the reference allele is ambiguous, or there is no read\ncoverage. It is designed to still call variants in these regions even without sufficient information\nsupport.\nThe results presented in the figure below (Figure R3) demonstrate that enabling this option enhances\nperformance on splice junction regions across various datasets. The average SNP performance near\nthe exon boundary has improved by 1.01% for Iso-Seq and by 2.56% for MAS-Seq, respectively.\nVGI\nerutaef\npueliP\nFigure R3. The performance comparison of Clair3-RNA with or without the enabling the\n“enable_padding_in_splice_junction_regions” option in various datasets.\n2. Benchmarking\n(a) The authors compare their method to LongcallR (cited via a github link). This method is\nnot published and there is no paper/preprint describing this method. The authors refer to the\ngithub repository (\"LongcallR employs a Bayesian statistical model for variant calling,\ndemanding higher sequencing coverage for accuracy, as suggested by the authors\") but I\nam not sure how readers can find this information. Also, the results reported on the github\npage for this tool for the accuracy on PacBio data are much better than the authors report.\nPerhaps, some explanation would be useful.\nRE3: We thank the reviewer for the comment. The algorithmic details of LongcallR were obtained\nthrough its presentation at the RECOMB-SEQ conference, where it was accepted as a talk, and\nthrough personal communication with the authors. For clarity, we have updated the manuscript’s\ndescription of LongcallR (line 170-171).\nThe results of LongcallR reported in our manuscript are consistent with the authors' findings. As the\nauthors stated in their GitHub repository, LongcallR requires that “All truth variants used for\nbenchmarking have at least 10X coverage in the input.”. This indicates that they only include variants\nwith DP≥10 in their benchmarking results in the GitHub repository. In contrast, our benchmarking\nevaluates all variants with read coverage ”DP≥4” by default.\nIn the authors’ GitHub README (https://github.com/huangnengCSU/longcallR?tab=readme-ov-\nfile#results), LongcallR achieved an F1-score of 94.38% on the PacBio MAS-Seq HG004 dataset. In\nour Supplementary Table 3, we reported a comparable performance of 93.18% F1-score for variants\nwith DP ≥ 10 and disregarding zygosity. We believe the authors’ decision to exclude variants with DP\n< 10 is primarily because it is challenging for the statistical model to handle complex cases in low-\ncoverage scenarios.\n(b) The authors compare their method to DeepVariant and Clair3 (DNA-based methods).\nHowever, there is a RNA-based version of the DeepVariant caller published recently in\nBioinformatics advances (ref #3).\nRE4: Until our first preprint release on bioRxiv (November 2024), DeepVariant had only provided a\npre-trained model for Illumina RNA-seq data, as described in Bioinformatics Advances. In December\n2024, DeepVariant released v1.8.0 to support variant calling model for PacBio Iso-Seq/MAS-Seq\n(https://github.com/google/deepvariant/releases/tag/v1.8.0). Hence, we were not able to benchmark\nthe DeepVariant PacBio RNA-seq model in the previous submission.\nAs shown in the table (Table R2) below, we assessed the performance of the DeepVariant Iso-Seq and\nMAS-Seq model (referred to as the “DeepVariant (RNA model)” in the table below) across the\ncorresponding datasets. The DeepVariant RNA model achieved SNP F1-scores of 91.47%, 92.77%,\nand 91.92% on the Iso-Seq HG004, Iso-Seq HG005, and MAS-Seq HG004 datasets, respectively.\nNotably, the DeepVariant RNA model significantly outperformed its DNA model, indicating the\nimportance of model training tailored for lrRNA-Seq. However, these results were ~3.5% lower in the\nF1-score compared to the best-performing Clair3-RNA model with phasing enabled.\nSNP (DP≥4 and AD≥2) Indel (DP≥4 and AD≥2)\nPlatform Dataset Sample Caller\nPrecision Recall F1-score Precision Recall F1-score\nHG004 Clair3-RNA with phasing 96.81% 94.76% 95.77% 83.49% 76.02% 79.58%\nHG004 Clair3-RNA 96.16% 93.73% 94.93% 80.74% 70.21% 75.11%\nHG004 DeepVariant (RNA model) 89.87% 93.12% 91.47% 78.89% 67.53% 72.77%\nIso-Seq HG004 DeepVariant (DNA model) 76.03% 43.18% 55.08% 55.63% 39.68% 46.32%\n(BCM) HG005 Clair3-RNA with phasing 97.13% 95.81% 96.47% 88.40% 82.99% 85.61%\nHG005 Clair3-RNA 96.48% 94.53% 95.50% 85.47% 77.69% 81.40%\nPacBio\nHG005 DeepVariant (RNA model) 91.79% 93.78% 92.77% 83.33% 76.47% 79.75%\nHG005 DeepVariant (DNA model) 72.74% 39.83% 51.48% 63.80% 39.08% 48.47%\nHG004 Clair3-RNA with phasing 95.81% 93.44% 94.61% 79.07% 67.28% 72.70%\nHG004 Clair3-RNA 95.23% 91.61% 93.38% 76.48% 60.71% 67.69%\nMAS-Seq\nHG004 DeepVariant (RNA model) 94.29% 89.68% 91.92% 83.28% 60.34% 69.98%\nHG004 DeepVariant (DNA model) 72.04% 28.26% 40.59% 43.15% 23.22% 30.19%\nTable R2. Performance comparison of DeepVariant RNA model and DNA model.\nA recent paper showed that \"transformation of alignment files improves performance of\nvariant callers for long-read RNA sequencing data\n([https://genomebiology.biomedcentral.com/articles/10.1186/s13059-023-02923-y)\". This\npaper is not cited in the manuscript. The implication of this paper is that the DeepVariant-\nRNA model can easily be applied to long-read RNA-seq data by transforming the input BAM\nfiles and using the DeepVariant-RNA model. I think that this would be a fair comparison.]\n(https://genomebiology.biomedcentral.com/articles/10.1186/s13059-023-02923-y)\nRE5: We thank the reviewer for the suggestion for benchmarking the method. The paper describes a\ntransformation and manipulation of RNA-seq alignments to adapt DNA variant callers for RNA-seq\ndata. Following the authors’ practice, we applied this transformation to the alignment files. The\nworkflow for the alignment transformation is detailed in the Supplementary Notes - Command line\nused Section. Then we ran Clair3 and DeepVariant on the transformed BAM files, as recommended\nby the authors. We followed the authors’ approach\n(https://github.com/vladimirsouza/lrRNAseqVariantCalling) to combine Clair3 pileup SNP and Indel\nVCFs to generate the Clair3-mix results.\nThe results are reported in the figure (Figure R4) and table (Table R3) below. We observed that the\nperformance of both Clair3 (Clair3-mix) and DeepVariant (DeepVariant -mix) improved when using\nthe transformed alignments compared to the original spliced alignments. This demonstrates that the\ntransformation workflow can enhance the performance of DNA germline variant callers when applied\nto RNA-seq data. However, Clair3-RNA still outperformed both Clair3-mix and DeepVariant-mix,\nunderscoring the advantage of designing an RNA-seq-specific caller over simply adapting DNA\ncallers to transformed BAM files. We have also added a description of the BAM transformation\nmethod in the Introduction section.\nFigure R4. Performance comparison between native lrRNA-seq variant callers and lrDNA-seq callers\nrepurposed for RNA-seq data.\nWith DP≥4 and AD≥2\nPlatform Dataset Sample alignment Caller SNP Indel\ntransformation Precision Recall F1-score Precision Recall F1-score\nHG004 √ DeepVariant-mix 83.47% 73.36% 78.09% 69.69% 50.46% 58.53%\nHG004 DeepVariant 76.03% 43.18% 55.08% 55.63% 39.68% 46.32%\nHG004 √ Clair3-mix 54.99% 91.24% 68.62% 49.53% 61.13% 54.72%\nHG004 Clair3 pileup 54.07% 93.08% 68.40% 45.92% 62.83% 53.06%\nHG004 Clair3-RNA 96.16% 93.73% 94.93% 80.74% 70.21% 75.11%\nIso-Seq HG004 Clair3-RNA with phasing 96.81% 94.76% 95.77% 83.49% 76.02% 79.58%\n(BCM) HG005 √ DeepVariant-mix 81.73% 76.74% 79.16% 76.09% 62.21% 68.45%\nHG005 DeepVariant 72.74% 39.83% 51.48% 63.80% 39.08% 48.47%\nHG005 √ Clair3-mix 53.23% 92.54% 67.58% 58.17% 72.14% 64.41%\nHG005 Clair3 pileup 52.93% 94.72% 67.91% 54.20% 74.30% 62.68%\nHG005 Clair3-RNA 96.48% 94.53% 95.50% 85.47% 77.69% 81.40%\nHG005 Clair3-RNA with phasing 97.13% 95.81% 96.47% 88.40% 82.99% 85.61%\nHG005 √ DeepVariant-mix 83.82% 44.81% 58.40% 65.48% 17.43% 27.53%\nHG005 DeepVariant 73.91% 12.71% 21.69% 41.33% 13.94% 20.85%\nHG005 √ Clair3-mix 44.28% 88.68% 59.07% 28.76% 54.46% 37.64%\nPacBio\nHG005 Clair3 pileup 41.42% 90.16% 56.77% 23.75% 55.74% 33.31%\nHG004 Clair3-RNA 94.09% 89.69% 91.84% 63.08% 65.02% 64.03%\nIso-Seq HG004 Clair3-RNA with phasing 95.77% 91.53% 93.60% 73.93% 70.61% 72.23%\n(Google) HG004 √ DeepVariant-mix 84.24% 47.46% 60.72% 70.70% 21.51% 32.98%\nHG004 DeepVariant 71.53% 11.79% 20.25% 44.52% 13.83% 21.11%\nHG004 √ Clair3-mix 43.17% 89.71% 58.29% 29.27% 63.56% 40.08%\nHG004 Clair3 pileup 40.74% 91.46% 56.38% 23.94% 65.12% 35.01%\nHG005 Clair3-RNA 94.27% 90.09% 92.13% 61.90% 70.49% 65.92%\nHG005 Clair3-RNA with phasing 95.88% 92.07% 93.93% 74.80% 75.58% 75.19%\nHG004 √ DeepVariant-mix 75.70% 57.66% 65.46% 60.40% 36.98% 45.87%\nHG004 DeepVariant 72.04% 28.26% 40.59% 43.15% 23.22% 30.19%\nMAS-Seq HG004 √ Clair3-mix 50.76% 91.21% 65.22% 36.14% 60.92% 45.37%\nHG004 Clair3 pileup 63.96% 87.66% 73.95% 40.67% 45.34% 42.88%\nHG004 Clair3-RNA 95.23% 91.61% 93.38% 76.48% 60.71% 67.69%\nHG004 Clair3-RNA with phasing 95.81% 93.44% 94.61% 79.07% 67.28% 72.70%\nHG004 √ DeepVariant-mix 88.07% 52.79% 66.01% 49.81% 27.59% 35.51%\nHG004 DeepVariant 86.53% 43.55% 57.94% 46.60% 26.79% 34.02%\nHG004 √ Clair3-mix 51.20% 80.13% 62.48% 18.07% 47.79% 26.23%\nHG004 Clair3 pileup 42.11% 85.10% 56.34% 27.27% 43.95% 33.65%\nHG004 Clair3-RNA 92.82% 89.26% 91.00% 83.61% 47.19% 60.33%\nHG004 Clair3-RNA with phasing 94.04% 92.35% 93.19% 87.36% 52.50% 65.59%\nONT\ndRNA004 HG005 √ DeepVariant-mix 88.70% 54.18% 67.27% 57.87% 32.80% 41.86%\nHG005 DeepVariant 87.14% 42.31% 56.96% 54.85% 29.04% 37.98%\nHG005 √ Clair3-mix 52.41% 83.75% 64.47% 20.13% 57.29% 29.79%\nHG005 Clair3 pileup 43.27% 87.72% 57.95% 33.30% 52.03% 40.61%\nHG005 Clair3-RNA 93.16% 90.34% 91.73% 85.64% 53.76% 66.05%\nHG005 Clair3-RNA with phasing 94.23% 92.90% 93.56% 88.55% 58.74% 70.63%\nTable R3. Performance comparison between native lrRNA-seq variant callers and lrDNA-seq callers\nrepurposed for RNA-seq data.\n3. The authors state \"However, there exists no deep learning-based caller designed\nspecifically for variant calling in lrRNA-seq yet.\" Following up on the previous point, I don't\nsee why the DeepVariant-RNA caller cannot be applied to long reads by applying the\ntransformation reported in the Genome Biology paper.\nRE6: We thank the reviewer for the comment. As mentioned in our response to RE5, the\ntransformation method is specifically designed for DNA-seq variant callers. Our benchmarking\nresults demonstrate that callers trained directly on lrRNA-seq datasets, such as DeepVariant Iso-\nSeq/MAS-Seq model and Clair3-RNA, outperform Clair3-mix and DeepVariant-mix when applied to\nlrRNA-seq data. This highlights the importance of designing lrRNA-seq-specific callers, as they are\nbetter suited to learn and adapt to the unique characteristics of RNA-seq data, such as RNA editing,\nRNA uneven coverage, zygosity flipping events, etc. Regarding the performance of DeepVariant\nRNA model, Clair3-RNA outperform DeepVariant-RNA in our benchmarks, the results are shown in the\ntable (Table R2) above. The transformation results have been added to the revision of our manuscript\n(line 238-246).\nIn addition, apart from the phasing information, it is not clear what difference is between the\nClair3 model for variant calling on long-read and short-read RNA-seq model? Can the\nsame model approach be applied to short read data?\nRE7: Clair3-RNA can support variant calling on short-read RNA-seq data if the model is trained on\nIllumina RNA-seq datasets. As short-reads have lower error rate and a different error profile\ncompared to long-reads, we anticipate some adaptations to optimize the performance for short-read\nanalysis. That includes tuning the minimum read mapping quality and base quality, minimum SNP\nand Indel VAF cutoff for training, and changing the neural network architectures, etc. Nevertheless,\ncalling germline variant using short-read is not a primary target of Clair3-RNA as we believe\nDeepVariant, GATK and Platypus have already achieved state-of-the-art performance for short-read\nRNA-seq germline variant calling.\nAs highlighted in DeepVariant RNA model’s paper (Cook et al., 2023), DeepVariant made use of the\nGenotype-Tissue Expression (GTEx) consortium short-read dataset for model training. This approach\ngenerated over 80 million training samples, which is ~300 times larger than that could be generated\nfrom the GIAB RNA-seq dataset, only in which lrRNA-seq data is available. Consequently,\nDeepVariant has achieved superior performance compared to models trained solely on the GIAB\nRNA-seq dataset. Following this proven practice, to further improve the performance of Clair3-RNA,\nwe will make use of any lrRNA-seq data made available by GTEx or other sources as soon as\navailable.\n4. Numbers providing the fraction of the genome/transcriptome that is callable utilizing RNA-\nseq reads (before and applying the coverage thresholds) should be provided in the results.\nRE8: We thank the reviewer for the valuable suggestion. We have calculated the genomic size for\ndifferent coverage settings and determined the fraction by dividing it by the reference genome size\nand the callable genomic regions.\nTo enhance clarity, we have updated the Results section with a description of these fractions (line\n250-251), helping readers better understand the coverage fractions. Additionally, the calculated\nfractions for all datasets are now included in Supplementary Table 7.\nPlatform Dataset Sample B ( a D se P ≥ p 1 a ) ir B ( a D se P ≥ p 4 a ) ir B ( a D se P ≥ p 8 a ) ir B (D as P e ≥ p 10 a ) ir (DP F ≥ a 4 c t / i D on P≥1)(DP F ≥ a 8 c t / i D on P≥1)(DP F ≥ a 1 c 0 ti o / D n P≥ (DP F ≥ a 1 c t / i w on hole (DP F ≥ a 4 c / t i w o h n ole (DP F ≥ a 8 c / t i w o h n ole(DP F ≥ a 10 c t / i o w n hole\n1) genome) genome) genome) genome)\nPacBio MAS-Seq HG004 732,094,285 454,991,885 388,341,698 368,586,367 62.15% 53.05% 50.35% 23.62% 14.68% 12.53% 11.89%\nPacBio MAS-Seq HG004 669,710,416 419,550,797 359,026,328 340,897,860 62.65% 53.61% 50.90% 21.60% 13.53% 11.58% 11.00%\nPacBio Iso-Seq (BCM) HG004 522,085,768 292,774,392 237,459,780 221,917,744 56.08% 45.48% 42.51% 16.84% 9.44% 7.66% 7.16%\nPacBio Iso-Seq (BCM) HG004 479,922,631 272,197,678 220,926,777 206,537,452 56.72% 46.03% 43.04% 15.48% 8.78% 7.13% 6.66%\nPacBio Iso-Seq (Google)HG004 416,457,292 212,747,323 176,321,158 166,446,157 51.09% 42.34% 39.97% 13.43% 6.86% 5.69% 5.37%\nPacBio Iso-Seq (Google)HG004 381,186,226 196,363,186 163,087,266 154,026,384 51.51% 42.78% 40.41% 12.30% 6.33% 5.26% 4.97%\nONT cDNA HG004 258,979,711 124,764,725 97,155,184 91,218,886 48.18% 37.51% 35.22% 8.35% 4.02% 3.13% 2.94%\nONT cDNA HG004 236,371,357 114,765,848 89,761,918 84,397,647 48.55% 37.97% 35.71% 7.63% 3.70% 2.90% 2.72%\nONT dRNA002 HG004 77,767,787 33,814,050 27,334,957 25,130,599 43.48% 35.15% 32.31% 2.51% 1.09% 0.88% 0.81%\nONT dRNA002 HG004 71,426,625 31,434,120 25,456,599 23,389,201 44.01% 35.64% 32.75% 2.30% 1.01% 0.82% 0.75%\nONT dRNA004 HG004 277,695,719 164,985,976 151,621,773 147,957,751 59.41% 54.60% 53.28% 8.96% 5.32% 4.89% 4.77%\nONT dRNA004 HG004 252,333,764 151,418,228 139,583,980 136,253,236 60.01% 55.32% 54.00% 8.14% 4.88% 4.50% 4.40%\nPacBio Iso-Seq (BCM) HG005 591,735,125 382,613,111 319,096,966 299,511,375 64.66% 53.93% 50.62% 19.09% 12.34% 10.29% 9.66%\nPacBio Iso-Seq (BCM) HG005 536,956,689 350,527,601 292,957,475 275,195,739 65.28% 54.56% 51.25% 17.32% 11.31% 9.45% 8.88%\nPacBio Iso-Seq (Google)HG005 408,821,925 216,989,506 178,078,942 167,110,256 53.08% 43.56% 40.88% 13.19% 7.00% 5.74% 5.39%\nPacBio Iso-Seq (Google)HG005 369,280,570 197,588,859 162,570,057 152,696,590 53.51% 44.02% 41.35% 11.91% 6.37% 5.24% 4.93%\nONT cDNA HG005 360,634,723 165,599,273 138,281,320 131,708,392 45.92% 38.34% 36.52% 11.63% 5.34% 4.46% 4.25%\nONT cDNA HG005 324,483,020 150,245,811 125,821,602 119,935,744 46.30% 38.78% 36.96% 10.47% 4.85% 4.06% 3.87%\nONT dRNA002 HG005 71,903,807 33,734,052 27,326,854 25,189,117 46.92% 38.00% 35.03% 2.32% 1.09% 0.88% 0.81%\nONT dRNA002 HG005 65,517,257 31,137,174 25,272,678 23,311,706 47.53% 38.57% 35.58% 2.11% 1.00% 0.82% 0.75%\nONT dRNA004 HG005 318,678,810 183,462,931 167,401,623 162,881,218 57.57% 52.53% 51.11% 10.28% 5.92% 5.40% 5.25%\nONT dRNA004 HG005 286,211,398 166,344,478 152,176,257 148,180,918 58.12% 53.17% 51.77% 9.23% 5.37% 4.91% 4.78%\nTable R4. The size and percentage of genomic regions at certain read coverages (DP) or above.\n5. The DeepVariant-RNA paper reported \"We find that DV RNA-seq [GTEx] GQ ≥ 18 is able\nto achieve a precision of 0.998 on SNPs and 0.989 on INDELs with similar levels of\nsensitivity we observed when deriving our cutoff (Table 3).\" In the text of the paper, authors\nprimarily used the F1-score to report the superior accuracy of their method. For RNA-seq\ndata, it would be useful to also provide some numbers on the recall that can be achieved at\na high precision value (similar to that provided by the DeepVariant-RNA paper).\nRE9: We thank the reviewer for the suggestion. We repeated the experiment in DeepVariant-RNA’s\npaper (regarding Table 3), when restricting to GQ ≥ 18, precision indeed improved from 0.978 to\n0.998, but the recall dropped significantly from 0.942 to 0.809, resulting in a decline of F1-score from\n0.960 to 0.893.\nTo illustrate the dynamics between precision and recall better in our study, we plotted the precision-\nrecall curve in Figure 3d, where each point represents the performance at a specific QUAL/GQ\ncutoff. We used the Area Under the Precision-Recall Curve (AUPRC) to summarize the performance\nof each caller across various cutoffs. Additionally, the table (Table R5) below showcases Clair3-\nRNA’s performance at multiple precision cutoffs (from 96% to 99%) and the corresponding GQ\ncutoffs for reference.\nDP≥4 and AD≥2 DP≥10 and AD≥4\nPrecision QUAL/GQ QUAL/GQ\nCaller Sample Dataset SNP SNP\ncutoff cutoff cutoff\nPrecision Recall F1-score Precision Recall F1-score\n96% 5.52 96.00% 93.42% 94.70% 5.52 96.00% 93.42% 94.70%\nPacBio\n97% 10.39 97.00% 87.62% 92.07% 10.39 97.00% 87.62% 92.07%\nIso-Seq\n98% 13.25 98.00% 74.45% 84.61% 13.25 98.00% 74.45% 84.61%\n(BCM)\n99% 16.8 99.00% 48.38% 65.00% 16.8 99.00% 48.38% 65.00%\n96% 7.68 96.00% 87.90% 91.77% 7.68 96.00% 87.90% 91.77%\nPacBio 97% 11.26 97.00% 77.59% 86.22% 11.26 97.00% 77.59% 86.22%\nHG004\nMAS-Seq 98% 12.69 98.00% 68.06% 80.33% 12.69 98.00% 68.06% 80.33%\n99% 16.12 99.00% 47.33% 64.04% 16.12 99.00% 47.33% 64.04%\n96% 12.62 96.00% 75.55% 84.56% 12.62 96.00% 75.55% 84.56%\nONT 97% 14.1 97.00% 66.63% 79.00% 14.1 97.00% 66.63% 79.00%\nClair3-RNA\ndRNA004 98% 15.83 98.00% 53.35% 69.09% 15.83 98.00% 53.35% 69.09%\n99% 17.84 99.00% 34.05% 50.67% 17.84 99.00% 34.05% 50.67%\n96% 4.79 96.20% 94.47% 95.33% 4.79 96.20% 94.47% 95.33%\nPacBio\n97% 7.72 97.00% 92.52% 94.71% 7.72 97.00% 92.52% 94.71%\nIso-Seq\n98% 12.46 98.00% 82.06% 89.32% 12.46 98.00% 82.06% 89.32%\n(BCM)\n99% 15.33 99.00% 63.14% 77.11% 15.33 99.00% 63.14% 77.11%\nHG005\n96% 12.25 96.00% 78.40% 86.31% 12.25 96.00% 78.40% 86.31%\nONT 97% 13.97 97.00% 67.92% 79.90% 13.97 97.00% 67.92% 79.90%\ndRNA004 98% 15.86 98.00% 52.82% 68.64% 15.86 98.00% 52.82% 68.64%\n99% 18.72 99.00% 24.37% 39.11% 18.72 99.00% 24.37% 39.11%\nTable R5. Impact of different precision cutoffs on the GQ cutoffs and performance of Clair3-RNA in\nvarious datasets.\nMinor points:\nDiscussion first line: Variant calling in RNA-seq datasets offers more comprehensive insights\ninto the transcriptome. Please provide references.\nRE10: We have added the reference accordingly.\nLine 61 in introduction: References for the higher error rate in lr-RNA-seq compared to DNA-\nseq would be helpful. Also, what is the reason for the higher error rate?\nRE11: We have added the reference accordingly. For long-read RNA sequencing, the raw signals\nproduced by nanopore and PacBio RNA sequencers are inherently more error-prone due to several\nfactors: 1) The sequencing chemistry, kits, and protocol settings for lrRNA-seq are less matured\ncompared to lrDNA-seq, which significantly impacts the accuracy of the raw signals. 2) Both PacBio\nand ONT basecallers are deep-learning-based methods that require extensive training materials.\nHowever, currently, the sample scale for long-read direct RNA sequencing is smaller compared to\nDNA sequencing, as RNA-seq products such as PacBio MAS-Seq and ONT dRNA004 kit were\ndeveloped more recently. 3) RNA sequences are more prone to degradation and are less stable than\nDNA, leading to higher sequencing error rates in RNA sequencing.\nThe authors use the term \"zygosity switching\" (first in the introduction and later as well). This\nis not a standard term unlike RNA editing.\nRE12: Thank the reviewer for the suggestion. We have renamed all “zygosity switching” to “zygosity\nflipping” and included a brief description (an observation in variant calling due to the expression of\nonly one haplotype) to the term to help readers better understand the meaning of the term.\nPlease refer to the DeepVariant-RNA paper when describing the RNA-editing database\nanalysis since a similar analysis was done.\nRE13: We have updated the manuscript accordingly.\nReviewer: #2\n----------------------------------------------------------\nIn this manuscript, the authors describe a computational method Clair3-RNA to call small\nvariants, like SNV and short indels, from long-read RNA-seq data, which is also the first\ndeep-learning-based method for this task. The authors have conducted very comprehensive\nevaluations and demonstrate that Clair3-RNA achieved substantially higher accuracy than\nother DNA-based approaches. In addition to the conventional analysis and techniques, the\nauthors also evaluated the RNA editing events, and utilized phasing information to improve\nthe method. The manuscript is well-written, where methods are clearly described and the\nevaluations are fair and comprehensive. Nevertheless, there are a few concerns the author\nmay address to strengthen the manuscript.\n1. For each position, the deep learning model in Clair3-RNA takes a flanking region to form\nthe pileup tensor. It seems that this flanking region is on the consecutive region on the DNA,\nso for the position near the exon boundary, the pileup may become empty and shallower for\nthe intronic region. Similar patterns may occur in alternative splicing events. If this is the\ncase, the authors may consider adding the evaluation of variant calling accuracy based on\nhow far the position is away from the closest exon boundary.\nRE14: We thank the reviewer for the valuable comments. Another reviewer had the same concern,\nand our new experiment show the variant calling performance of various distances to the exon\nboundary is shown in RE2. We found that the performance decreased gradually when a variant is\napproaching the exon boundaries. We also discussed how a new option called “--\nenable_padding_in_splice_junction_regions\" inspired by reviewers’ comments has further improved\nthe variant calling performance near exon boundary.\n2. Since Clair3 for variant calling from DNA data is developed by the same lab, is it possible\nto retrain or fine-tune Clair3 with the RNA-seq data’s pileup, and then compare with Clair3-\nRNA with the adjusted Clair3? This could demonstrate the improvement from Clair3-RNA is\ndue to suitable structures and other feature engineering.\nRE15: We thank the reviewer for the comment. We have added to the revision new results about fine-\ntuning Clair3 using RNA-seq data, following the training instructions and procedures outlined for\nDNA-seq in the Clair3 documentation (https://github.com/HKU-\nBAL/Clair3/blob/main/docs/pileup_training.md). For the model fine-tuning, we used the pre-trained\nDNA-seq model with a smaller learning rate (default: 1e-4) to train the neural network. We held out\n10% of the training data as a validation set and selected the epoch with the best performance in the\nvalidation dataset for testing.\nFor the PacBio Iso-Seq and MAS-Seq datasets, we used the default PacBio DNA model\n(https://www.bio8.cs.hku.hk/clair3/clair3_models/hifi.tar.gz). For the ONT dRNA004 dataset, we used\nthe latest Q20+ R10.4.1 chemistry pre-trained model “r1041_e82_400bps_sup_v410” from ONT\nRerio (https://github.com/nanoporetech/rerio?tab=readme-ov-file#clair3-models)\nThe results in the figure (Figure R5) below demonstrate that while the fine-tuned Clair3 model\nperformed better than the Clair3 DNA model, Clair3-RNA still performed the best, surpassing the\nfine-tuned model by 11.69%, 7.89%, and 13.60% in average F1-score in the PacBio Iso-Seq, MAS-\nSeq, and ONT dRNA004 datasets, respectively.\nFigure R5. Performance comparison of Clair3 DNA model and Clair3 fine-tuned model in various\ndatasets.\nDP≥4 and AD≥2\nDataset Sample Caller SNP Indel\nPrecision Recall F1-score Precision Recall F1-score\nClair3 fine-tuned from RNA-seq dataset 88.96% 82.40% 85.55% 67.67% 48.27% 56.35%\nClair3 DNA-seq model 68.73% 96.22% 80.18% 57.17% 69.35% 62.67%\nHG004\nClair3-RNA 96.16% 93.73% 94.93% 80.74% 70.21% 75.11%\nPacBio\nClair3-RNA with phasing enabled 96.81% 94.76% 95.77% 83.49% 76.02% 79.58%\nIso-Seq\nClair3 fine-tuned from RNA-seq dataset 89.73% 85.13% 87.37% 75.63% 59.90% 66.85%\n(BCM)\nClair3 DNA-seq model 66.80% 97.15% 79.16% 64.31% 76.18% 69.75%\nHG005\nClair3-RNA 96.48% 94.53% 95.50% 85.47% 77.69% 81.40%\nClair3-RNA with phasing enabled 97.13% 95.81% 96.47% 88.40% 82.99% 85.61%\nClair3 fine-tuned from RNA-seq dataset 84.92% 73.75% 78.94% 54.95% 38.65% 45.38%\nClair3 DNA-seq model 54.65% 94.23% 69.18% 35.56% 67.84% 46.66%\nHG004\nClair3-RNA 94.09% 89.69% 91.84% 63.08% 65.02% 64.03%\nPacBio\nClair3-RNA with phasing enabled 95.77% 91.53% 93.60% 73.93% 70.61% 72.23%\nIso-Seq\nClair3 fine-tuned from RNA-seq dataset 73.90% 89.09% 80.79% 39.07% 46.34% 42.40%\n(Google)\nClair3 DNA-seq model 70.33% 85.58% 77.21% 50.96% 49.50% 50.22%\nHG005\nClair3-RNA 93.16% 90.34% 91.73% 85.64% 53.76% 66.05%\nClair3-RNA with phasing enabled 94.23% 92.90% 93.56% 88.55% 58.74% 70.63%\nClair3 fine-tuned from RNA-seq dataset 91.89% 82.10% 86.72% 68.61% 47.79% 56.34%\nPacBio Clair3 DNA-seq model 75.38% 92.77% 83.17% 74.24% 45.07% 56.09%\nHG004\nMAS-Seq Clair3-RNA 95.23% 91.61% 93.38% 76.48% 60.71% 67.69%\nClair3-RNA with phasing enabled 95.81% 93.44% 94.61% 79.07% 67.28% 72.70%\nClair3 fine-tuned from RNA-seq dataset 73.10% 87.27% 79.56% 33.04% 38.66% 35.63%\nClair3 DNA-seq model 66.93% 84.89% 74.85% 41.23% 45.32% 43.18%\nHG004\nClair3-RNA 92.82% 89.26% 91.00% 83.61% 47.19% 60.33%\nONT Clair3-RNA with phasing enabled 94.04% 92.35% 93.19% 87.36% 52.50% 65.59%\ndRNA004 Clair3 fine-tuned from RNA-seq dataset 85.45% 75.86% 80.37% 59.39% 46.24% 52.00%\nClair3 DNA-seq model 53.25% 95.07% 68.27% 36.22% 73.19% 48.46%\nHG005\nClair3-RNA 94.27% 90.09% 92.13% 61.90% 70.49% 65.92%\nClair3-RNA with phasing enabled 95.88% 92.07% 93.93% 74.80% 75.58% 75.19%\nTable R6. Performance of Clair3 DNA model and Clair3 fine-tuned model in various datasets.\nSimiarly, a recent paper: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-\n023-02923-y claims that by processing the alignment file properly, the DNA-based variant\ncaller approach may achieve good accuracy for PacBio long-read IsoSeq data.\nRE16: We thank the reviewer for mentioning the work. Another reviewer had the same comment,\nand we have replied with benchmarking results in RE5. We have benchmarked the DNA-seq callers\nDeepVariant and Clair3 following the practices outlined in the paper. The results indicate that,\nalthough these methods repurposed the DNA-seq callers for RNA-seq after manipulating the\nalignment, the callers that natively designed for lrRNA-seq such as Clair3-RNA and DeepVariant-\nRNA performed much better.\n3. When Clair3-RNA utilizes the REDIportal database for RNA editing site tagging, does the\nstrand information take into account? For example, whether A-to-G or T-to-C alterations are\nboth allowed, or one of them will be eliminated if the read alignments clearly indicate the\npileup region corresponds to a specific strand RNA?\nRE17: We thank the reviewer for the question. Clair3-RNA tags a variant as RNA editing using\nREDIportal if both the called reference allele and alternative allele match exactly the corresponding\nrecord in the REDIportal database. Worth noticing, the REDIportal database categorizes all A-to-G\nediting events on the forward strand (denoted by “+” in the 5th column) and T-to-C editing events on\nthe reverse strand (denoted by “-” in the 5th column). Inspired by the reviewer’s question, we tried to\nmodify the RNA editing tagging logic that requires strand matching when using REDIportal for\ntagging. As shown in the table (Table R7) below, incorporating strand information resulted in a slight\ndecrease in precision and a slight increase in recall across all datasets. However, the overall F1-\nscore remained comparable to the results without considering strand information, indicating that only\na very small fraction of RNA editing events were affected by strand matching.\nConsider strand DP≥4 and AD≥2\nCaller Platform Dataset Sample matching in RNA SNP\nediting tagging TRUTH.FP TRUTH.FN TRUTH.TP QUERY.TP Precision Recall F1-score\nHG004 √ 6,236 10,210 152,814 152,822 96.08% 93.74% 94.89%\nIso-Seq HG004 6,100 10,214 152,810 152,818 96.16% 93.73% 94.93%\n(BCM) HG005 √ 8,377 12,972 224,196 224,215 96.40% 94.53% 95.46%\nPacBio\nHG005 8,175 12,972 224,196 224,215 96.48% 94.53% 95.50%\nHG004 √ 13,879 25,214 275,377 275,395 95.20% 91.61% 93.37%\nClair3-RNA MAS-Seq\nHG004 13,806 25,222 275,369 275,387 95.23% 91.61% 93.38%\nHG004 √ 6,900 10,706 89,022 89,023 92.81% 89.26% 91.00%\nHG004 6,887 10,711 89,017 89,018 92.82% 89.26% 91.00%\nONT\ndRNA004 HG004 √ 7,366 10,698 100,096 100,096 93.15% 90.34% 91.72%\nHG004 7,349 10,701 100,093 100,093 93.16% 90.34% 91.73%\nTable R7. Performance comparison of Clair3-RNA with or without considering strand matching in\nRNA editing tagging.\n4. Two minor concerns:\nFigure 8 is referred to between Figure 3 and Figure 4 in the text. The authors may consider\nreordering the figures or moving Figure 8 to the supplementary file, so the figure reference\ncan be in order.\nRE18: We have updated the manuscript accordingly.\nIn the speed comparison, why would Clair3-RNA be so much faster than Clair3?\nRE19: Clair3-RNA ran faster than Clair3 mainly for two reasons:\n1) Clair3-RNA used the pileup input representation, which is much smaller than Clair3’s full-\nalignment image-alike representation. The input size of the pileup is 33 flanking positions times 18\npileup features, and it is much smaller than the Clair3 full-alignment network, where the input is\n33x8x89 (594 integers versus 23,496 integers). The full-alignment is 5-8 times slower than the pileup\nmodel in network inference. 2) Some parameter optimizations have been used in Clair3-RNA, to avoid\nintroducing massive candidates into the neural network, such as the minimum coverage, the minimum\nallele depth, and the Indel candidate selection strategies for lrRNA-seq data, thus reducing the\nnumber of variant candidates that need to be processed by neural network.\nReviewer #2 (Remarks on code availability):\nDocumentation is comprehensive and easy to follow.\nRE20: We thank the reviewer for his/her efforts on checking the code and documentation.\nReviewer: #3\n----------------------------------------------------------\nThe manuscript introduces Clair3-RNA, a new deep learning method for small variant\ndetection from long-read RNA-seq data (lrRNA-seq). Clair3-RNA achieves strong\nperformance in several ONT and PacBio benchmarks, substantially outperforming other\nbenchmarked tools. The manuscript is very well written and includes numerous well-\ndesigned figures and tables. Notably, the authors performed extensive and informative\nexperiments to examine the impact of key data attributes on performance (e.g., genomic\ncontext, coverage, zygosity switching), which were also thoughtfully incorporated into\nmethod design.\nMajor concerns:\n* The HG004 and HG005 benchmark results appear to be only reported for the entire\ngenome in each sample. While the authors mention that they left out chromosome 20 during\nmodel training on the HG002 genome, I could not find the results specific only to\nchromosome 20 for HG004 and HG005. Since the HG002 sample (son) is related to the\nHG004 sample (mother), the model has been trained on variants which also appear in the\ntest set; it is hence important to test performance specifically on variants unseen during\ntraining.\nRE21: We thank the reviewer for the suggestion. The SNP and Indel performance of each dataset on\nchromosome 20 are provided in the table (Table R8 and Table R9) below. Across the HG002, HG004,\nand HG005 datasets, Clair3-RNA achieved an F1-score of at least 91.80% on PacBio Iso-Seq and\nMAS-Seq datasets and 90.88% on the ONT dRNA004 dataset, respectively. The results indicate that\nthe performance on chromosome 20 is consistent with the WGS results. Meanwhile, while HG004 is\nbiologically related to HG002, HG005 is part of a different family (Chinese trio family HG005/6/7).\nAnd we did not observe any performance abnormality when applying the pre-trained model trained\non HG002 to HG005.\nOnly for the sake of argument, the results of HG002 using the model trained on HG002 are also\nshown in the table (Table R10) below. The results have shown no performance increase in HG002\nwhen compared to the performance of HG004 and HG005. This aligns with the observations from the\nother tools in the Clair series that the model is not overfitting to specific samples.\nRegarding zygosity on chr20 Disregarding zygosity on chr20\nRead Allele SNP performance SNP performance SNP performance SNP performance\nPlatform Dataset Sample\nCoverage depth (apply phasing ✓) (apply phasing ✘) (apply phasing ✓) (apply phasing ✘)\nPrecision Recall F1-scorePrecision Recall F1-scorePrecision Recall F1-scorePrecision Recall F1-score\nHG002 93.74% 94.20% 93.97% 92.92% 93.06% 92.99% 96.29% 96.75% 96.52% 95.67% 95.80% 95.73%\nIso-Seq\nHG004 96.76% 94.42% 95.58% 95.93% 93.27% 94.58% 99.70% 97.28% 98.48% 99.12% 96.37% 97.72%\n(BCM)\nPacBio HG005 97.18% 95.46% 96.31% 96.74% 94.10% 95.40% 99.00% 97.25% 98.12% 98.83% 96.13% 97.46%\nHG002 93.41% 92.91% 93.16% 93.07% 90.57% 91.80% 96.65% 96.11% 96.38% 96.69% 94.08% 95.37%\nDP≥4 AD≥2 MAS-Seq\nHG004 95.30% 92.08% 93.66% 94.59% 89.62% 92.04% 99.22% 95.85% 97.51% 98.59% 93.40% 95.92%\nHG002 93.41% 92.75% 93.08% 92.38% 89.43% 90.88% 95.74% 95.07% 95.40% 95.18% 92.15% 93.64%\nONT drna004 HG004 94.13% 92.49% 93.30% 93.54% 89.36% 91.40% 96.07% 94.39% 95.23% 95.79% 91.51% 93.60%\nHG005 94.69% 93.52% 94.10% 94.08% 91.10% 92.57% 96.19% 95.01% 95.60% 96.00% 92.96% 94.46%\nHG002 95.71% 95.89% 95.80% 94.80% 94.50% 94.65% 96.44% 96.62% 96.53% 95.76% 95.46% 95.61%\nIso-Seq\nHG004 98.83% 96.00% 97.40% 98.15% 94.54% 96.31% 99.89% 97.03% 98.44% 99.21% 95.57% 97.36%\n(BCM)\nPacBio HG005 99.36% 96.86% 98.10% 99.12% 95.97% 97.52% 99.60% 97.09% 98.33% 99.40% 96.24% 97.79%\nHG002 95.08% 95.24% 95.16% 95.08% 92.92% 93.99% 96.47% 96.62% 96.55% 96.76% 94.56% 95.65%\nDP≥8 AD≥2 MAS-Seq\nHG004 97.43% 94.28% 95.83% 96.89% 91.54% 94.14% 99.65% 96.43% 98.01% 99.23% 93.76% 96.42%\nHG002 95.57% 95.25% 95.41% 94.50% 93.05% 93.77% 96.97% 96.64% 96.80% 96.10% 94.63% 95.36%\nONT drna004 HG004 95.95% 94.58% 95.26% 95.31% 91.57% 93.40% 96.76% 95.38% 96.07% 96.74% 92.93% 94.80%\nHG005 96.10% 95.95% 96.02% 95.25% 94.13% 94.69% 96.90% 96.75% 96.82% 96.55% 95.41% 95.98%\nHG002 95.58% 96.22% 95.90% 94.75% 95.06% 94.90% 96.20% 96.84% 96.52% 95.52% 95.83% 95.67%\nIso-Seq\nHG004 99.03% 95.92% 97.45% 98.27% 94.65% 96.42% 99.93% 96.79% 98.34% 99.24% 95.59% 97.38%\n(BCM)\nPacBio HG005 99.39% 96.97% 98.17% 99.34% 96.06% 97.67% 99.58% 97.16% 98.35% 99.53% 96.24% 97.86%\nHG002 95.45% 95.77% 95.61% 95.53% 93.47% 94.49% 96.55% 96.86% 96.70% 96.98% 94.89% 95.92%\nDP≥10 AD≥2 MAS-Seq\nHG004 97.78% 94.58% 96.15% 97.19% 91.79% 94.41% 99.78% 96.51% 98.12% 99.40% 93.88% 96.56%\nHG002 95.81% 95.55% 95.68% 95.01% 93.94% 94.47% 97.15% 96.89% 97.02% 96.58% 95.50% 96.04%\nONT drna004 HG004 96.05% 95.20% 95.62% 95.57% 92.62% 94.07% 96.82% 95.96% 96.39% 96.81% 93.82% 95.29%\nHG005 96.32% 96.79% 96.55% 95.68% 95.45% 95.57% 96.92% 97.39% 97.16% 96.60% 96.36% 96.48%\nHG002 96.13% 96.69% 96.41% 95.17% 95.70% 95.44% 96.74% 97.30% 97.02% 95.93% 96.47% 96.20%\nIso-Seq\nHG004 99.01% 96.42% 97.70% 98.65% 95.66% 97.13% 99.93% 97.32% 98.61% 99.65% 96.63% 98.11%\n(BCM)\nPacBio HG005 99.47% 97.35% 98.40% 99.51% 96.92% 98.20% 99.66% 97.54% 98.59% 99.71% 97.11% 98.39%\nHG002 95.70% 96.50% 96.10% 95.96% 94.52% 95.24% 96.81% 97.60% 97.20% 97.43% 95.97% 96.70%\nDP≥10 AD≥4 MAS-Seq\nHG004 97.80% 95.64% 96.71% 97.38% 93.41% 95.35% 99.86% 97.65% 98.74% 99.65% 95.59% 97.58%\nHG002 95.92% 96.19% 96.05% 95.25% 95.09% 95.17% 97.23% 97.49% 97.36% 96.78% 96.62% 96.70%\nONT drna004 HG004 96.10% 96.16% 96.13% 95.77% 94.27% 95.02% 96.88% 96.94% 96.91% 97.03% 95.51% 96.26%\nHG005 96.38% 97.28% 96.83% 95.87% 96.29% 96.08% 97.00% 97.90% 97.45% 96.80% 97.22% 97.01%\nTable R8. SNP performance of Clair3-RNA on chromosome 20 in various datasets.\nRegarding zygosity on chr20 Disregarding zygosity on chr20\nRead Allele Indel performance Indel performance Indel performance Indel performance\nPlatform Dataset Sample\nCoverage depth (apply phasing ✓) (apply phasing ✘) (apply phasing ✓) (apply phasing ✘)\nPrecision Recall F1-scorePrecision Recall F1-scorePrecision Recall F1-scorePrecision Recall F1-score\nIso-Seq HG002 81.51% 74.18% 77.67% 76.48% 67.89% 71.93% 90.05% 82.04% 85.86% 87.52% 77.80% 82.37%\n(BCM) HG004 85.78% 77.54% 81.45% 80.32% 67.51% 73.36% 92.03% 83.20% 87.39% 91.44% 76.85% 83.51%\nPacBio HG005 86.10% 80.44% 83.17% 83.02% 74.20% 78.36% 92.96% 86.85% 89.80% 91.32% 81.62% 86.20%\nDP≥4 AD≥2 MAS-Seq HG002 78.37% 69.54% 73.69% 75.81% 61.93% 68.17% 87.49% 77.77% 82.35% 88.67% 72.48% 79.76%\nHG004 78.06% 65.46% 71.21% 74.10% 58.33% 65.28% 87.28% 73.31% 79.69% 86.43% 68.14% 76.20%\nHG002 85.26% 51.70% 64.37% 79.86% 48.19% 60.11% 89.12% 54.04% 67.28% 86.22% 52.03% 64.89%\nONT drna004 HG004 90.45% 56.21% 69.34% 82.55% 49.58% 61.95% 92.27% 57.34% 70.73% 89.62% 53.82% 67.26%\nHG005 85.55% 64.04% 73.24% 87.39% 58.77% 70.28% 89.45% 66.96% 76.59% 91.74% 61.70% 73.78%\nIso-Seq HG002 84.92% 77.35% 80.96% 78.27% 71.55% 74.76% 90.90% 82.80% 86.66% 87.07% 79.61% 83.17%\n(BCM) HG004 89.63% 81.42% 85.33% 81.71% 69.08% 74.87% 93.33% 84.80% 88.86% 91.44% 77.30% 83.78%\nPacBio HG005 90.75% 84.18% 87.34% 84.68% 75.00% 79.55% 95.09% 88.20% 91.52% 91.29% 80.85% 85.75%\nDP≥8 AD≥2 MAS-Seq HG002 80.51% 70.91% 75.40% 77.43% 63.38% 69.70% 87.41% 77.09% 81.93% 88.19% 72.24% 79.42%\nHG004 81.71% 66.57% 73.37% 75.13% 59.17% 66.20% 88.89% 72.45% 79.83% 85.29% 67.17% 75.15%\nHG002 86.67% 50.78% 64.04% 80.70% 48.17% 60.33% 89.33% 52.34% 66.01% 85.53% 51.05% 63.93%\nONT drna004 HG004 89.20% 56.07% 68.86% 81.82% 51.80% 63.44% 90.91% 57.14% 70.18% 88.07% 55.76% 68.28%\nHG005 87.06% 65.06% 74.47% 87.23% 59.42% 70.69% 89.55% 66.91% 76.60% 90.96% 61.96% 73.71%\nIso-Seq HG002 86.02% 79.07% 82.40% 79.23% 72.53% 75.73% 91.06% 83.71% 87.23% 86.31% 79.03% 82.51%\n(BCM) HG004 91.04% 80.33% 85.35% 80.95% 68.55% 74.24% 95.28% 84.10% 89.34% 91.43% 77.42% 83.84%\nPacBio HG005 91.07% 84.44% 87.63% 83.09% 73.86% 78.20% 94.29% 87.42% 90.72% 90.07% 80.07% 84.78%\nDP≥10 AD≥2 MAS-Seq HG002 81.22% 71.21% 75.89% 78.25% 63.54% 70.13% 87.63% 76.89% 81.91% 88.31% 71.71% 79.15%\nHG004 83.30% 65.19% 73.14% 76.67% 59.45% 66.97% 90.61% 70.94% 79.57% 86.67% 67.21% 75.71%\nHG002 85.78% 50.58% 63.64% 79.61% 47.95% 59.85% 88.73% 52.31% 65.82% 84.47% 50.88% 63.50%\nONT drna004 HG004 88.55% 55.89% 68.53% 81.71% 51.54% 63.21% 90.36% 57.03% 69.93% 87.20% 55.00% 67.45%\nHG005 86.89% 64.63% 74.13% 87.13% 59.13% 70.45% 89.07% 66.26% 75.99% 90.06% 61.11% 72.81%\nIso-Seq HG002 87.42% 83.85% 85.60% 80.60% 77.33% 78.93% 91.84% 88.09% 89.93% 87.29% 83.76% 85.49%\n(BCM) HG004 93.03% 84.93% 88.80% 83.08% 73.25% 77.86% 97.51% 89.04% 93.08% 94.03% 82.89% 88.11%\nPacBio HG005 93.39% 88.89% 91.08% 85.27% 80.29% 82.71% 96.50% 91.85% 94.12% 92.25% 86.86% 89.47%\nDP≥10 AD≥4 MAS-Seq HG002 83.71% 74.85% 79.03% 80.61% 67.50% 73.47% 89.72% 80.29% 84.74% 90.16% 75.50% 82.18%\nHG004 86.43% 70.81% 77.85% 79.22% 64.89% 71.35% 93.44% 76.58% 84.17% 88.53% 72.52% 79.73%\nHG002 86.07% 52.74% 65.41% 79.80% 50.00% 61.48% 89.05% 54.57% 67.67% 84.24% 52.78% 64.90%\nONT drna004 H\nH\nG\nG\n0\n0\n0\n0\n4\n5\n8\n8\n8\n7\n.\n.\n9\n7\n6\n8\n%\n%\n5\n6\n8\n6\n.\n.\n0\n6\n0\n7\n%\n%\n7\n7\n0\n5\n.\n.\n2\n7\n2\n8\n%\n%\n8\n8\n1\n7\n.\n.\n7\n5\n1\n0\n%\n%\n5\n6\n4\n0\n.\n.\n2\n4\n5\n9\n%\n%\n6\n7\n5\n1\n.\n.\n2\n5\n1\n3\n%\n%\n9\n9\n0\n0\n.\n.\n8\n0\n0\n0\n%\n%\n5\n6\n9\n8\n.\n.\n2\n3\n0\n5\n%\n%\n7\n7\n1\n7\n.\n.\n6\n7\n7\n0\n%\n%\n8\n9\n7\n0\n.\n.\n2\n4\n0\n8\n%\n%\n5\n6\n7\n2\n.\n.\n8\n5\n9\n5\n%\n%\n6\n7\n9\n3\n.\n.\n5\n9\n9\n7\n%\n%\nTable R9. Indel performance of Clair3-RNA on chromosome 20 in various datasets.\nSNP Indel\nRead Allele\nSample Dataset Regarding zygosity Disregarding zygosity Regarding zygosity Disregarding zygosity\nCoverage depth\nPrecision Recall F1-scorePrecision Recall F1-scorePrecision Recall F1-scorePrecision Recall F1-score\nIso-Seq(BCM) 95.34% 95.52% 95.43% 97.40% 97.58% 97.49% 85.57% 76.13% 80.57% 92.83% 82.60% 87.41%\nDP≥4 AD≥2 MAS-Seq 95.67% 94.01% 94.84% 97.89% 96.18% 97.03% 83.44% 69.37% 75.76% 92.37% 76.81% 83.87%\ndrna004 93.10% 90.08% 91.56% 95.72% 92.61% 94.14% 85.95% 49.93% 63.16% 89.90% 52.23% 66.07%\nIso-Seq(BCM) 96.91% 97.33% 97.12% 97.61% 98.04% 97.82% 88.66% 80.53% 84.40% 93.54% 84.98% 89.06%\nDP≥8 AD≥2 MAS-Seq 97.41% 96.08% 96.74% 98.29% 96.95% 97.61% 85.36% 71.65% 77.91% 92.85% 77.95% 84.75%\nHG002 drna004 94.56% 93.33% 93.94% 96.16% 94.90% 95.53% 86.25% 50.55% 63.74% 89.34% 52.36% 66.03%\n(WGS) Iso-Seq(BCM) 97.17% 97.63% 97.40% 97.75% 98.21% 97.98% 89.42% 81.70% 85.39% 93.81% 85.72% 89.59%\nDP≥10 AD≥2 MAS-Seq 97.70% 96.51% 97.10% 98.42% 97.23% 97.82% 86.05% 72.42% 78.65% 93.14% 78.40% 85.14%\ndrna004 94.81% 94.12% 94.46% 96.31% 95.61% 95.96% 86.39% 50.60% 63.82% 89.32% 52.32% 65.99%\nIso-Seq(BCM) 97.50% 98.23% 97.86% 98.10% 98.83% 98.46% 91.21% 86.40% 88.74% 95.20% 90.19% 92.63%\nDP≥10 AD≥4 MAS-Seq 97.93% 97.31% 97.62% 98.66% 98.04% 98.35% 87.95% 76.19% 81.65% 94.50% 81.87% 87.73%\ndrna004 95.01% 95.26% 95.14% 96.53% 96.78% 96.66% 86.98% 53.76% 66.45% 89.71% 55.45% 68.54%\nTable R10. SNP and Indel performance of Clair3-RNA on HG002 WGS.\n* Additional metrics describing the generated training data are needed to evaluate the\napproach (i.e. how balanced/diverse was the training data), including: how many true\nvariants, artifacts, and editing sites were included; what types of artifacts are expected in this\ndata and were any artifact sub-categories specifically identified/selected; what was the\ndistribution of allele frequencies in the training data.\nRE22: We thank the reviewer for the suggestions. We have updated the statistics of the training\nsamples, including the distribution of variant types and output task categories.in Supplementary\nFigure 1. We have depicted the overall distribution of the training data across four categories:\nhomozygous germline variants, heterozygous germline variants, sequencing artifacts (homozygous\nreference), and RNA editing sites. Additionally, we further classified the artifacts into ACGT,\ninsertion, and deletion artifacts based on the signal of the candidates with the maximum alternate\nread support.\nWe also categorized the training dataset according to the 21-genotype output task, assigning each\nvariant to one of the 21 genotypes based on the reference allele, alternate allele, and zygosity.\nFor RNA editing sites, they were assigned as AA or TT (indicating a reference allele A for A-to-G editing\nor T for T-to-C editing on the reverse strand). The allele frequency distribution of the different\ntraining datasets is also illustrated in Supplementary Figure 1.\na\nb\nc\na\nb\nc\nFigure R6. The figure shows the statistics of the training samples, including the distribution of variant\ntypes, and output task categories, for (a) PacBio Iso-Seq HG002 and Mas-Seq HG002, and (b) ONT\ncDNA and ONT dRNA004 datasets. (c) shows the allelic fraction distribution of different training\ndatasets.\n* In the paragraph starting with line 638, the authors mention that they explored including\nRNA edit site prediction as a separate task or category within the zygosity task but these\napproaches did not produce accurate results (why?). However, when integrating editing sites\ninto the homozygous reference category, the model performed well (why?). This discrepancy\nin performance could point to a limitation of the input representation or the training dataset\ncomposition and should be further examined. \nRE23: We thank the reviewer for the comment. A major challenge is that RNA editing sites used for\ntraining are not perfect, making training a separate model task for a comprehensive RNA editing site\ndetection challenging. For instance, as depicted in Figure 6a, editing sites are collected under\nspecific read coverage and allele depth conditions, which may result in the remaining editing sites\nbeing assigned to the artifacts category, thereby affecting model training. During model training,\nidentifying all RNA editing sites with lower read coverage (e.g., A-to-G or T-to-C allele depth <8)\nremains challenging, as these less-confident editing sites may also include sequencing artifacts.\nAssigning editing sites with AD >8 as RNA editing and those with AD <8 as artifacts could lead to\nconfusion during model training.\nIn contrast, by assigning all editing sites to homozygous reference, the neural network only needs to\ndistinguish true germline variants from non-variants (can be either RNA editing or sequencing\nartifacts). Empirically, assigning all editing sites to homozygous reference has yield better results.\nHow many edit sites versus reference sites versus true SNP sites were included during\ntraining and evaluation?\nGiven the features provided to the model as input, what distinguishing information is currently\navailable to the model to differentiate between RNA edit sites and true variants?\nRE24: The distribution of RNA editing in model training is depicted in RE22. We believe the Clair3-\nRNA’s model has learned the latent features to tell a true variant from RNA editing sites. As shown in\nFigure 6b and the Result - A-to-G and T-to-C false positives analysis section, 67.07%/93.27% of the\nediting sites were identified by the neural network, and 29.87%/26.94% in the remainder were further\nfiltered by REDIportal in turn for ONT and PacBio platform, respectively.\n* The paper includes comparisons to only a single method designed specifically for lrRNA-\nseq variant calling (namely, LongcallR). A few other methods should be cited and\nbenchmarked to enable a better performance assessment — e.g., CTAT-Mutations,\nisoLASER, LongcallR-nn, or the recent lrRNA-seq model of DeepVariant.\nRE25: We thank the reviewer for raising this point. To ensure a comprehensive performance\nassessment, we have added benchmarking results for these methods. One exception is that we were\nunable to run isoLASER following the GitHub instructions using either conda or singularity\ninstallation. The authors mentioned on GitHub that they plan to update the source code. We have\nraised this issue with the authors but have not yet received a response\n(https://github.com/gxiaolab/isoLASER/issues/5). Other users have also reported similar issues when\nattempting to run isoLASER (https://github.com/gxiaolab/isoLASER/issues/4 and\nhttps://github.com/gxiaolab/isoLASER/issues/3). For the other methods, we have included the\ncommand lines used in the Supplementary Methods - Command line usedsection.\nThe results shown in the table (Table R11) below demonstrate that Clair3-RNAoutperforms CTAT-\nMutations, LongcallR-nn, and DeepVariant lrRNA-seq models across multiple datasets on both\nPacBio and ONT platforms.\nSNP (DP≥4 and AD≥2) Indel (DP≥4 and AD≥2)\nPlatform Dataset Sample Caller\nPrecision Recall F1-score Precision Recall F1-score\nHG004 Clair3-RNA with phasing 96.81% 94.76% 95.77% 83.49% 76.02% 79.58%\nHG004 Clair3-RNA 96.16% 93.73% 94.93% 80.74% 70.21% 75.11%\nHG004 DeepVariant (RNA model) 89.87% 93.12% 91.47% 78.89% 67.53% 72.77%\nHG004 CTAT-Mutations 80.99% 93.22% 86.67% 38.79% 72.15% 50.45%\nIso-Seq HG004 LongcallR-nn 96.10% 62.88% 76.02% - - -\n(BCM) HG005 Clair3-RNA with phasing 97.13% 95.81% 96.47% 88.40% 82.99% 85.61%\nHG005 Clair3-RNA 96.48% 94.53% 95.50% 85.47% 77.69% 81.40%\nPacBio HG005 DeepVariant (RNA model) 91.79% 93.78% 92.77% 83.33% 76.47% 79.75%\nHG005 CTAT-Mutations 78.20% 95.22% 85.87% 41.38% 79.67% 54.47%\nHG005 LongcallR-nn 95.91% 68.05% 79.61% - - -\nHG004 Clair3-RNA with phasing 95.81% 93.44% 94.61% 79.07% 67.28% 72.70%\nHG004 Clair3-RNA 95.23% 91.61% 93.38% 76.48% 60.71% 67.69%\nMAS-Seq HG004 DeepVariant (RNA model) 94.29% 89.68% 91.92% 83.28% 60.34% 69.98%\nHG004 CTAT-Mutations 70.84% 89.32% 79.01% 13.58% 62.14% 22.29%\nHG004 LongcallR-nn 94.29% 66.20% 77.79% - - -\nHG004 Clair3-RNA with phasing 94.04% 92.35% 93.19% 87.36% 52.50% 65.59%\nHG004 Clair3-RNA 92.82% 89.26% 91.00% 83.61% 47.19% 60.33%\nHG004 LongcallR-nn 93.35% 65.67% 77.10% - - -\nONT dRNA004\nHG005 Clair3-RNA with phasing 94.23% 92.90% 93.56% 88.55% 58.74% 70.63%\nHG005 Clair3-RNA 93.16% 90.34% 91.73% 85.64% 53.76% 66.05%\nHG005 LongcallR-nn 93.88% 65.83% 77.40% - - -\nTable R11. Performance comparison of different callers.\nMinor:\n* In addition to the reported parallel runtime performance, runtime should also be reported\nfor each method in single-threaded mode (since each method may employ different parallel\nprocessing schemes).\nRE26: The runtime for each method in different datasets using a single thread is reported in\nSupplementary Table 11. For Clair3-RNA, Clair3, and LongcallR, the “--threads=1” option was set.\nFor DeepVariant, the “--num_shards=1” option was used.\nClair3 RNA Clair3 DeepVariant LongcallR\nTesting Peak Peak Peak Peak\nRuntime Runtime Runtime\ndataset memory memory memory Runtime (min) memory\n(min) (min) (min)\n(GB) (GB) (GB) (GB)\nPacBio Iso-Seq 303 0.6 383 3.1 442 9.3 120 3.8\nPacBio MAS-Seq 1390 0.7 1471 23.3 1007 9.4 1717 7.0\nONT dRNA004 860 0.7 858 15.3 619 5.1 490 4.8\nONT cDNA 960 0.6 1795 18.5 1040 5.1 3748 4.9\nTable R12. Runtime and memory usage of different callers using a single thread.\nReviewer #3 (Remarks on code availability):\nThe code provides extensive documentation with detailed installation/usage instructions; I\nwas able to run the code (using docker).\nRE27: We thank the reviewer for their valuable time in testing our method."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have performed additional benchmarking and evaluations to address the comments from the previous review. It is great to see that the splice junction analysis has allowed them to improve the accuracy of the method further. I have no further comments.", "reviewer_2": "The authors have nicely addressed all my concerns.", "reviewer_3": "I thank the authors for providing useful additional information and benchmarks in response to the reviewers; however, several major concerns from my previous review have not been substantially addressed:\n* I could not find updates to the manuscript regarding the concern that the HG004 testing dataset is significantly similar to the HG002 training dataset (as they correspond to the mother and son genomes, respectively). As such, key model evaluation is still performed on a test dataset, which significantly overlaps with the training data — this is not addressed/discussed in the paper. Even if the model generalizes to other genomes, this evaluation approach does not follow standard ML practice. The authors could update the manuscript motivating the use of HG004 for evaluation (e.g. if no other datasets were available) and clarify to the readers that these genomes are related, including a breakdown of the number of shared examples between HG002 and HG004, as well as between HG002 and HG005 (although HG002 and HG005 are not related, they may still share variant and artifact sites). The model could then be evaluated only on the examples that are not included in the training dataset (from all the chromosomes). Based on the chr20 tables, showing results just on chr20 would likely miss important information. E.g., the performance on HG002 is consistently significantly lower — what could explain this result?\n* The new Supplementary Figure 1 is very informative. However, there does appear to be a lack of balance across several training categories (e.g. allele fraction, artifact indels vs true indels) — it would be useful for authors to comment on the performance implications/limitations of this training dataset design.\n* For the model to learn the latent features to distinguish a true variant from an RNA editing site, there should be some distinguishing information provided in the input — what input features can be used to differentiate these sites?\n* The latest benchmarks should be included in the manuscript (e.g. with CTAT and DeepVariant in table R11). The performance shown for DeepVariant is significantly lower than what is reported in their Github case study for HG004. What is the reason for this difference? E.g., if a method is not designed to be used in a particular regime (e.g. on some variants included in the evaluation), performance could be shown as a function of different regimes."}, "author_response": "Responses to\nComments of reviewers\nClair3-RNA: A deep learning-based small variant caller for long-read RNA sequencing\ndata\nNCOMMS-24-75685A\nReviewer: #3\n----------------------------------------------------------\nI could not find updates to the manuscript regarding the concern that the HG004 testing\ndataset is significantly similar to the HG002 training dataset (as they correspond to the\nmother and son genomes, respectively). As such, key model evaluation is still performed on\na test dataset, which significantly overlaps with the training data — this is not\naddressed/discussed in the paper. Even if the model generalizes to other genomes, this\nevaluation approach does not follow standard ML practice. The authors could update the\nmanuscript motivating the use of HG004 for evaluation (e.g. if no other datasets were\navailable) and clarify to the readers that these genomes are related, including a breakdown\nof the number of shared examples between HG002 and HG004, as well as between HG002\nand HG005 (although HG002 and HG005 are not related, they may still share variant and\nartifact sites). The model could then be evaluated only on the examples that are not included\nin the training dataset (from all the chromosomes).\nRE1: We sincerely thank the reviewer for their valuable comments for us to improve our manuscript.\nWe have updated the manuscript in the Result section on the motivation of using HG004 and HG005\nsamples for benchmarking. (line 124-128)\n“For benchmarking, we selected two GIAB samples: HG004 (the mother of HG002) and HG005\n(from an unrelated Chinese trio family). This dual-sample approach enables comprehensive\nevaluation while controlling for potential biases from biological relationship, as HG004 shares\ngenetic background with HG002.”\nFollowing the reviewer’s suggestion, we conducted additional benchmarking to evaluate variants and\nartifacts in HG004 and HG005 while explicitly excluding those present in HG002, which were used\nduring model training, to prevent any potential evaluation bias. As shown in Table R1 below, this\nexclusion removed ~350k - ~651k candidates from HG004 or HG005 benchmarking, depending on\nthe read coverage of the sequencing datasets. The numbers confirm the reviewer’s observation that\neven biologically unrelated samples can share certain variants and artifacts. After this adjustment, as\nshown in Table R2 below, Clair3-RNA achieved with average SNP F1-scores of 94.50%, 90.80%, and\n90.43% in PacBio Iso-Seq, MAS-Seq, and ONT drna004 datasets, respectively. The model showed\nfurther improvement when phasing was enabled, achieving 95.74%, 93.10%, and 92.89% on the same\ndatasets. After removing HG002 variants, the precision remains comparable to the original\nHG004/HG005 benchmarks in all callable regions (average 95.31% versus 94.77% before phasing,\nand 96.24% versus 95.61% after phasing). The recall decreased by an average of 4.14% and 2.72%\nbefore and after phasing—which is expected due to an average ~43% reduction in benchmarked true\npositives without a substantial proportional change in false negatives. The results indicate that\nClair3-RNA did not overfit specific samples. We have added a “Performance analysis after removing\nshared variants between samples” section in the Result section to report the results. (line 497-510)\nWe would like to thank the reviewer once again for the valuable suggestion. We believe this critical\nanalysis provides more comprehensive evidence of the robustness of our method.\nTotal Candidates in Total Total SNP variants in Indel variants in\ncandidates specific sample SNP variants Indel variants specific sample specific sample\nPlatform Dataset Sample\n(DP≥4 && (HG004 or HG005) (DP≥4 && (DP≥4 && (HG004 or HG005) (HG004 or HG005)\nAF≥0.05) and not in HG002 AF≥0.05) AF≥0.05) and not in HG002 and not in HG002\nHG002 3,957,587 - 434,400 65,233 - -\nIso-Seq (BCM) HG004 1,354,990 1,003,846 163,024 24,634 80,263 13,613\nPacBio HG005 1,965,315 1,514,321 237,168 31,688 158,078 21,442\nHG002 6,384,773 - 781,823 119,398 - -\nMAS-Seq\nHG004 2,739,639 2,088,447 300,591 45,857 141,620 25,681\nHG002 3,582,623 - 128,755 18,559 - -\nONT dRNA004 HG004 2,895,297 2,318,438 99,728 13,835 52,383 9,371\nHG005 3,212,196 2,608,619 110,794 13,525 75,672 9,819\nTable R1. Statistics of shared variant and candidate sites between HG004/HG005 and HG002.\nRead Allele Variant\nDataset Sample Caller Precision Recall F1-score TRUTH.FP TRUTH.FNTRUTH.TP QUERY.TP\ncoverage depth type\nClair3-RNA with phasing 97.56% 92.74% 95.09% 1,860 5,825 74,438 74,453\nHG004\nIso-Seq Clair3-RNA 96.75% 91.08% 93.83% 2,448 7,151 72,973 72,981\n(BCM) Clair3-RNA with phasing 97.67% 95.14% 96.39% 3,596 7,679 150,399 150,421\nHG005\nClair3-RNA 96.97% 93.42% 95.17% 4,610 10,394 147,679 147,698\nClair3-RNA with phasing 96.12% 90.26% 93.10% 5,166 13,791 127,829 127,851\nSNP MAS-Seq HG004\nClair3-RNA 95.16% 86.83% 90.80% 6,312 18,818 124,058 124,076\nClair3-RNA with phasing 94.62% 90.06% 92.28% 2,682 5,208 47,175 47,177\nHG004\nClair3-RNA 93.36% 85.70% 89.36% 3,179 7,460 44,696 44,697\ndRNA004\nClair3-RNA with phasing 95.25% 91.81% 93.50% 3,467 6,195 69,477 69,476\nHG005\nClair3-RNA 94.29% 88.84% 91.49% 4,059 8,420 67,047 67,047\nDP≥4 AD≥2\nClair3-RNA with phasing 86.28% 69.40% 76.93% 1,505 4,165 9,448 9,466\nHG004\nIso-Seq Clair3-RNA 83.95% 62.07% 71.37% 1,599 5,085 8,321 8,365\n(BCM) Clair3-RNA with phasing 90.76% 82.28% 86.31% 1,798 3,799 17,643 17,665\nHG005\nClair3-RNA 88.83% 76.85% 82.41% 2,045 4,889 16,227 16,267\nClair3-RNA with phasing 82.40% 57.18% 67.51% 3,142 10,996 14,685 14,708\nIndel MAS-Seq HG004\nClair3-RNA 80.27% 48.96% 60.82% 3,112 13,182 12,643 12,659\nClair3-RNA with phasing 87.83% 41.03% 55.93% 534 5,526 3,845 3,855\nHG004\nClair3-RNA 83.57% 33.93% 48.26% 615 6,088 3,126 3,128\ndRNA004\nClair3-RNA with phasing 90.19% 54.02% 67.57% 578 4,515 5,304 5,315\nHG005\nClair3-RNA 87.37% 47.95% 61.92% 671 5,035 4,638 4,640\nTable R2. Performance evaluation in HG004 and HG005 datasets after excluding HG002-shared\nvariants and artifacts.\nBased on the chr20 tables, showing results just on chr20 would likely miss important\ninformation. E.g., the performance on HG002 is consistently significantly lower — what could\nexplain this result?\nRE2: We thank the reviewer for their question regarding performance variation across different\nsamples. These differences are expected due to inherent variability in read coverage and variant\nallelic fraction distributions, which are influenced by disparities in gene expression levels and allele-\nspecific expression. To quantify this, we analyzed read coverage distributions across datasets (Table\nR3 below) and observed that HG002 exhibits a higher proportion of low-coverage (DP≤10) regions\n(average 3.77%) compared to HG004 and HG005 (average 2.11%) in all three platforms. This\ndiscrepancy leads to a greater number of low-coverage candidates being included and called in\nmodel inference for HG002. For example, HG002 Iso-Seq had 53,633 candidates (1,042 SNPs and 77\nIndels were called) with DP ≤ 10, substantially more than HG004 (23,469 candidates; 816 SNPs, 48\nIndels) and HG005 (30,057 candidates; 956 SNPs, 59 Indels). Lower read coverage typically results\nin inferior variant calling performance, as it becomes more difficult to differentiate true variants from\nsequencing artifacts. As a result of these additional low-coverage candidates, Clair3-RNA called\nmore variants in HG002, ultimately contributing to its marginally lower overall F1-score compared\nto HG004 and HG005.\nDetected SNP Detected Indel Faction Faction Faction\nCandidates Candidates Candidates Base pair Base pair Base pair\nDataset Sample Variant Variant (DP≤6 / (DP≤8 / (DP≤10 /\n(4≤DP≤6) (4≤DP≤8) (4≤DP≤10) (1≤DP≤6) (1≤DP≤8) (1≤DP≤10)\n(4≤DP≤10) (4≤DP≤10) WGS) WGS) WGS)\nHG002 1,042 77 29,414 46,353 53,633 1,977,600 2,548,524 2,898,918 3.07% 3.95% 4.50%\nPacBio\nIso-Seq (BCM) HG004 816 48 14,215 21,178 23,469 1,106,909 1,377,991 1,533,038 1.72% 2.14% 2.38%\nchr20\nHG005 956 59 16,922 26,366 30,057 1,169,341 1,525,442 1,757,309 1.81% 2.37% 2.73%\nPacBio HG002 5,420 566 42,743 69,728 84,081 2,326,071 3,060,277 3,671,183 3.61% 4.75% 5.70%\nMAS-Seq\nchr20 HG004 3,012 323 26,402 41,356 48,125 1,493,716 1,968,293 2,297,052 2.32% 3.05% 3.56%\nHG002 3,875 493 32,802 52,827 63,390 536,124 648,421 716,242 0.83% 1.01% 1.11%\nONT\nHG004 1,851 222 25,549 40,881 49,148 431,252 496,952 549,349 0.67% 0.77% 0.85%\ndRNA004\nHG005 2,065 258 30,480 49,057 58,143 508,964 605,037 674,617 0.79% 0.94% 1.05%\nTable R3. The size and percentage of genomic regions and candidates at certain low read coverages\n(DP) across different samples.\nThe new Supplementary Figure 1 is very informative. However, there does appear to be a\nlack of balance across several training categories (e.g. allele fraction, artifact indels vs true\nindels) — it would be useful for authors to comment on the performance\nimplications/limitations of this training dataset design.\nRE3: We thank the reviewer for the valuable suggestions. We have updated the manuscript to include\na detailed explanation of our use of Focal loss to address the data and category imbalance issue. (line\n811-821)\n“The variant and artifact identification in our study was significantly affected by the coverage and\nerror rate distributions in the training datasets. As shown in Supplementary Figure 1, artifact\ncandidates substantially outnumbered SNPs and Indels during model training, and multi-allelic\nvariants were far less frequent than homozygous and heterozygous variants. To address data and\ncategory imbalance issues, we implemented Focal loss [1] during model training - an approach first\nintroduced in our previous work [2]. Focal loss improves upon standard cross-entropy by\nincorporating a modulating factor that automatically down-weights well-classified examples while\nmaintaining focus on challenging, underrepresented cases. This approach proved particularly\neffective for our application, as it ensured proper attention to rare but critical variant categories\ndespite the overwhelming prevalence of artifacts in the training data, ultimately leading to more\nstable model convergence and improved performance.”\nWe hope this explanation helps readers better understand the data and category imbalance\nchallenges inherent in model training and our Focal loss solution on it.\nFor the model to learn the latent features to distinguish a true variant from an RNA editing\nsite, there should be some distinguishing information provided in the input — what input\nfeatures can be used to differentiate these sites?\nRE4: We thank the reviewer for the question. Clair3-RNA utilizes the pileup features including the\nread depth, allele frequency, strand distribution, and the flanking reference sequences for variant\ndetection. These features have been successfully demonstrated to identify A-to-I RNA editing using\nRNA-Seq from a single sample in another study [3].\nTo systematically evaluate which features contribute most to distinguishing RNA editing sites from\ntrue variants, we performed supervised UMAP analysis across three datasets (HG002 PacBio Iso-\nSeq, PacBio MAS-Seq, and ONT dRNA004) using four tensor configurations:\n1. Full-context pileup tensor: the pileup tensor used by Clair3-RNA.\n2. Tensor with center position only: only the pileup information for the central variant/editing\nsite was used, masking flanking sequence information.\n3. Tensor without the Indel information: For all the positions, the Indel features are masked.\n4. Tensor after strand unification: All reverse strands were converted to the forward strand,\nremoving strand-specific information.\nAs shown in the UMAP visualization (Figure R1), we observe clear separation between RNA editing\nsites and true variants when using the full-context tensor, while no distinct clustering emerges when\nanalyzing only the center position. This finding demonstrates that flanking sequence information is\nessential for discriminating between RNA editing events and genuine variants. For the Indel\ninformation and strand information, less significant impacts were observed in the UMAP analysis.\nAnother notable observation is that ONT dRNA004 data exhibit a strand-specific pattern (a barbell-\nlike cluster and clusters at two areas for the default group versus an oval-like cluster for the strand\nshift group). This highlights another advantage of ONT direct RNA sequencing, as it retains strand\ninformation that aligns with the transcript orientation.\nIn conclusion, our methods rely more on the flanking sequences context to differentiate RNA editing\nwith variants.\nFigure R1. UMAP analysis of germline variants versus RNA editing sites in multiple datasets.\nThe latest benchmarks should be included in the manuscript (e.g. with CTAT and\nDeepVariant in table R11).\nRE5: We appreciate the reviewer for the suggestion. We have included the benchmarking result with\nCTAT, LongcallR-nn, and DeepVariant RNA model in the Results section (line 249-252) and\nSupplementary Table 7.\nThe performance shown for DeepVariant is significantly lower than what is reported in their\nGithub case study for HG004. What is the reason for this difference? E.g., if a method is not\ndesigned to be used in a particular regime (e.g. on some variants included in the evaluation),\nperformance could be shown as a function of different regimes.\nRE6: We thank the reviewer for their comments. The results presented in Table R11 are consistent\nwith those reported in the DeepVariant MAS-Seq case study, where the evaluation was restricted to a\nmini-region of chr20 (using the “--target-regions” option in hap.py) with a minimum read coverage\nof 10× and limited to exonic regions (1,002 SNPs and 135 Indels). In contrast, our benchmark\nemployed a minimum coverage of 4× and encompassed broader genomic regions (including 6,689\nSNPs and 1,019 Indels). Since lower-coverage regions inherently present greater challenges for\nvariant detection—as demonstrated in Figure 5a—our evaluation naturally yields a lower F1-score\ncompared to DeepVariant’s targeted analysis. We attribute this performance discrepancy primarily to\nthe inclusion of lower-coverage data in our benchmarking.\nTo further validate our findings, we replicated their analysis using the same mini BED region and\ncalling parameters. As shown in Table R4 below, Clair3-RNA and DeepVariant RNA model achieved\nSNP F1-scores of 95.60% and 95.62% and Indel F1-score 81.57% and 79.97%, respectively,\ndemonstrating performance comparable to their reported results. Only a minor discrepancy was\nfound (2 FPs for SNPs and 1 FP for Indels) due to “FlagCorrection” processing in their BAM file.\nAdditionally, we evaluated both models across chr20 callable regions using different DP and AD\nthresholds. The result showed that Clair3-RNA demonstrated superior performance, outperforming\nthe DeepVariant RNA model by an average 3.26% F1-score in SNP and 5.51% F1-score in Indel.\nVariant Read Allele\nDataset Sample Caller Precision Recall F1-score TRUTH.FP TRUTH.FN TRUTH.TP QUERY.TP\ntype coverage depth\nClair3-RNA with phasing 95.30% 92.07% 93.66% 313 546 6,343 6,346\nDP≥4 AD≥2\nDeepVariant (RNA model) 93.10% 87.22% 90.06% 434 857 5,849 5,852\nClair3-RNA with phasing 97.43% 94.28% 95.83% 110 253 4,168 4,169\nDP≥8 AD≥2\nDeepVariant (RNA model) 95.67% 89.44% 92.45% 170 444 3,759 3,760\nClair3-RNA with phasing 97.81% 94.58% 96.16% 79 202 3,522 3,523\nSNP DP≥10 AD≥2\nDeepVariant (RNA model) 96.45% 89.68% 92.94% 116 362 3,147 3,148\nClair3-RNA with phasing 97.80% 95.64% 96.71% 77 156 3,424 3,425\nDP≥10 AD≥4\nDeepVariant (RNA model) 96.99% 90.91% 93.85% 95 306 3,059 3,060\nPacBio DeepVariant BED with DP≥10 Clair3-RNA with phasing 99.35% 92.12% 95.60% 6 79 923 923\nMAS-Seq and exnomic region DeepVariant (RNA model) 98.93% 92.52% 95.62% 10 75 927 927\nHG004 DP≥4 AD≥2 Clair3-RNA with phasing 78.06% 65.46% 71.21% 188 352 667 669\nchr20 DeepVariant (RNA model) 80.76% 57.95% 67.48% 142 423 583 596\nClair3-RNA with phasing 81.85% 66.67% 73.48% 106 238 476 478\nDP≥8 AD≥2\nDeepVariant (RNA model) 83.23% 57.91% 68.30% 85 298 410 422\nClair3-RNA with phasing 83.47% 65.30% 73.27% 79 211 397 399\nIndel DP≥10 AD≥2\nDeepVariant (RNA model) 83.76% 56.21% 67.27% 69 268 344 356\nClair3-RNA with phasing 86.62% 70.94% 78.00% 61 161 393 395\nDP≥10 AD≥4\nDeepVariant (RNA model) 86.87% 59.82% 70.85% 52 223 332 344\nDeepVariant BED with DP≥10 Clair3-RNA with phasing 86.67% 77.04% 81.57% 16 31 104 104\nand exnomic region DeepVariant (RNA model) 87.93% 73.33% 79.97% 14 36 99 102\nTable R4. Performance comparison in PacBio MAS-Seq HG004 chr20.\nReferences:\n[1] Lin T Y, Goyal P, Girshick R, et al. Focal loss for dense object detection[C]//Proceedings of the\nIEEE international conference on computer vision. 2017: 2980-2988.\n[2] Luo R, Wong C L, Wong Y S, et al. Exploring the limit of using a deep neural network on pileup\ndata for germline variant calling[J]. Nature Machine Intelligence, 2020, 2(4): 220-227.\n[3] Liu Z, Quinones-Valdez G, Fu T, et al. L-GIREMI uncovers RNA editing sites in long-read RNA-\nseq[J]. Genome Biology, 2023, 24(1): 171."}, {"round": 3, "reviewer_comments": {"reviewer_3": "Unfortunately the model evaluation strategy presented in the latest version of the manuscript continues to suffer from data leakage (as described in the previous rounds of review) — since information from the test set (for which key results/figures are presented in the main text) is still partially shared with the training set.\nThe authors have pushed back on correcting this issue stating, in the latest reply, that it is “impractical” to perform the benchmarking on samples unseen during training since (1) they don’t have access to the training datasets/rulesets of other methods and (2) since other papers didn’t avoid data leakage and have been published. While I agree that it is challenging to correctly evaluate other tools (since their training datasets are not available), I think the authors should still evaluate their own method correctly (while noting that other tools’ performance may be biased in the text -- something the reviewers can take into account). Secondly, the fact that several other tools have been published with a flawed evaluation strategy should not serve as the reason to continue this practice; and there is an increasing number of papers that do correctly take data leakage into account. Finally, the authors still have a very large number of samples left after removing training examples from their test set, so they wouldn’t need to retrain their models; they also already have several supplementary tables which could be leveraged as their main results.\nI think it is very important to address this issue to both (1) accurately evaluate the presented model and (2) foster an appropriate level of scientific rigor and reproducibility in the rapidly growing field of ML-based genome analysis."}, "author_response": "Responses to\nComments of reviewers\nClair3-RNA: A deep learning-based small variant caller for long-read RNA sequencing\ndata\nReviewer 3:\nUnfortunately the model evaluation strategy presented in the latest version of the\nmanuscript continues to su6er from data leakage (as described in the previous rounds of\nreview) — since information from the test set (for which key results/figures are presented in\nthe main text) is still partially shared with the training set.\nRE: Please allow us to explain what data leakage is, and why our case is not a data leakage.\nData leakage happens when a model uses information during training that would not be\navailable at the time of prediction. This can lead to overly optimistic performance metrics\nduring validation, as the model has access to data it shouldn't have. To give an easy-to-\nunderstand example, say we trained a model with picture A (a dog, a labrador) as input, and\nprediction ‘dog’ as output, data leakage means in testing, picture A is used again. When\npicture B, which is a corgi that is diJerent from the labrador in picture A, is used for testing,\nit is considered appropriate, and it is not a data leakage.\nReviewer 3’s claim is that “a genome position used in training is not allowed to be used in\ntesting, even if the human sample used for training is totally diJerent from the sample used\nfor testing (in our manuscript, HG002 for training, HG004 and HG005 for testing, HG005 is\nbiologically unrelated to HG002)”. Let us paraphrase Reviewer 3’s claim using the dog\npicture example: “dog used in training is not allowed to be used in testing, even if the input\nfor testing is a corgi, while no corgi but only labrador was used for training”. It is not hard to\nunderstand that the testing input picture of a corgi is totally diJerent from the training input\npicture of a labrador. Similarly, in genetic variant calling, the testing input of sample HG005\nof a genome position, is diJerent from the training input of sample HG002 at the same\ngenome position in every aspect, including read alignments, coverage, base quality, etc.\nThe genome position information itself is not any part of the input to the model for variant\ncalling. Simply put, the input tensor of the same genome position of HG005 (a corgi) and\nHG002 (a labrador) are totally diJerent. The only thing the two diJerent input tensors share\nin common is that they were used for the same prediction task – is it a variant or not? Just\nlike the corgi picture and the labrador picture – is it a dog or not? To stress our point again,\nthere is no experiment in our manuscript that qualifies data leakage.\nAccording to the 1k genome project, any two individuals, even ethnically diJerent, share at\nleast 30% of their genetic variants. With more individuals used in variant calling model\ntraining, which has been one of the primary goals since the beginning of using deep-\nlearning models for variant calling, much fewer genome positions would be found not used\nin training. With more samples added to model training, according to Reviewer 3, the\ntesting dataset needs to be continuously changing. Is a continuously changing testing\ndataset an acceptable practice in machine learning? No. DiJerent teams will use a\ndiJerent set of training samples, according to Reviewer 3, diJerent teams will have\ndiJerent genome positions to exclude, thus using diJerent testing datasets. Is the use of\ndiJerent testing datasets an acceptable practice in machine learning? No. Since the\ngenome has only 3 billion positions and every position has a chance to mutate, with a large\nenough number of samples, all 3 billion can be included in training. Then there will be no\nposition qualified for testing.\nIn cancer genomics, if panel sequencing, amplicon sequencing, or adaptive sequencing is\nused, variant calling (also called genotyping in clinical genomics) is usually restricted to\nonly a few genes of interest, or even only a few variants of interest. For method\ndevelopment, the norm is, all genome positions of these variants of interest will be used for\ntraining, and all genome positions of a diJerent batch of samples will be tested. Then it is\nhard to benchmark any of the published variant calling methods for cancer genomics.\nIn fact, the phrase “data leakage” was not used by the reviewer in all previous rounds of\ncomments, or otherwise we would explain and defend strongly. Data leakage is in machine\nlearning basics, and it is unlikely that a field (genetic variant calling) that is already with\nhundreds of publications in renowned venues and serving critical applications including\nclinical diagnosis and pathogen detection, is still suJering from data leakage. If the\naccusation of ‘data leakage’ is not because of a misunderstanding, we sincerely suggest to\nrevise the variant calling benchmarking standards established by NIST (National Institute of\nStandards and Technology) in the papers “Best Practices for Benchmarking Germline Small\nVariant Calls in Human Genomes” (Nature Biotechnology, 2019) and “PrecisionFDA Truth\nChallenge V2: Calling variants from short and long reads in diJicult-to-map regions” (Cell\nGenomics, 2022). Figure R1, excerpted from (Cell Genomics, 2022), shows the training and\ntesting scheme of PrecisionFDA v2. In this scheme, variants from HG002 were used for\nmodel training, while those from HG003 and HG004 were used for testing. The\nbenchmarking practice proposed by NIST has been broadly used by established variant\ncallers including but not limited to, GATK, FreeBayes, Samtools, Mutect2, Strelka2,\nlongshot, Neusomatic, the Clair series, the DeepVariant series, and NanoCaller. A practice\ncan hardly go wrong with hundreds of papers published in renowned venues, thousands of\ncitations for each of quite a number of these papers, and over millions of accumulated\ndownloads, as shown in public statistics for both the Clair series and DeepVariant series\nvariant callers.\nThe authors have pushed back on correcting this issue stating, in the latest reply, that it is\n“impractical” to perform the benchmarking on samples unseen during training since (1)\nthey don’t have access to the training datasets/rulesets of other methods and (2) since\nother papers didn’t avoid data leakage and have been published. While I agree that it is\nchallenging to correctly evaluate other tools (since their training datasets are not available),\nI think the authors should still evaluate their own method correctly (while noting that other\ntools’ performance may be biased in the text -- something the reviewers can take into\naccount). Secondly, the fact that several other tools have been published with a flawed\nevaluation strategy should not serve as the reason to continue this practice; and there is an\nincreasing number of papers that do correctly take data leakage into account. Finally, the\nauthors still have a very large number of samples left after removing training examples from\ntheir test set, so they wouldn’t need to retrain their models; they also already have several\nsupplementary tables which could be leveraged as their main results.\nRE: “Data leakage” was not a comment in any previous rounds of review. Instead, what\nreviewer 3 asked in the previous round was to explain why there is a drop in recall. We\nbelieve we have explained clearly in detail in our previous reply, and we benchmarked an\nadditional tool called longcallR with the same testing datasets as an external validation for\nour explanation. As we explained, in fact there exists no data leakage, thus there is no need\nto perform the benchmarking only at genome positions not used in model training. Again,\nwe would like to clarify that it is absolutely scientifically correct that a variant from the\nsame sample at the same genome position used in training shall not be used in testing. The\nexperiments in our manuscript fully followed the rule, which is also the norm of the Clair\nseries variant callers, and other published deep-learning-based variant callers. But not\nallowing the same genome position for testing, even if the samples are diJerent and thus\nthe alignments, coverage, base qualities, and everything in the input, is wrong. Corgi is not\na labrador, and the purpose of machine learning is to allow a model trained on only\nlabradors to also identify corgis as dogs.\nWe want to point out that we have not denied Reviewer 3’s request to benchmark only at\ngenome positions not used in training. In fact, in the previous revision, 1) we have added\nwhat reviewer 3 has asked in the Results section under “Performance analysis after\nremoving shared variants between samples”, 2) Table R3 included “HG002-unseen”\n(benchmarking HG004 and HG005 only at genome positions having no variant in HG002 for\nmodel training) results for PacBio and ONT datasets on both Clair3 and LongcallR, and 3)\nSupplementary Table 12 was added with both the norm and “HG002-unseen” results for\nClair3-RNA. The modifications have been highlighted in the manuscript.\nWe still believe Reviewer 3’s claim that our method has a ‘data leakage’ problem is because\nof misunderstanding. We sincerely apologize for any confusion we have caused.\nFigure R1. Training and testing scheme of PrecisionFDA Truth Challenge V2.\nThe paper states that “Participants were provided sequencing reads (FASTQ files) from Illumina, PacBio HiFi,\nand ONT for the GIAB Ashkenazi trio (HG002, HG003, and HG004). Participants uploaded VCF files for each\nindividual before the end of the challenge, and then the new benchmarks for HG003 and HG004 were made\npublic.”\nSource from the paper “PrecisionFDA Truth Challenge V2: Calling variants from short and long reads in\ndiTicult-to-map regions” in Cell Genomics, 2022."}, {"round": 4, "reviewer_comments": {"reviewer_1": "It is important to demonstrate that the method is not achieving better results through overfitting. The experiments that the authors have performed (chromosome 20 only benchmarking in Supp. Table 5 and excluding shared candidates in Supp Table 12) are sufficient to indicate that the method is not overfitting the training data.\nOne comment about the new analysis reported on lines 496-509: If I understood correctly, the authors removed both true positives (variants) and false positives (artifacts) from benchmarking. If so, the title of the section should be revised to reflect this.\nAlso, a preprint for the longcallR method that the authors have used for comparison/benchmarking has come out earlier this year. It would be good to cite it (https://www.biorxiv.org/content/10.1101/2025.05.26.656191v1) rather than just the github link.", "reviewer_2": "Data leakage is a valid concern for machine learning based approach. Theoretically, the model of Clair3-RNA is not affected from data leakage. The pileup features depend on the sequencing read coverage, and this varies from sample to sample especially for RNA-seq data, where the gene/transcript expression can be substantially different. In other words, for the same position, the training sample and test samples would have different tensors in most cases. The authors did a good job to further address the concern by holding out chr20 during training (in the first version of the manuscript). The other reviewers raised a good question on the performance on chr20 in the test data, and the authors have shown the accuracy is comparable with whole genome data and outperformed competing methods in the revision. The authors also explicitly evaluated the accuracy on the variants that do not overlap between HG002 and the test samples, and showed comparable accuracy. The authors tried to explain the accuracy difference, but the argument might be inaccurate. In particular, the authors wrote “…which is expected due to an average ~43% reduction in benchmarked true positives without a substantial proportional change in false negatives” may not be the reason for the sensitivity decrease. For example, true positive and false negative sum up to the all true (variant), so if the true positive number are reduced more than expected after removing overlapped variant from the data set, then the false negative should be reduced less than expected. Therefore, the term “proportional change in false negative” is somewhat vague and difficult to follow.\nFrom my understanding, suppose the X denote the overlapped variant between HG002 and the test sample, and the recall for that part is r(X). The new evaluated recall for the non-overlapped variant should be (TP-r(X)*X)/(T-X), which will be smaller than the original recall TP/T as long as r(X) as is greater than TP/P. In other words, the decrease of recall implies that the recall for the overlapped variants (X) are higher than other variants. By plugging the number based on Table S12a, for PacBio BCM HG004, T=163024, TP=154482 (inferred from overall recall 94.76%), X=82761, the recall r(X)=96.72% based on the recall on the non-overlapped variant is 92.74%. Nevertheless, the precision and recall for P-X set is still very high, suggesting the data leakage may not be a major issue. The explanation of losing recall/F-1 might not be just due to numeric reasons like the changing of FP and FN, it could be biological reasons. One potential reason could be that the shared variant are common variants and they might be far from the exon-intron boundaries so they don’t affect the alternative splicing much and can be a retained as a common variants in the population. In sum, the authors may look into the reason for the recall change a bit more (or I made a mistake in the derivation) if explanations are mentioned in the manuscript.", "reviewer_4": "I have reviewed the original manuscript and the author's response to reviewer 3. I did not see the full reviewer's comments though, since I guess the authors only cited part of reviewer 3's full comments in their response.\nNevertheless, any two genomes will share a substantial proportion of genetic variants, so it is unavoidable that a genetic variant used in training will show up in testing again. I believe the reviewer is probably concerned with the fact that certain features of the variants used in training may artificially make the prediction performance higher if the same sets of variants are also in testing data on a different sample. This is a reasonable concern, in fact, many variant callers evaluate performed on specific subsets of variants (such as common variants vs rare variants, as rare variants are unlikely to be in training set). Some others also perform training on chr1-20 and testing on chr21-22.\nFrom the author's response, it seems that the reviewer has already requested to benchmark only at genome positions not used in training in the first round of review. This is reasonable, and I think the authors have already addressed this concern (“Performance analysis after removing shared variants between samples”).\nIn summary, I believe that the authors provided reasonable responses to address the reviewer's comments on the possibility that some information from training is leaked to testing."}, "author_response": "Responses to\nComments of reviewers\nClair3-RNA: A deep learning-based small variant caller for long-read RNA sequencing\ndata\nNCOMMS-24-75685D\nReviewer: #1\n----------------------------------------------------------\nIt is important to demonstrate that the method is not achieving better results through\noverfitting. The experiments that the authors have performed (chromosome 20 only\nbenchmarking in Supp. Table 5 and excluding shared candidates in Supp Table 12) are\nsufficient to indicate that the method is not overfitting the training data.\nOne comment about the new analysis reported on lines 496-509: If I understood correctly,\nthe authors removed both true positives (variants) and false positives (artifacts) from\nbenchmarking. If so, the title of the section should be revised to reflect this.\nAlso, a preprint for the longcallR method that the authors have used for\ncomparison/benchmarking has come out earlier this year. It would be good to cite it\n(https://www.biorxiv.org/content/10.1101/2025.05.26.656191v1) rather than just the github\nlink.\nRE: We sincerely thank the reviewer for the valuable comments and suggestions. We have\nupdated the citations for both longcallR and longcallR-nn to reference the newly available\npreprint (https://www.biorxiv.org/content/10.1101/2025.05.26.656191v1\n) instead of the GitHub link.\nRegarding the title of the section, we appreciate the reviewer’s clarification and have revised\nit into “Performance analysis after removing shared variants and artifacts between\nsamples” to accurately reflect that both true positives (variants) and false positives\n(artifacts) were excluded from benchmarking in this analysis.\nReviewer: #2\n----------------------------------------------------------\nData leakage is a valid concern for machine learning based approach. Theoretically, the\nmodel of Clair3-RNA is not affected from data leakage. The pileup features depend on the\nsequencing read coverage, and this varies from sample to sample especially for RNA-seq\ndata, where the gene/transcript expression can be substantially different. In other words, for\nthe same position, the training sample and test samples would have different tensors in most\ncases. The authors did a good job to further address the concern by holding out chr20 during\ntraining (in the first version of the manuscript). The other reviewers raised a good question\non the performance on chr20 in the test data, and the authors have shown the accuracy is\ncomparable with whole genome data and outperformed competing methods in the revision.\nThe authors also explicitly evaluated the accuracy on the variants that do not overlap\nbetween HG002 and the test samples, and showed comparable accuracy. The authors tried\nto explain the accuracy difference, but the argument might be inaccurate. In particular, the\nauthors wrote “…which is expected due to an average ~43% reduction in benchmarked true\npositives without a substantial proportional change in false negatives” may not be the reason\nfor the sensitivity decrease. For example, true positive and false negative sum up to the all\ntrue (variant), so if the true positive number are reduced more than expected after removing\noverlapped variant from the data set, then the false negative should be reduced less than\nexpected. Therefore, the term “proportional change in false negative” is somewhat vague\nand difficult to follow.\nFrom my understanding, suppose the X denote the overlapped variant between HG002 and\nthe test sample, and the recall for that part is r(X). The new evaluated recall for the\nnon-overlapped variant should be (TP-r(X)*X)/(T-X), which will be smaller than the original\nrecall TP/T as long as r(X) as is greater than TP/P. In other words, the decrease of recall\nimplies that the recall for the overlapped variants (X) are higher than other variants. By\nplugging the number based on Table S12a, for PacBio BCM HG004, T=163024, TP=154482\n(inferred from overall recall 94.76%), X=82761, the recall r(X)=96.72% based on the recall\non the non-overlapped variant is 92.74%. Nevertheless, the precision and recall for P-X set\nis still very high, suggesting the data leakage may not be a major issue. The explanation of\nlosing recall/F-1 might not be just due to numeric reasons like the changing of FP and FN, it\ncould be biological reasons. One potential reason could be that the shared variant are\ncommon variants and they might be far from the exon-intron boundaries so they don’t affect\nthe alternative splicing much and can be a retained as a common variants in the population.\nIn sum, the authors may look into the reason for the recall change a bit more (or I made a\nmistake in the derivation) if explanations are mentioned in the manuscript.\nRE: We sincerely thank the reviewer for the detailed and thoughtful analysis regarding the\nobserved decrease in recall. We appreciate the reviewer’s clarification and agree that the\nprevious explanation could be made clearer. In the revised manuscript, we have refined the\ndescription to better reflect the underlying reasoning:\n“The recall decreased by an average of 4.14% and 2.72% before and after phasing,\nrespectively. This reduction is largely attributable to the exclusion of variants overlapping\nwith the HG002 benchmark, which resulted in an average ~43% decrease in the number of\nbenchmarked true positives. These overlapping variants are likely enriched for\nwell-supported, common variants that are easier to detect; therefore, their removal naturally\nleads to a modest reduction in recall. Importantly, the precision and recall for the\nnon-overlapping variants remain high, indicating that Clair3-RNA did not overfit to specific\nsamples.”\nReviewer: #4\n----------------------------------------------------------\nI have reviewed the original manuscript and the author’s response to reviewer 3. I did not\nsee the full reviewer’s comments though, since I guess the authors only cited part of\nreviewer 3's full comments in their response.\nNevertheless, any two genomes will share a substantial proportion of genetic variants, so it\nis unavoidable that a genetic variant used in training will show up in testing again. I believe\nthe reviewer is probably concerned with the fact that certain features of the variants used in\ntraining may artificially make the prediction performance higher if the same sets of variants\nare also in testing data on a different sample. This is a reasonable concern, in fact, many\nvariant callers evaluate performed on specific subsets of variants (such as common variants\nvs rare variants, as rare variants are unlikely to be in training set). Some others also perform\ntraining on chr1-20 and testing on chr21-22.\nFrom the author’s response, it seems that the reviewer has already requested to benchmark\nonly at genome positions not used in training in the first round of review. This is reasonable,\nand I think the authors have already addressed this concern (“Performance analysis after\nremoving shared variants between samples”).\nIn summary, I believe that the authors provided reasonable responses to address the\nreviewer’s comments on the possibility that some information from training is leaked to\ntesting.\nRE:We sincerely thank the reviewer for the thoughtful comments and for recognizing our\ne our\nefforts to address potential data leakage concerns."}]

Decision letter