Terminal modifications independent cell-free RNA sequencing enables sensitive early cancer detection and classification

Domain: earth_science
Published: 2024-01-02
Editorial outcome: minor_revision

Abstract
Cell-free RNAs (cfRNAs) offer an opportunity to detect diseases from a transcriptomic perspective, however, existing techniques have fallen short in generating a comprehensive cell-free transcriptome profile. We develop a sensitive library preparation method that is robust down to 100 µl input plasma to analyze cfRNAs independent of their 5’-end modifications. We show that it outperforms adapter ligation-based method in detecting a greater number of cfRNA species. We perform transcriptome-wide characterizations in 165 lung cancer, 30 breast cancer, 37 colorectal cancer, 55 gastric cancer, 15 liver cancer, and 133 cancer-free participants and demonstrate its ability to identify transcriptomic changes occurring in early-stage tumors. We also leverage machine learning analyses on the differentially expressed cfRNA signatures and reveal their robust performance in cancer detection and classification. Our work sets the stage for in-depth study of the cfRNA repertoire and highlights the value of cfRNAs as cancer biomarkers in clinical applications.

Introduction
Early diagnosis of malignancy offers the greatest opportunity to improve long-term patient survival by increasing the likelihood of cure1. Utilization of liquid biopsies for cancer detection is rapidly gaining prominence in clinical practice2. Although cell-free DNA-based assays interrogating somatic mutations3, fragmentation patterns4,5, and methylation signatures6–8 have shown high specificity in detecting and localizing cancer, complementing them with transcriptomic information will improve their sensitivity in diagnosing cancer at its earliest stage. Cell-free RNAs (cfRNAs) represent a landscape of extracellular transcripts and thus provide extensive clinical value9,10. Alterations in cellular RNA expression are a dynamic process that can serve as an indicator of disease11. Overexpression of specific tumor-associated transcripts can result in amplification of tumor-derived cfRNA signals in the blood12. The release of cfRNAs into the bloodstream is not solely through cell death but can also occur through dynamic signaling processes of extracellular vesicles13, expanding the opportunities to detect cancer in individuals with limited levels of circulating tumor DNA.

Over the past decades, the characterization of cfRNAs has mostly focused on microRNAs (miRNAs)14–16, which have limited tissue or disease specificity17. Recently, there has been a rising interest in the study of messenger RNAs (mRNAs), which have been shown to possess tissue and cancer subtype specificity12 as well as cell types of origin specificity18. mRNAs can also distinguish pre-malignant conditions of cancer19 and evaluate pathophysiological alterations in the brain20. Moreover, mRNA signatures have demonstrated value in obstetric care through early prediction of preeclampsia21–24. However, miRNA and mRNA account for only a small fraction of the cfRNA repertoire. The potential value of other kinds of cfRNAs in plasma remains elusive.

There are two major barriers to comprehensive cfRNA profiling. First, large volumes of plasma (up to 2 ml) are required because the recovery rate of cfRNAs is typically low. Second, cfRNAs are severely fragmented, resulting in diverse terminal modifications25. However, conventional adapter ligation-based library preparation workflows are only sensitive to RNA fragments with a 5’-phosphate (5’-P) and 3’-hydroxyl (3’-OH) group26,27. Therefore, innovative cfRNA profiling technologies that can efficiently handle low abundance and highly fragmented samples are urgently needed.

Here, we develop Splint Ligation and Phosphate-independent RNA Sequencing (SLiPiR-seq), a sensitive library preparation method that remains reliable with as little as 100 µl of input plasma and allows profiling of cfRNAs independent of their 5’-end phosphorylation, coupled with a bioinformatics analysis pipeline (Supplementary Fig. 1) for the discovery and absolute quantification of non-canonical small RNA species in plasma. To verify the clinical feasibility of SLiPiR-seq, we conduct transcriptome-wide characterizations in patients with various cancer types and cancer-free donors, demonstrating its ability to detect transcriptomic changes occurring in early-stage tumors. We further carry out machine learning analyses and establish multiple models to detect and localize different types of cancers, highlighting the strong potential of cfRNA signatures for future clinical applications.

Results
Technology optimizations of SLiPiR-seq
To perform library preparation for highly fragmented cfRNAs, we optimized the S-Poly(T) Plus method that was previously developed for high-throughput small RNA profiling28–30. In brief, the input cfRNA was subjected to 3’-end polyadenylations along with reverse transcription (RT) in a one-step reaction (Fig. 1a). A custom synthetic RT primer consisting of oligo(dT), sample barcode, and sequencing adapter sequences was designed. cDNA with a 5’-adapter was obtained after the polyadenylations/RT reaction. The excess primer was then depleted using exonuclease I (ExoI). This is a critical step as any remaining RT primer can lead to adapter-primer ligation (~160 nt) (Fig. 1b), which would significantly impact the ratio of informative reads and the number of RNAs detected (Fig. 1c, Supplementary Fig. 2a–c). Consistently, similar results were observed upon the addition of excessive RT primer (Fig. 1d, Supplementary Fig. 2d, e). Subsequently, a double-stranded adapter with 3’-overhanding degenerate nucleotides was annealed to the 3’ end of cDNA by splint ligation (Fig. 1a). The randomized bases on the degenerate extension can transiently hybridize to cDNA31,32. The blocking strand of the adapter was then removed by excision of uracil residues using USER enzyme, followed by library amplification and sequencing. Removal of the blocking strand can significantly improve library amplification efficiency by facilitating primer binding (Fig. 1e, Supplementary Fig. 3a–d). To reduce RNA/DNA loss during library preparation, we further optimized SLiPiR-seq protocol by making the entire workflow compatible within a single tube.

To generate stable and reproducible transcriptome profiles, we extensively benchmarked SLiPiR-seq regarding plasma input volumes (from 12.5 µl to 400 µl) and the impact of various pre-analytical conditions (i.e., blood standing time before plasma separation and number of freeze/thaw cycles of the plasma). We found that the lower limit of plasma input for reliable SLiPiR-seq results is 100 µl (Fig. 1f, Supplementary Fig. 3e, f), significantly lower than the minimum requirement for most commercially available small RNA library preparation kits. The ratio of unspecified reads increased as the starting plasma volume decreased (Fig. 1g). We observed no significant difference in the recovery of total cfRNAs among plasma samples separated after 3, 6, and 9 hours of blood collection (3h vs. 6h p=0.291; 3h vs. 9h p=0.189) (Supplementary Fig. 4a). Furthermore, our results showed that cfRNAs remain stable upon one freeze/thaw cycles of plasma (p=0.237) (Supplementary Fig. 4b). However, the informative reads were significantly reduced upon more freeze/thaw cycles (p<0.001). Correlation analysis suggested that long noncoding RNAs (lncRNAs), mRNAs, miRNAs, small nuclear RNAs (snRNAs), small nucleolar RNAs (snoRNAs), and transfer RNA-derived small RNAs (tsRNAs) are stable in plasma under different pre-analytical conditions (Supplementary Fig. 4c, d). Taken together, we determined to separate the plasma within 3 hours of blood collection and to preserve 200 µl of aliquoted plasma by freezing it until SLiPiR-seq was performed to ensure one freeze/thaw cycle.

Technology assessments of SLiPiR-seq
To assess whether the sequencing results of SLiPiR-seq can accurately represent the molecular characteristics of the plasma transcriptome, we compared it to quantitative real-time PCR (qPCR), the gold standard technique for RNA quantification. We randomly selected varying expression levels of cfRNAs detected by SLiPiR-seq for comparison. Primer design was based on our previously described S-Poly(T) Plus strategy28,29. Notably, we observed a high concordance (R=0.86) between read counts profiled by SLiPiR-seq and Ct values measured by qPCR (Fig. 2a, Supplementary Fig. 5, Supplementary Data 1), suggesting that SLiPiR-seq can accurately profile the plasma transcriptome.

Taking advantage of the splint ligation, we hypothesized that SLiPiR-seq could be used to profile cfRNAs independently of their 5’-end phosphorylation, thereby providing a broader coverage of the cfRNA repertoire. To test this hypothesis, we compared SLiPiR-seq with an adapter ligation-based method (NEBNext Small RNA Library Prep Set, NEBNext) using synthetic small RNAs with or without a 5’-P group. Our findings revealed that SLiPiR-seq can reliably construct a sequencing library regardless of the presence or absence of a 5’-P group. However, adapter ligation of NEBNext failed when using RNAs lacking a 5’-P group (Fig. 2b). Consistently, sequencing read counts generated by SLiPiR-seq showed no difference between the two groups given the same amount of library input, whereas NEBNext exhibited sensitivity only when the 5’-end phosphorylation was present (Fig. 2c, Supplementary Data 2).

We subsequently compared the two library preparation methods to evaluate their capacity to detect cfRNAs in human plasma samples. SLiPiR-seq detected 3.68 times more total RNA species than NEBNext (65204 vs. 17696) with 20 million clean reads (Fig. 2d). Specifically, 17932 mRNA and 12236 lncRNA species were detected by SLiPiR-seq, which were 4.37 and 14.34 times more than NEBNext, respectively. More importantly, by performing a cross-analysis with tissue-specific genes of various tissues from the Human Protein Atlas database, we discovered that a considerable number of mRNAs detected by SLiPiR-seq showed tissue-specific elevated expression patterns (Fig. 2e, Supplementary Data 2). Compared to NEBNext, SLiPiR-seq displayed a substantially increased proportion of mRNAs (7.9% vs. 1.9%) and lncRNAs (8.4% vs. 0.7%) (Fig. 2f). In addition, the fragment size of mRNAs (NEBNext: α=3.32, β=0.10; SLiPiR-seq: α=10.22, β=0.24, fitted to Gamma distribution) and lncRNAs (NEBNext: α=6.29, β=0.21; SLiPiR-seq: α=12.07, β=0.35) exhibited smoother distribution and broader coverage in SLiPiR-seq, as suggested from the kurtosis values (mRNA: 13.22 for NEBNext vs. 11.25 for SLiPiR-seq; lncRNA: 30.58 for NEBNext vs. 14.92 for SLiPiR-seq) (Fig. 2g). This finding was further validated by analyzing the fragment length of the twenty most highly expressed mRNA and lncRNA species (Supplementary Fig. 6a, b). Exonic mRNA and lncRNA sequences detected by both methods were derived from the sense orientation rather than the antisense orientation (Supplementary Fig. 6c), suggesting successful recovery of transcript. Our findings also indicated that SLiPiR-seq detected a broader coverage of the 3’ end of mRNAs than NEBNext (Supplementary Fig. 6d, e).

miRNAs emerged as the predominant RNA type (60.7%) in NEBNext results (Fig. 2f), consistent with the fact that plasma miRNAs are in the canonical form with 5’-P and 3’-OH termini33. Although miRNAs accounted for a lower proportion (10.3%) of total cfRNAs in SLiPiR-seq, the expression levels of miRNAs detected by both methods were strongly correlated (R=0.938, Supplementary Fig. 7). We also observed low correlations for the expression levels of other RNA types between the two methods, because SLiPiR-seq identified a greater number of RNA species than NEBNext. Together, these observations suggest that SLiPiR-seq holds great promise to uncover the biological implications of diseases in plasma, given its enhanced detection of more diverse and integral RNA species associated with gene expression and regulation.

Characterizations of tsRNAs, rsRNAs, and ysRNAs
Intrigued by the abundant expression of tsRNAs (14.4%), ribosomal RNA-derived small RNAs (rsRNAs, 6.0%), and Y RNA-derived small RNAs (ysRNAs, 9.6%) revealed by SLiPiR-seq (Fig. 2f) and recognizing the absence of publicly available reference genome for rsRNA and ysRNA, we purposed to establish a reference genome for cfRNAs aligned to rRNAs and YRNAs. We mapped the small RNAs to the parent rRNAs (28S, 18S, 5.8S, and 5S rRNA) and YRNAs (RNY1, RNY3, RNY4, and RNY5) and extracted unique sequences from the outputs of alignment to ensure a complete annotation of each rsRNA and ysRNA (Fig. 3a). From all the samples included in this study, we detected a total of 45397 rsRNA and 2664 ysRNA unique sequences (Supplementary Data 3). Together with 28,824 unique tsRNA sequences acquired from an online tsRNA database34, we used a custom Python script to call the read counts of tsRNAs, rsRNAs, and ysRNAs. This script only tallies a fragment if it aligns perfectly with a sequence in the FASTA file (Fig. 3a).

Next, we analyzed the differences between SLiPiR-seq and NEBNext in characterizing tsRNAs, rsRNAs, and ysRNAs. Compared to NEBNext, the number of amino acids corresponding to the parent tRNAs of tsRNAs was significantly increased in SLiPiR-seq (Fig. 3b). Furthermore, our results showed that SLiPiR-seq identified a higher overall sequencing depth and broader coverage of the 3’ end than NEBNext, because RNA fragments resulting from 3’ end cleavage are more likely to lack a 5’-P group. For example, SLiPiR-seq displayed increased read counts for tsRNAs starting at the 3’ end and ending at the anticodon position (3’-tRHs) or ending before the anticodon position (3’-tRFs) (Fig. 3c, d)35. This finding was particularly evident for the 3’ end of 5.8S rRNA, RNY4, and RNY5, which were undetectable using NEBNext (Fig. 3e, f). These findings highlight the potential of SLiPiR-seq in discovering underexplored small RNAs in blood.

Identification of cfRNA signatures for lung cancer using SLiPiR-seq
To investigate whether SLiPiR-seq can detect differentially expressed (DE) cfRNAs in cancer, we retrospectively profiled a discovery cohort containing plasma samples from 139 lung cancer (LC) patients (cases) and 106 cancer-free individuals (controls). Transcriptome-wide comparison of normalized sequencing reads between the two cohorts showed a high Pearson’s correlation (R=0.987, p<2.2×10−16), indicating that SLiPiR-seq has robust reproducibility (Fig. 4a, Supplementary Fig. 8). We identified 17622 DE cfRNAs (11550 gain and 6072 loss, adjusted p value <0.1) between cases and controls based on a negative binomial model of fragment counts (Fig. 4b, Supplementary Fig. 9, Supplementary Data 4). The majority of DE cfRNAs have increased expression in cancer patients, consistent with the fact that more cfRNAs are released into the circulation in the context of cancer. The three most abundant categories of DE cfRNAs are rsRNAs (73.4%), tsRNAs (8.3%), and piwi-interacting RNAs (piRNAs) (8.2%) (Fig. 4c).

To enhance computational efficiency, we reduced our dataset by filtering out non-DE cfRNAs, followed by removal of cfRNAs with low expression levels (mean read counts <10) and minimal differences between cases and controls (log fold change <0.8). Next, we compared the sum of reads of the remaining DE cfRNAs between early-stage LC patients (stage I and II, N=81), late-stage LC patients (stage III and IV, N=57), and controls. We observed a significant increase in the cumulative expression level of DE cfRNAs in early-stage LC patients compared to controls (p=7.10×10−13), and no difference (p=0.321) between early and late-stage LC patients (Fig. 4d, Supplementary Fig. 10). Hierarchical cluster analysis using the top fifteen DE cfRNAs of each type revealed a higher level of DE cfRNAs in cancer patients than in cancer-free controls (Fig. 4e). As expected, no evident differences in the expression profiles of these cfRNAs between patients with different cancer stages were observed in the heatmap. These results indicate that SLiPiR-seq can detect aberrant transcriptional events that occur early in tumor pathogenesis, enabling SLiPiR-seq as a promising tool for early cancer screening.

Construction of lung cancer detection models
Given that SLiPiR-seq successfully identified tumor-derived DE cfRNAs, we then investigated whether the utility of plasma transcriptome profiles could detect LC. We leveraged three different strategies to screen for candidate features: (1) filter method-“TopN”; (2) wrapper method-“Boruta”; (3) intrinsic method-“LASSO” (see Methods for details). To explore the potential of different kinds of RNAs in cancer detection, we carried out a series of machine-learning analyses. Three different algorithms (logistic regression (LR), random forest (RF), and support vector machine (SVM)) were applied to ensure the reliability of selected candidate features. We randomly partitioned the discovery cohort into 80% as a training set to train a diagnostic model and 20% as a test set to evaluate model performance. The sample partitioning and train-test processes were repeated 100 times for all feature selection and machine-learning analyses to avoid sampling bias. Consequently, we obtained 27 sets of candidate features, all of which can precisely detect LC in the training and test sets of the discovery cohort (Fig. 5a, Supplementary Table 1a, b). Since the candidate features identified by the LASSO method achieved the highest overall area under the receiver operating characteristic curve (AUC) in the test sets, we used the models trained with these features for subsequent validations.

To verify the performance of these models in detecting early-stage LC, we profiled an independent validation cohort of 26 stage I LC patients and 27 cancer-free individuals using SLiPiR-seq (Fig. 5a, Supplementary Table 1c). Among all candidate cfRNAs, miRNA features detected LC with the highest accuracy (LR AUC=0.905 [interquartile range (IQR) 0.895–0.912]) (Fig. 5b). mRNA features also showed robust LC detection performance (LR AUC=0.846 [0.823–0.860]). Interestingly, the seldom studied snRNA (SVM AUC=0.903 [0.893–0.911]), snoRNA (LR AUC=0.788 [0.772–0.798]), and tsRNA (SVM AUC=0.741 [0.721–0.765]) features revealed considerably great ability in detecting LC as well (Fig. 5b). LR and SVM models training with rsRNA (SVM AUC=0.819 [0.786–0.843]) or ysRNA (SVM AUC=0.793 [0.747–0.829]) features showed high AUC in LC detection. However, models trained by the RF algorithm shared no agreement with these results (Supplementary Table 1c). Although the selected lncRNA and piRNA candidates demonstrated high AUC in the discovery cohort, none of them revealed LC detection capability in the validation cohort. Fifteen representative candidate cfRNAs showed concordant expression levels in both the discovery and validation cohorts (Supplementary Fig. 11). These observations indicate that SLiPiR-seq is promising in identifying not only commonly studied (miRNAs and mRNAs) but also rarely reported (snoRNAs, snRNAs, and tsRNAs) cfRNAs as cancer biomarkers.

Combination of multiple types of cfRNAs achieved higher accuracy
Intrigued by the fact that cancer is a highly heterogeneous disease and that SLiPiR-seq has detected multifarious types of DE cfRNAs, we wondered if combining different RNA categories might improve model accuracy. To this end, we tested all possible combinations of the five types of candidate cfRNAs (aka. 29 mRNAs, 26 miRNAs, 10 snRNAs, 19 snoRNAs, and 38 tsRNAs, Supplementary Data 5) that have successfully undergone the ordeal of all three machine learning algorithms, resulting a total of 31 different combinations. We found that the combination of multiple RNA types provided more accurate results than using a single RNA type in most cases (Supplementary Fig. 12, Supplementary Data 6), suggesting that the combination strategy is beneficial for model performance.

To identify the best combination, we focused on the LR algorithm for calculating AUC, model sensitivity and specificity, and cancer risk scores in subsequent analyses. The top three most accurate combinations ranked by the median of AUC in the test set and validation cohort were “mi+sn+sno” (AUC=0.979), “m+mi+sn+sno” (AUC=0.970) and “mi+sn+sno+ts” (AUC=0.970) (Fig. 5c). Next, we calculated the “risk scores” (aka. model-predicted probabilities) for all study subjects using the top panel “mi+sn+sno”. Of the 137 individuals in the discovery cohort who were assigned a high median risk score (≥0.5), 134 were true LC patients in the real world (97.81% specificity) (Supplementary Fig. 13a). Meanwhile, 134 of 139 LC patients in the discovery cohort were correctly diagnosed by this panel (96.40% sensitivity). As expected, we observed no evident differences in the model-predicted AUC and risk scores between early- and late-stage cancer patients (Fig. 5d, Supplementary Fig. 13b, c). In the validation cohort, all LC patients were diagnosed as high-risk individuals (100% sensitivity). However, the specificity of this panel is low in the validation cohort, 7 out of 33 individuals diagnosed as high-risk were not LC patients in the real world (78.79% specificity) (Supplementary Fig. 13d).

In the context of cancer screening, specificity often assumes a heightened importance than sensitivity. By ranking the specificities of all combinations (Supplementary Data 6), we identified an optimal panel for lung cancer screening in the general population. The “m+sn+sno+ts” panel shows 100% specificity and 99.28% sensitivity in the discovery cohort, and 95.24% specificity and 76.92% sensitivity in the validation cohort (Fig. 5e, Supplementary Fig. 13e). Prioritizing specificity can mitigate the multifaceted consequences of false-positive results, ensuring a more efficient and patient-centered approach to healthcare. These findings demonstrate the robust performance of cfRNA signatures in early-stage LC detection achieved by the integration of multiple types of cfRNAs.

Establishment of cancer classification models
To ascertain the clinical feasibility of SLiPiR-seq, we further explored its ability to identify cfRNA signatures for cancer classification. We expanded our study to encompass a pan-cancer patient cohort, comprising 30 breast cancer (BRC) patients, 37 colorectal cancer (CRC) patients, 55 gastric cancer (GC) patients, and 15 liver cancer (HCC) patients. The high Pearson correlation coefficients of the transcriptome profiles from distinct subject cohorts demonstrate that SLiPiR-seq can produce stable results under varying temporal and spatial conditions (Supplementary Table 2). We performed differential expression analyses for each cancer type in a one-class versus-other-classes paradigm (Supplementary Data 4), followed by LASSO LR to select discriminative features from the five RNA categories of interest (aka. mRNA, miRNA, snRNA, snoRNA, and tsRNA). Consequently, we identified five cancer type-specific panels (aka. 21-cfRNA BRC panel, 33-cfRNA CRC panel, 36-cfRNA GC panel, 33-cfRNA HCC panel, and 30-cfRNA LC panel) (Fig. 6a, Supplementary Fig. 14, Supplementary Data 5). Visualization using t-distributed stochastic neighbor embedding (t-SNE) revealed clear separation based on cancer type (Fig. 6b).

We also combined all cancer patients as a group named “cancers” to compare with cancer-free donors and identified 65 cfRNAs that were commonly upregulated in all five cancer types studied. Subsequently, we conducted one-class versus-other-classes machine learning analyses on the six sets of candidates and observed high AUC values for the held-out test set samples across all classes (Fig. 6c, Supplementary Fig. 15).

We then assessed cancer risk scores as well as five cancer type-specific risk scores for all subjects. Notably, cancer-free individuals (NOR) exhibited low median risk scores in all six panels (Fig. 6d). All cancer patients exhibited high cancer risk scores in the 65-cfRNA common cancer panel. Furthermore, patients with a particular cancer type had high-risk scores only in the corresponding cancer type-specific panel but had low risk scores in the non-corresponding panels (Fig. 6d, Supplementary Fig. 16, Supplementary Data 7). The aforementioned independent validation cohort (N=53) was used to verify the cancer classification performance of these panels. Remarkably, LC patients exhibited high median cancer risk scores (cancer 65 cfRNAs panel) and LC risk scores (LC 30 cfRNAs panel) while maintaining low risk scores in the BRC, CRC, GC, and HCC models (Fig. 6e). Cancer-free individuals had low median risk scores in all six models (Fig. 6e). Of the 27 cancer-free individuals, seven were erroneously assigned high cancer risk scores. Interestingly, only one was misclassified as GC but the remaining six cancer-free individuals had low risk scores in all five cancer type-specific models (Supplementary Data 7). This observation suggests that a two-step verification process could enhance the specificity of the cancer detection assay. A diagnosis should only be made when both the common cancer risk scores and at least one of the cancer type-specific risk scores are high. However, impaired sensitivity is expected as such a strategy will misclassify some real LC patients with only one positive risk scores. Collectively, these results demonstrate that cfRNA signatures identified by SLiPiR-seq are promising in the classification of different cancer types, motivating us to validate the clinical utility of SLiPiR-seq in larger sample cohorts in future investigations.

Methods
Sample acquisition and clinical cohorts
Pooled plasma samples of healthy individuals were used in the technology optimization part. In the clinical feasibility test part, study subjects were recruited from People’s Hospital of Bao’an Shenzhen (N=111 cancer-free individuals, NOR_SZBA), The Fifth People’s Hospital of Suzhou (N=27 cancer-free individuals, NOR_SZDW), The Second People’s Hospital of Shenzhen (N=140 lung cancer patients, LC_SZDE), Peking University Shenzhen Hospital (N=28 lung cancer patients, LC_SZBU), General Hospital of Ningxia Medical University (N=30 breast cancer patients, BRC_NXYK; N=37 colorectal cancer patients, CRC_NXYK; N=55 gastric cancer patients, GC_NXYK; and N=16 liver cancer patients, HCC_NXYK). Inclusion criteria for cancer patients were: (1) patient was diagnosed with cancer through histopathology; (2) patient has not undergone tumor resection surgery or any other form of therapy. The clinical characteristics of all study subjects are available in Supplementary Table 3. The study protocol was reviewed and approved by the Institutional Review Board of all participating sites. Informed written consent was obtained from each participant prior to sample collection.

Whole blood processing
For all cohorts, 2 ml of whole blood samples were collected in EDTA anticoagulated vacutainers. Blood was drawn in the morning from participants who had fasted overnight. Plasma separation was done within 3 hours of whole blood collection. Blood samples were first centrifuged at 1600×g for 10 min at 4°C. The plasma layer was transferred to a clean tube and further centrifuged at 16,000×g for 10 min at 4°C to ensure complete depletion of cell debris. The double-spun plasma was split into 200 µl aliquots and stored at −80°C until RNA isolation.

cfRNA extraction
Frozen plasma samples were thawed on ice prior to cfRNA extraction. 200 µl of plasma samples were subjected to cfRNA extraction using the Apostle MiniMaxTM High-Efficiency cfRNA Isolation Kit (Apostle), following the manufacturer’s protocol with minor modifications. The experiments were conducted through manual operation on the bench up until the step involving the addition of 6.5 µl of Binding Enhancer to the 200 µl supernatant. The mixture was transferred to a 96-well plate with the prepared binding/nanoparticle solution, cfRNA Wash Solution, and 80% Ethanol (two columns) on different columns (eight samples in parallel). The subsequent steps, including binding and three rounds of magnetic nanoparticles washing, were performed using an automated magnetic bead handling device (TANBead Nucleic Acid Extraction System-Maelstrom 8 Autostage). After the final washing step, the liquid phase was transferred from the automation system to a new 1.5 ml tube, followed by the removal of any residual liquid in the tube and air drying on a magnetic rack. Lastly, cfRNA was eluted from the nanoparticles using 10 µl of RNase-free water. Due to the trace amounts of extracted cfRNA, it was directly employed for library preparation.

SLiPiR-seq library preparation
Polyadenylations (poly A) and reverse transcription (RT) reactions were performed simultaneously in a 10 µl reaction system within 0.2 ml PCR tubes. This system consisted of 6.75 µl of cfRNA, 1.25 µl of 4× polyA/RT buffer29, 1 µl of 5 nM RT primer (GenScript Gene Synthesis, 5’-TACGAGATGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTNNNNNNNNTTTTTTTTTTTTTTTTTTTVN-3’), and 1 µl of premixed poly-A polymerase (Enzymatics) and M-MLV reverse transcriptase (Fapon Biotech). The tube was incubated at 37°C for 30 min in a thermocycler. The 8 nt of ‘NNNNNNNN’ in the RT primer sequence represent unique sample barcodes. To deplete any remaining RT primer, 1 µl of Exonuclease I (New England BioLabs) was added after the polyA/RT reaction, followed by incubation at 37°C for 30 min and 80°C for 20 min. Single-stranded cDNA was obtained through denaturation at 95°C for 5 min.

The 11 µl RT reaction mix was then combined with a 9 µl master mix for splint ligation. The master mix comprised 1 µl of T4 DNA ligase (New England BioLabs), 1 µl of 10× T4 DNA Ligase buffer (New England BioLabs), 2 µl of 50% Polyethylene glycol 4000 (Thermo), 2 µl of 10 mM ATP (Enzymatics), 1 µl of 10 µM premixed synthetic adapter (GenScript Gene Synthesis, 5’-CACCTCTCTAUACACUCTTUCCCUACACGACGCTCTUCCGATCUNNNNNN-3’ and 5’-AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTATAGAGAGGTG-3’), and 2 µl of double-distilled water. Splinted ligation occurred by incubating at 20°C for 1 hour and 65°C for 10 min. To digest the blocking strand of the adapter, 1 µl of USER enzyme (New England BioLabs) was added, followed by incubation at 37°C for 15 min.

The libraries were amplified using 2× KAPA HiFi HotStart ReadyMix (Roche Diagnostics) in a 50 µl reaction system with synthetic forward primer (5’-AATGATACGGCGACCACCGAGATCTACACCTCTCTATACACTCTT-3’) and reverse primer (5’-CAAGCAGAAGACGGCATACGAGATGTGACTGGAGTT-3’). Amplification conditions were set as follows: activation at 95°C for 3 min, followed by sixteen cycles of 98°C for 20 s, 60°C for 30 s and 72°C for 30 s, and a final extension of 72°C for 1 min. Amplified libraries were quantified by Qubit assay with the Equalbit 1× dsDNA HS Assay Kit (Vazyme). If the concentration of a sample was less than 10 ng/µl, additional amplification cycles were applied. The amplified libraries were purified using 1.8 times the volume of AMPure XP beads (Beckman Coulter), followed by quantification of the final libraries. To reduce RNA/DNA loss, the entire SLiPiR-seq library preparation workflow was performed within one tube. Libraries with different sample barcodes were pooled together and sequenced on the Illumina NovaSeq 6000 platform (HaploX, Shenzhen) to at least 10 million paired-end reads per sample.

Raw reads processing
The adapter, sample index, and poly-A sequences at the 3’ end of R1 raw sequencing reads were cut away using Cutadapt46. Low-quality reads (q<20) and reads less than 15 nt in length in the R1 file were trimmed away using Trimmomatic47. The first 8 nt of R2 raw sequencing reads (aka sample index), and their identities were extracted using awk. The trimmed R1 fastq files were then split to several subfiles based on the sample indices using a customized python script. The R1 subfiles were used for subsequent reads calling. The R2 files were discarded because their insert reads are of low quality due to the foregoing poly-T sequences. The data processing pipeline was accelerated by executing jobs in parallel in Linux shell using GNU Parallel48.

Establishment of the reference genome and reads calling for tsRNAs, rsRNAs, and ysRNAs
Since there is no existing publicly available reference genome for rsRNAs and ysRNAs, the trimmed clean reads were aligned to the full-length Homosapiens ribosomal RNAs (28S, 18S, 5.8S, and 5S rRNA) and YRNAs (RNY1, RNY3, RNY4, and RNY5) sequences obtained from NCBI with zero mismatch tolerance (--score-min C,0) using Bowtie249. Small RNA fragments mapped to rRNAs and YRNAs were considered as rsRNA and ysRNA, respectively. Ordered ID were named to all unique rsRNAs and ysRNAs, and then their sequences were archived in a FASTA format file. The reference genome of tsRNA was obtained from MINTbase (https://cm.jefferson.edu/MINTbase/)34. Due to the high similarity of the sequences of parent tRNAs, alignment to a database of validated tsRNA sequences allows for more accurate quantification than alignment to a database of full-length tRNAs. The read counts of tsRNAs, rsRNAs, and ysRNAs were called using a custom python script that only counts when a fragment is 100% matched to one of the sequences in the FASTA file.

Reads calling for miRNAs and piRNAs
Reference sequences of miRNAs and piRNAs were acquired from miRbase (https://www.mirbase.org/)50 and piRNABank (http://pirnabank.ibab.ac.in/)51, respectively. The clean reads were first mapped to these references using Bowtie2. Read counts of miRNAs and piRNAs were then called from the output SAM file using awk. Fragments less than 19 nt in length were further trimmed away before miRNA read calling to avoid misclassifications.

Reads calling for lncRNAs, mRNAs, snRNAs, and snoRNAs
The clean reads were mapped to the GRCh38 (hg38) reference genome using Bowtie2. The expression level of mRNAs, lncRNAs, snRNAs, and snoRNAs were called with the GENCODE annotation (release 41, https://www.gencodegenes.org/human/release_41.html) using featureCounts52 from the Subread package. Additional trimming of fragments less than a length of 23 nt was performed before mRNA and lncRNA read calling to reduce the multiple alignment rate.

Quality controls of sequencing results
The called read counts of each sample were merged into an expression matrix using R, followed by reads per million (RPM) normalization and log transformation. The number of total reads, clean reads (≥15 nt), and mappable reads (mapped to GRCh38) were calculated. The sum of reads for each of the nine RNA types, and their percentage ratios were calculated. The number of RNA species detected (RPM > 0) for each RNA type was also calculated. Exclusion criteria of low-quality samples were: (1) Clean reads ratio < 20%; (2) Clean reads < 2 million; (3) rsRNA ratio > 30%; (4) the sum of lncRNA and mRNA ratio > 30%. 9 out of 444 samples were excluded based on these criteria. The reproducibility of SLiPiR-seq was measured by transcriptome-wide Pearson correlation using the normalized read counts of each detected cfRNA.

Technology optimizations
Plasma samples from multiple healthy donors were pooled together for the optimizations of SLiPiR-seq. The SLiPiR-seq library construction was carried out using either the full protocol described above, or a modified protocol with a single variable, including the replacement of ExoI or USER enzyme with distilled water, different final concentrations (1.25 nM, 2.5 nM, 5 nM, 10 nM, 20 nM, 40 nM, and 80 nM) of RT primers, different input plasma volumes (12.5 µl, 25 µl, 50 µl, 100 µl, 200 µl, and 400 µl), different blood standing time before plasma separation (3h, 6h and 9h), and different freeze-thaw cycles of plasma (FT0, FT1, FT2, and FT3). The same number of total reads was randomly subsampled from the sequencing results of different conditions for subsequent reads calling.

Quantitative real-time PCR
Single-stranded cDNA was obtained through the same procedures as described in the SLiPiR-seq library preparation section. qPCR was performed based on the S-Poly(T) Plus method: 1 µl of cDNA template (1:25, 1:10, 1:10, 1:20, 1:2, 1:20, 1:2, 1:20, and 1:10 dilution of original cDNA for lncRNA, miRNA, piRNA, mRNA, snRNA, rsRNA, snoRNA, ysRNA, and tsRNA, respectively), 0.2 µl of AceTaq DNA Polymerase (Vazyme), 2 µl of 10× AceTaq Buffer (Vazyme), 0.4 µl of 10 mM dNTP Mix (Vazyme), 0.2 µl of 100× ROX (Sigma-Aldrich), 0.5 µl of 5× SYBR Green (Roche Diagnostics), 5 µl of 1 µM universal reverse primer (5’-TACGAGATGTGACTGGAGTT-3’), 5 µl of small RNA specific forward primer (Supplementary Data 1), and add double-distilled water up to 20 µl. Amplification conditions were set as follows: activation at 95°C for 5 min, followed by forty cycles of 95°C for 10 s, 60°C for 30 s. Assays were performed on ABI StepOneplus a real-time PCR system (Applied Biosystems).

5’-phosphate modification detection tests
RNA oligos synthesized with or without 5’-phosphate modification (Supplementary Data 2) (GenScript Gene Synthesis) were resuspended at a concentration of 100 µM in TE buffer. Equimolar of RNA oligo was mixed at a final concentration of 10 µM. 1 µl of the mixture was used for library preparation using NEBNext Small RNA Library Prep Set (New England Biolabs) according to the manufacturer’s protocol and SLiPiR-seq as described above.

Calculation of coverage and fragment length distribution
cfRNA was extracted from the pooled plasma for the comparison of NEBNext and SLiPiR-seq. 20 million clean reads were randomly subsampled from the sequencing results of both methods and were used for subsequent reads calling. Fragment length distributions were calculated based on the size and depth of sense-aligned reads. Exon and intron coordinates were extracted from mRNA and lncRNA annotations in GENCODE (release 41). Aligned reads obtained from sam files were used to determine “sense” or “antisense” in exon and intron regions. Gene body coverage of both NEBNext and SLiPiR-seq was evaluated using the ‘geneBody_coverage.py’ function in the RSeQC package with corresponding bam files and GRCh38 reference genome53.

Cross-analysis of tissue-specific genes
Tissue-specific genes were obtained from the Human Protein Atlas database (https://www.proteinatlas.org/humanproteome/tissue) by selecting the ‘The Tissue-specific Proteome’ in the ‘Tissue’ tag. 131 breast, 919 intestine, 441 kidney, 956 liver, 127 lung, 145 ovary, 303 pancreas, 127 prostate, and 307 stomach tissue elevated genes were selected to intersect with mRNAs detected by NEBNext and SLiPiR-seq (Supplementary Data 2).

Sample grouping and partitioning
In the lung cancer study, samples from NOR_SZBA and LC_SZDE were grouped as a discovery cohort (N=245). In the pan-cancer study, samples from NOR_SZBA, LC_SZDE, BRC_NXYK, CRC_NXYK, GC_NXYK and HCC_NXYK were grouped as a discovery cohort (N=382). Samples from NOR_SZDW and LC_SZBU were grouped as a validation cohort (N=53) to validate models established in both the lung cancer detection study and the pan-cancer classification study. The discovery cohort was randomly partitioned into training and test sets in an 80–20% manner using the createDataPartition function in the R package caret54. The training set was used to train models, and the held-out test set was used to assess the models. The independent validation cohort represents unbiased measures of model accuracy, as samples in this cohort were not used for either candidate feature selection or model training and tuning.

Differential expression analyses
The raw count matrices were fitted to a negative binomial distribution using the R package DESeq255. DE cfRNAs were identified at a significance level of 0.1 Benjamini–Hochberg false discovery rate (BH-FDR, Wald test). Differential expression analyses were implemented on (1) SLiPiR-seq cfRNA data of lung cancer patients versus cancer-free controls in the discovery cohort; (2) SliPiR-seq cfRNA data of both the discovery and pan-cancer cohort in a one-class versus-other-classes paradigm. Volcano plots were made by the EnhancedVolcano R package. Pie chart was used to visualize the distribution of different RNA types in the DE cfRNAs. Inclusion criteria for initial candidates are: (1) Upregulated in cancer (fold change > 0); (2) moderate expression level (mean read counts across all samples > 10); (3) large differences between cases and controls (log fold change > 0.8).

Candidate feature selection
Three different strategies were applied to select candidate features for lung cancer detection: (1) Filter method—‘TopN’, based on the rank of significance level; (2) Wrapper method—‘Boruta’, a relevant feature selection wrapper algorithm; (3) Intrinsic method—‘LASSO’, based on least absolute shrinkage and selection operator (LASSO) regularized binomial generalized linear model (GLM). All strategies were applied to each type of cfRNA individually.

In the self-developed filter method, different combinations of the top N cfRNAs ranked by significance level in the differential expression test were examined. Candidate features were selected by first combining the most and the second most significant cfRNAs (N=2), and then sequentially adding the next significant one after another. For example, N=3 means the combinations of top three significant cfRNAs. The best value of N was determined by the mean Area Under the receiver operating characteristic Curve (AUC) in the test set across three machine learning algorithms (see the next section for details). Feature that negatively affected the mean testing set AUC were further filtered out. For example, if AUC of the top three cfRNAs was smaller than AUC of the top two cfRNAs, the third most significant cfRNA was removed from the combinations. The remaining features were verified in the independent validation cohort.

In the wrapper method, Boruta algorithm was applied using the R package Boruta. The variable importance measure was based on RF. The maximal number of important source runs was set to 500. The importance measurement process was repeated 100 times in the randomly partitioned training set. Features evaluated as high importance at least 10 times of the 100 repeats were considered as candidates for verification in the held-out test set and validation cohort.

In the intrinsic method, LASSO regularization of binomial GLM was constructed using the R package glmnet56. ‘Lambda.1se’ determined by 10-fold cross-validation was used as it gives a less complicated model while maintaining a low deviance. LASSO regression coefficient of each cfRNA was calculated for 100 times in the randomly partitioned training set. Features evaluated with a non-zero coefficient at least 10 times of the 100 repeats were considered as candidates. The candidates were further verified in the test set and validation cohort. In the cancer classification analysis, ‘LASSO’ method was used to select candidate features from the six sets of one-class versus-other-classes DE cfRNAs. Over 100 iterations, features evaluated with a non-zero coefficient at least 10 times were considered as candidates.

Machine learning analyses
Three different machine learning algorithms, including logistic regression (LR), random forest (RF), and support vector machine (SVM), were applied in this study using the R packages glmnet, randomForest57, and caret, respectively. For the selection of candidate features, LASSO regularized (alpha=1) LR was used. For the verification of selected candidate features, Ridge regularized (alpha=0) LR was used. 10-fold cross-validation was done to determine the minimum lambda penalty value in all LR models. Random forest models were grown by 500 decision trees with default parameters. Likewise, Support Vector Machines with Linear Kernel (‘svmLinear’) was done using default parameters. All sample partitions and model train-test processes were repeated 100 times to prevent random sampling biases. The performance of models was evaluated by AUC and predicted risk scores of individuals. Receiver operating characteristic (ROC) curve was plotted by the R package pROC. AUC was calculated by the R package ROCR. Risk score was the probability of individual sample predicted by the model. Sensitivity was the ratio of LC patients with a high median risk score to all LC patients. Specificity was the ratio of LC patients who received a high median risk score to all individuals with a high median risk score.

Statistics and reproducibility
All experiments for the technology optimizations of SLiPiR-seq were repeated at least three times to ensure reproducibility. All sample partitions and model training and testing processes were repeated 100 times to avoid random sampling biases. Unless otherwise stated, comparison of means was performed with a two-sided, unpaired Welch’s t-test. Differential expression analyses were performed with Wald test and multiple testing correction (BH-FDR). No statistical method was used to predetermine the sample size. Nine out of 444 samples were excluded due to low sequencing data quality (see the exclusion criteria above). Samples were randomized and blinded during cfRNA extraction, library preparation, and sequencing. Samples were unblinded during data processing and analysis.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "expertise in RNA subtypes\nThe authors tried to develop a sensitive library preparation method that could analyze the profiling\nof cfRNAs independent of their 5’-end modifications in plasma. They have conducted\ntranscriptome-wide characterizations in patients with various cancer types and cancer-free donors,\ncoupled with a bioinformatics analysis pipeline and machine learning analyses, demonstrating its\nability to detect transcriptomic changes occurring in early-stage tumors.\nSome of my concerns are listed as follows:\n1, In the methods of reads calling for tsRNAs, rsRNAs and ysRNAs, the full-length Homo sapiens\nribosomal RNAs and Y RNAs sequences from NCBI were used for alignment, however, reference\ngenome of tsRNA were obtained from the MINTbase. why not use the Homo sapiens transfer RNAs\nsequences, or the reference genome of tRNA that from wellknown GtRNAdb database\n(http://gtrnadb.ucsc.edu/).\n2, Line 380-383, in the characterizations of tsRNAs, rsRNAs and ysRNAs, fragments with identical\nsequences for the first 15 nucleotides from the 5' end were treated as one RNA species. The\nsequences of the subspecies with the longest fragment size were archived into the reference\ngenome under construction.\nIt is interesting, read counts of all subspecies was summed to represent the counts of one RNA\nspecies. But, sometimes, the length of subspecies fragments varies between 16bp and 50bp in one\nRNA species (Supplementary Data 3), this issue is worth thinking about deeply.\n3, In the Fig.3d-f, the visualization of read coverage profiles of diverse cfRNA, the vertical axis\nindicate the RPM values. However, there are large variation in the range of scale values, with\ndifferent index number of x10, it is a little confusing for me. What does it mean differently? Could\nyou use a unified index number?\n4, In the legend, there is no clear definition and description for some abbreviations, such as “R1”\nand “R2” in Fig.S1, “5’-tRHs”, “i-tRFs” and “3’-tRFs” in Fig.3c, “AUC” in Fig.5a. The aim of the\nlegend should be to describe the key messages of the figure, but the figure should also be\ndiscussed in the text. The legend itself should be succinct, while still explaining all symbols and\nabbreviations. May be, you should try to describe the mean of abbreviations in the legend of figure\nor in the supplementary materials.\n5, There are many descriptions for the categories of tsRNA, such as “5’-tRHs”, “i-tRFs”, “3’-tRFs”\nand “3’-tRHs” in the manuscript, and other abbreviations in previous papers. Meanwhile, the\nauthor should discuss the advantage and scientific merit of tsRNAs over other previously explored\nbiomarkers. Some relevant citation could be found, \"DOI:10.1155/2022/8708312 \".\n6, Minor comments: The authors should revise English writing carefully and eliminate some errors\nin the paper to make the paper easier to read.", "reviewer_2": "expertise in cfRNA detection methods (wet lab)\nThe research manuscript titled \"Terminal modifications independent cell-free RNA sequencing\nenables sensitive early cancer detection and classification\" by Jun Wang and colleagues presents a\nnovel and insightful approach for future investigations of small RNAs in cancer. The manuscript\nexplores the potential of small RNAs as diagnostic biomarkers in various cancer. MiRNAs as\npotential diagnostic biomarkers have been intensively explored in various diseases, including\ncancer. Most of the commercially available small RNA library preparation kits are adaptor ligation-\nbased which need the presence of a 3' hydroxyl and 5' phosphate group. As miRNAs are typically\n5’ phosphorylated and 3’ hydroxylated, therefore, these kits were mostly used for miRNA\nsequencing. To address this obstacle, the authors developed the innovative Splint Ligation and\nPhosphate-Independent RNA Sequencing (SLiPiR-seq) method, which allows for the construction of\nlibraries using small RNAs regardless of their 5’ phosphorylation status. The applicability of SLiPiR-\nseq was further validated using synthetic small RNAs, both with and without a 5’-P group, and was\ncompared with a commercially available adapter ligation-based kit (NEBNext Small RNA Library\nPrep Set) in detecting cfRNAs in human plasma samples. The ability of SLiPiR-seq in detecting a\nbroader coverage of the cfRNA repertoire will enable investigation of the functional roles and\ndiagnostic potential of these small RNAs in various disease. Furthermore, the authors should be\ncommented for their excellent work in establishing a robust bioinformatics analysis pipeline and\nimplementing machine learning models to detect and localize different types of cancers using the\nsmall RNA sequencing data. In conclusion, I recommend the publication of this manuscript in the\njournal.\nAs SLiPiR-seq literally could sequence any fragmented RNA, it will bring challenge in following\nanalysis of the sequenced data. However, these challenges could also be opportunity in developing\nnew analysis protocols and identifying new diagnostic markers using these under investigated\ncfRNAs. Additionally, some minor issues in grammar need to be addressed before final acceptance\nas shown below.", "reviewer_3": "expertise in cfRNA bioinformatics\n1. What are the noteworthy results?\nThe authors described a novel method for analyzing cell-free RNAs (cfRNAs) in plasma. The major\nnovelty of the article is:\n• The optimisation of SLiPiR-seq, a library preparation method that is sensitive to cfRNAs with\ndifferent 5’-end modifications and requires only 100 µl of input plasma.\n• The transcriptome-wide characterization of cfRNAs in 165 lung cancer, 30 breast cancer, 37\ncolorectal cancer, 55 gastric cancer, 15 liver cancer, and 133 cancer-free participants, revealing\ntranscriptomic changes in early-stage tumors.\n• The application of machine learning analyses on the differentially expressed cfRNA signatures,\ndemonstrating their robust performance in cancer detection and classification.\n2. Will the work be of significance to the field and related fields? How does it compare to the\nestablished literature? If the work is not original, please provide relevant references.\nTo the best of my knowledge, the author is the first group that demonstrated SLiPiR-seq can\ncapture more RNA species than conventional adapter ligation-based methods. The authors also\nprovided extensive technical characterisation of SLiPiR-seq data through clear visualisations of\nsequencing metrics and comparisons with canonical methods.\nHowever, the SLiPiR-seq methodology described in this study is not significantly different from the\nmethod previously described by Maguire et al. Both methods use splint ligation as a key step to\nligate adapters to RNA fragments without requiring a 5’-phosphate group. Maguire’s team\ncompared SLiPiR-seq against Illumina's TruSeq, New England Biolab's NEBNext, Perkin Elmer's\nNEXTflex methods, while the current study only compared SLiPiR-seq against NEBNext.\nMaguire, S., Lohman, G. J., & Guan, S. (2020). A low-bias and sensitive small RNA library\npreparation method using randomized splint ligation. Nucleic acids research, 48(14), e80-e80.\nThe method developed by Maguire et al was based on the following works.\nNilsen T. W. (2013). Splinted ligation method to detect small RNAs. Cold Spring Harbor protocols,\n2013(1), pdb.prot072611. https://doi.org/10.1101/pdb.prot072611\nJin, J., Vaud, S., Zhelkovsky, A. M., Posfai, J., & McReynolds, L. A. (2016). Sensitive and specific\nmiRNA detection method using SplintR Ligase. Nucleic acids research, 44(13), e116.\nhttps://doi.org/10.1093/nar/gkw399\nProper discussion and citation of the above works would allow readers differentiate this work from\nestablished literature. Besides, the protocol optimisations made by the authors, presumably\nallowing the method to be more sensitive to low quantity of cfRNA, should be clearly highlighted in\nthe main text.\n3. Does the work support the conclusions and claims, or is additional evidence needed?\nIn line 115, the authors noted that “the fragment size of msRNAs and lsRNAs exhibited smoother\ndistribution and broader coverage (23-60 nt)”. Could the authors apply appropriate statistical tests\nto substantiate the claim?\nIn line 170, the authors noted “a significant increase in the cumulative expression level in early\nstage LC patients compared to controls”. Could this observation due to confounding factors, bias in\nnormalisation or multi-mapped reads? A recent study found that common normalization methods\nvary across different datasets.\nDüren, Y., Lederer, J., & Qin, L. X. (2022). Depth normalization of small RNA sequencing: using\ndata and biology to select a suitable method. Nucleic acids research, 50(10), e56.\nhttps://doi.org/10.1093/nar/gkac064\nIn line 185-261, the authors attempted to use SLiPiR-seq data to develop a cfRNA signature for\nearly detection of cancer. However, given by the relatively small sample size and diverse cancer\ntypes, it is unclear if the study has enough statistical power to substantiate the claim in line 260.\n“Collectively, these results demonstrate that cfRNA signatures identified by SLiPiR-seq can\nprecisely discriminate different cancer types”.\n4. Are there any flaws in the data analysis, interpretation and conclusions? Do these prohibit\npublication or require revision?\nIn Fig. 2a, it looks like the dots corresponding to miRNAs are all below the trend line. Could it\nmean that miRNA expression is most likely underestimated by SLiPiR-seq?\nIn line 140, the authors classified tsRNAs, rsRNAs and ysRNAs species based on the first 15\nnucleotides. I would expect more comprehensive comparison of RNA species using established\ndatabases such as Rfam and RNAcentral.\nIn Fig. 4e, the legends are unclear. The colour blocks below the dendrogram, on the left of the\nheatmap, at the bottom of the heatmap should be clearly labelled. Description of blocks in grey is\nnot defined.\nIn line 168, should absolute log2 fold change < 1 be used instead?\nThe surprisingly high sensitivity of 96.15% and specificity of 100% may indicate overfitting under\nhigh dimensionality of the features. While it may not be feasible to have independent training,\ntesting and validation samples, the use of regularisation techniques and cross-validation may\nalleviate overfitting of features.\n5. Is the methodology sound? Does the work meet the expected standards in your field?\nYes in general, but please note the flaws above.\n6. Is there enough detail provided in the methods for the work to be reproduced?\nYes"}, "author_response": "RESPONSE TO REVIEWERS' COMMENTS\nReviewer #1\nQuestion1: In the methods of reads calling for tsRNAs, rsRNAs and ysRNAs, the full-\nlength Homo sapiens ribosomal RNAs and Y RNAs sequences from NCBI were used for\nalignment, however, reference genome of tsRNA were obtained from the MINTbase. why\nnot use the Homo sapiens transfer RNAs sequences, or the reference genome of tRNA that\nfrom wellknown GtRNAdb database (http://gtrnadb.ucsc.edu/).\nAnswer: We appreciate the reviewer raising the question about our rationale for using\nMINTbase rather than GtRNAdb or Homo sapiens transfer RNA sequences for characterization\nof tsRNA in plasma. To provide robust quantification and classification of tsRNA profiles, we\nopted to use MINTbase over GtRNAdb for the following reasons:\n1) When aligning to rRNA and Y RNA reference genomes, approximately 95% of mapped reads\nwere unique alignments. However, due to the highly similar sequences for the ~600 tRNAs,\nmapping to GtRNAdb resulted in only ~1% uniquely aligned reads, with ~14% multiple\nmapped reads (Table 1). The high rate of multi-mapping with GtRNAdb precludes reliable\nquantification and classification of individual tRNA species.\nTable 1 Reads mapped to the GtRNAdb tRNA reference genome (hg38-tRNAs.fa).\nTotal clean reads Unique mapped reads Multiple mapped reads\nSLiPiR-seq 20M ~ 0.23M (1.13%) ~ 2.84M (14.19%)\nNEBNext 20M ~ 0.002M (0.01%) ~ 0.08M (0.40%)\n2) MINTbase contains over 28,824 non-redundant, validated tsRNA sequences with\ncomprehensive annotations (length, category, amino acid, anticodon, source, etc). This allowed\nrobust classification of plasma tsRNA profiles not possible by alignment to full-length tRNAs\nin GtRNAdb.\nPliatsika, V. et al. MINTbase v2.0: a comprehensive database for tRNA-derived fragments that\nincludes nuclear and mitochondrial fragments from all The Cancer Genome Atlas projects.\nNucleic Acids Res 46, D152-D159, doi:10.1093/nar/gkx1075 (2018).\n3) Critically, when using MINTbase for tsRNA quantification, we were able to classify nearly\nall (14% vs 15%) of the reads aligning to the tRNA genome (hg38-tRNAs.fa). This indicates\nthat MINTbase captures the majority of authentic tsRNA sequences in plasma.\nIn summary, due to the sequence similarity among tRNAs, alignment to a database of validated\ntsRNA sequences in MINTbase enables more accurate quantification and classification\ncompared to alignment to full-length tRNAs. The exhaustively curated MINTbase tsRNA\nreference provides a more suitable foundation for characterization of the plasma tsRNA\nrepertoire. We agree with the reviewer that justification of our database choice will add rigor to\nour tsRNA quantification bioinformatic analysis, and we added “Due to the high similarity of\nthe sequences of parent tRNAs, alignment to a database of validated tsRNA sequences allows\nfor more accurate quantification than alignment to a database of full-length tRNAs.” to the\nmethod section to validate our choice.\nQuestion2: Line 380-383, in the characterizations of tsRNAs, rsRNAs and ysRNAs,\nfragments with identical sequences for the first 15 nucleotides from the 5' end were treated\nas one RNA species. The sequences of the subspecies with the longest fragment size were\narchived into the reference genome under construction.\nIt is interesting, read counts of all subspecies was summed to represent the counts of one\nRNA species. But, sometimes, the length of subspecies fragments varies between 16bp and\n50bp in one RNA species (Supplementary Data 3), this issue is worth thinking about\ndeeply.\nAnswer: We agree with the reviewer’s concern about the classification of RNA species. We\nhave removed the analytical step that classifies fragments with identical sequences for the first\n15 nucleotides from the 5' end as one RNA species. To avoid too many RNA species, we also\ntweaked the codes for bowite2 alignment to tolerate zero mismatch. As a result, we detected a\ntotal of 45397 rsRNA and 2664 ysRNA unique sequences.\nHowever, the revision on the classification of RNA species inevitably altered other results\ndemonstrated in the paper, resulting in extensive modifications. For example, the number of\ndetected RNA species for tyr-sRNA in fig. 1c and fig. 2d were updated. Fig. 3a was also updated.\nAfter revision, a total of 152858 cfRNA species were used for the differential expression\nanalysis, which is 1.5 times more than before. This will affect the normalization of the cfRNA\nread count matrix and the calculation of adjusted p values during the DESeq2 analysis.\nTherefore, all results from fig. 4 to fig. 6 and their supplementary figures were inevitably\nchanged. After revision, rsRNA accounts for 73.4% of all differentially expressed cfRNA and\novershadows the performance of other RNA types. To compensate for this issue, we adjusted\nthe inclusion criteria of candidate cfRNA selection from log2 fold change>1 to >0.8, which\nallows more cfRNAs of other types to be applied for the subsequent feature selection test.\nOverall, although most results were updated, the inclusion of more tyr-sRNA species does not\nalter the findings and conclusions of our study. Interestingly, LR and SVM models training with\n50 rsRNA (SVM AUC= 0.819 [0.786-0.843]) or 25 ysRNA (SVM AUC= 0.793 [0.747-0.829])\nfeatures showed high AUC in LC detection. However, models trained by the RF algorithm\nshared no agreement with these results (Supplementary Fig. 11c). Therefore, rsRNA and\nysRNA were not used in the combination test. Again, the five RNA types of interest (mRNAs,\nmiRNAs, snRNAs, snoRNAs, and tsRNAs) were used in the combination test and the cancer\nclassification test.\nQuestion3: In the Fig.3d-f, the visualization of read coverage profiles of diverse cfRNA,\nthe vertical axis indicate the RPM values. However, there are large variation in the range\nof scale values, with different index number of x10, it is a little confusing for me. What\ndoes it mean differently? Could you use a unified index number?\nAnswer: We appreciate the reviewer raising this point about the inconsistent scaling between\npanels in Figures 3d-f, which was unintentionally confusing. To provide a clear comparison of\nthe read coverage profiles between SLiPiR-seq and NEBNext for the different RNA types, we\nhave revised the figures (shown below) to display the data for both methods in each panel with\na unified y-axis scale. The updated figures now clearly show the difference in read coverage\ndepth and breadth between the two methods on the same scale. For example, the much greater\nsequencing depth achieved by SLiPiR-seq compared to NEBNext for 5S rRNA, 5.8S rRNA,\n18S rRNA, 28S rRNA, genomic tRNA, mitochondrial tRNA and RNY5 is now evident in the\nrevised figure with the unified scale.\nQuestion4: In the legend, there is no clear definition and description for some\nabbreviations, such as “R1” and “R2” in Fig.S1, “5’-tRHs”, “i-tRFs” and “3’-tRFs” in\nFig.3c, “AUC” in Fig.5a. The aim of the legend should be to describe the key messages of\nthe figure, but the figure should also be discussed in the text. The legend itself should be\nsuccinct, while still explaining all symbols and abbreviations. May be, you should try to\ndescribe the mean of abbreviations in the legend of figure or in the supplementary\nmaterials.\nAnswer: We appreciate the reviewer raising this important point about clearly defining\nabbreviations and symbols used in our figures and legends. As suggested, we have thoroughly\nreviewed all figure legends to ensure abbreviations and acronyms are defined at first use.\nFor example, in the caption of Figure S1, we added “R1 and R2 indicate the forward reads and\nreverse reads from paired-end sequencing, respectively.”. In the caption of Figure 3c, we added\n“tRHs refers to tRNA-halves and tRFs refers to tRNA-derived fragments. i-tRFs refers to\ninternal tRNA-derived fragments.”. In the caption of Figure 5a, we defined AUC as Area Under\nthe receiver operating characteristic Curve.\nQuestion5: There are many descriptions for the categories of tsRNA, such as “5’-tRHs”,\n“i-tRFs”, “3’-tRFs” and “3’-tRHs” in the manuscript, and other abbreviations in previous\npapers. Meanwhile, the author should discuss the advantage and scientific merit of\ntsRNAs over other previously explored biomarkers. Some relevant citation could be found,\n\"DOI:10.1155/2022/8708312 \".\nAnswer: We added a paragraph discussing the advantage and scientific merit of tsRNAs in the\ndiscussion section. “Recent studies have reported that tsRNAs can regulate cancer progression\nat the post-transcriptional level through multiple mechanisms and are thus considered as critical\nregulators and biomarkers of cancer. When we sorted the average specificities of different\ncfRNA combinations in the discovery and validation cohort, we found that tsRNAs were\npresent in all top five combinations (\"m+sn+sno+ts”, “m+mi+ts”, “m+sn+ts”, “sn+ts”, “m+ts”)\n(Supplementary Data 6). This finding suggested that the inclusion of tsRNAs in the\ncombination panel significantly enhances the specificity of the machine learning models in lung\ncancer detection. This heightened specificity is of paramount importance in early cancer\nscreening, where the cost and psychological implications of false positives can be substantial.”\nFu, M. et al. Emerging roles of tRNA-derived fragments in cancer. Mol Cancer 22, 30,\ndoi:10.1186/s12943-023-01739-5 (2023).\nLiu, X. et al. A pro-metastatic tRNA fragment drives Nucleolin oligomerization and\nstabilization of its bound metabolic mRNAs. Mol Cell 82, 2604-2617 e2608,\ndoi:10.1016/j.molcel.2022.05.008 (2022).\nRen, J. et al. The tRNA-Cys-GCA Derived tsRNAs Suppress Tumor Progression of Gliomas\nvia Regulating VAV2. Dis Markers 2022, 8708312, doi:10.1155/2022/8708312 (2022).\nQuestion6: Minor comments: The authors should revise English writing carefully and\neliminate some errors in the paper to make the paper easier to read.\nAnswer: Thank you for your careful review and for pointing out the grammatical errors in our\nmanuscript. To address these issues, we have utilized the DeepL tool\n(https://www.deepl.com/write) to assist in refining the English writing of our manuscript. We\nbelieve that this tool has significantly enhanced the clarity and coherence of our text. We kindly\nrequest that you refer to the tracked changes version of the manuscript, where all the\nmodifications have been highlighted.\nOther revisions we made:\n1. Renamed all ‘msRNA’ to ‘mRNA’ and all ‘lsRNA’ to ‘lncRNA’ to avoid confusion.\n2. Added ‘Supplementary Data 6- Summary of different cfRNA combinations’ to show\nAUC, sensitivities, specificities of different combinations.\n3. Replaced fig. 5c. To avoid confusion between sensitivity at 100% specificity and the\nsensitivity/ sensitivity calculated from risk scores, we determined to use AUC as\nindicator of model accuracy throughout the paper instead of using sensitivity at 100%\nspecificity.\n4. Reordered supplementary fig. 12-14. Replaced the cumulative read count boxplot in\nsupplementary fig. 14 to risk score boxplot of individuals. The presentation of the\ncumulative read count boxplot is less meaningful than the read count boxplot of\nindividual cfRNA in supplementary fig. 12, and their interpretations are duplicative.\nReviewer #2\nAs SLiPiR-seq literally could sequence any fragmented RNA, it will bring challenge in\nfollowing analysis of the sequenced data. However, these challenges could also be\nopportunity in developing new analysis protocols and identifying new diagnostic markers\nusing these under investigated cfRNAs. Additionally, some minor issues in grammar need\nto be addressed before final acceptance as shown below.\nWe greatly appreciate the reviewer recognizing the novelty and potential impact of the SLiPiR-\nseq method. We agree that developing new analysis approaches for this novel cfRNA\nsequencing data represents an exciting area for future work, with potential to uncover\npreviously underexplored RNA biomarker candidates.\nWe sincerely appreciate the reviewer highlighting the minor grammar issues present in the\ninitial draft. We have thoroughly proofread the manuscript and addressed these concerns as\nsuggested.\nAll changes were incorporated in the revised manuscript with track changes.\nOther revisions we made:\n1. Renamed all ‘msRNA’ to ‘mRNA’ and all ‘lsRNA’ to ‘lncRNA’ to avoid confusion.\n2. Added ‘Supplementary Data 6- Summary of different cfRNA combinations’ to show\nAUC, sensitivities, specificities of different combinations.\n3. Replaced fig. 5c. To avoid confusion between sensitivity at 100% specificity and the\nsensitivity/ sensitivity calculated from risk scores, we determined to use AUC as\nindicator of model accuracy throughout the paper instead of using sensitivity at 100%\nspecificity.\n4. Reordered supplementary fig. 12-14. Replaced the cumulative read count boxplot in\nsupplementary fig. 14 to risk score boxplot of individuals. The presentation of the\ncumulative read count boxplot is less meaningful than the read count boxplot of\nindividual cfRNA in supplementary fig. 12, and their interpretations are duplicative.\nReviewer #3\nQuestion 1: To the best of my knowledge, the author is the first group that demonstrated\nSLiPiR-seq can capture more RNA species than conventional adapter ligation-based methods.\nThe authors also provided extensive technical characterisation of SLiPiR-seq data through\nclear visualisations of sequencing metrics and comparisons with canonical methods.\nHowever, the SLiPiR-seq methodology described in this study is not significantly different\nfrom the method previously described by Maguire et al. Both methods use splint ligation as a\nkey step to ligate adapters to RNA fragments without requiring a 5’-phosphate group.\nMaguire’s team compared SLiPiR-seq against Illumina's TruSeq, New England Biolab's\nNEBNext, Perkin Elmer's NEXTflex methods, while the current study only compared\nSLiPiR-seq against NEBNext.\nMaguire, S., Lohman, G. J., & Guan, S. (2020). A low-bias and sensitive small RNA library\npreparation method using randomized splint ligation. Nucleic acids research, 48(14), e80-e80.\nThe method developed by Maguire et al was based on the following works.\nNilsen T. W. (2013). Splinted ligation method to detect small RNAs. Cold Spring Harbor\nprotocols, 2013(1), pdb.prot072611. https://doi.org/10.1101/pdb.prot072611\nJin, J., Vaud, S., Zhelkovsky, A. M., Posfai, J., & McReynolds, L. A. (2016). Sensitive and\nspecific miRNA detection method using SplintR Ligase. Nucleic acids research, 44(13), e116.\nhttps://doi.org/10.1093/nar/gkw399\nProper discussion and citation of the above works would allow readers differentiate this work\nfrom established literature. Besides, the protocol optimisations made by the authors,\npresumably allowing the method to be more sensitive to low quantity of cfRNA, should be\nclearly highlighted in the main text.\nAnswer: We appreciate the reviewer for raising the pertinent point regarding the use of splint\nligation for library preparation in previous studies. While both SLiPiR-seq and Maguire et al.’s\nmethod utilize splint ligation as a central step for adapter ligation, it is crucial to delineate the\nspecific distinguishing features between these two approaches.\nFirstly, in the method proposed by Maguire et al., they optimized their workflow to ligate a 5'-splint\nadapter containing six degenerate nucleotides to the 5' end of RNA prior to reverse transcription.\nThis approach effectively enhances ligation efficiency compared to traditional single-strand RNA-\nRNA base adapter ligation, thereby increasing sensitivity, especially for low-input RNA libraries\n(shown below). However, it's noteworthy that even in this method, 5'-phosphate modification of\nRNA remains a prerequisite for successful adapter ligation.\nIn contrast, SLiPiR-seq introduces several distinctive attributes (shown below):\n1) Reverse transcription prior to ligation: In SLiPiR-seq, reverse transcription is performed as the\ninitial step. Important sequencing elements such as sample barcode and the P7 primer sequences\nare added to the first-strand cDNA after reverse transcription.\n2) DNA-based splint adapter: SLiPiR-seq employs a DNA-based splint adapter that is capable of\nligating to the 3' end of the first-strand cDNA, corresponding to the 5' end of the RNA before\nreverse transcription. Consequently, it circumvents the necessity for a 5'-phosphate modification\non the RNA, ensuring compatibility with RNA molecules lacking this modification.\nTo address the reviewer's valid concern and highlight these distinctions, we have incorporated a\ndiscussion within the manuscript, supplemented by appropriate citations of the relevant works by\nMaguire et al. (2020) and Nilsen (2013).\n“The use of splint ligation for small RNA library preparation has been reported in previous studies.\nIt is crucial to delineate the specific distinguishing features between SLiPiR-seq and the previous\napproach. Maguire et al. employed a novel RNA splint ligation adapter containing six degenerate\nnucleotides at the 5' end of the RNA. This approach effectively increases ligation efficiency\ncompared to traditional single-stranded RNA-RNA base adapter ligation, thereby increasing\nsensitivity, especially for low-input RNA libraries. However, it's noteworthy that even with this\nmethod, 5'-phosphate modification of the RNA remains a prerequisite for successful adapter\nligation. In SLiPiR-seq, reverse transcription is performed prior to adapter ligation. Important\nsequencing elements such as sample barcode and the P7 primer sequences are added to the first-\nstrand cDNA after reverse transcription. SLiPiR-seq employs a DNA-based splint adapter that is\ncapable of ligating to the 3' end of the cDNA, corresponding to the 5' end of the RNA. Consequently,\nthis bypasses the necessity for a 5'-phosphate modification on the RNA.”\nComparison of the splint ligation strategies utilized in the method developed by Maguire et al.\nand SLiPiR-seq.\nQuestion 2: In line 115, the authors noted that “the fragment size of msRNAs and lsRNAs\nexhibited smoother distribution and broader coverage (23-60 nt)”. Could the authors apply\nappropriate statistical tests to substantiate the claim?\nAnswer: We appreciate the reviewer raising this important point about substantiating our claim\nregarding the fragment size distribution and coverage of mRNAs and lncRNAs. To statistically\nsupport this observation, we conducted further analysis using the median and median absolute\ndeviation (MAD) to characterize the distribution of fragment sizes between the NEBNext and\nSLiPiR-seq datasets as shown below (Table 1).\nTable1 Statistical comparison of mRNA and lncRNA fragment size profiles detected by NEBNext\nand SLiPiR-seq.\nNEBNext SLiPiR-seq\nmRNA 31bp ± 0bp 36bp ± 14.8bp\nlncRNA 31bp ± 8.9bp 33bp ± 11.9bp\nIn summary, the larger median fragment size values combined with the greater MAD provide\nstatistical evidence that SLiPiR-seq detected mRNAs and lncRNAs with a smoother distribution\nand broader size range compared to NEBNext. We have incorporated these quantitative results in\nthe revised manuscript to support our initial qualitative observation.\nQuestion 3: In line 170, the authors noted “a significant increase in the cumulative expression\nlevel in early stage LC patients compared to controls”. Could this observation due to\nconfounding factors, bias in normalisation or multi-mapped reads? A recent study found that\ncommon normalization methods vary across different datasets.\nDüren, Y., Lederer, J., & Qin, L. X. (2022). Depth normalization of small RNA sequencing:\nusing data and biology to select a suitable method. Nucleic acids research, 50(10), e56.\nhttps://doi.org/10.1093/nar/gkac064\nAnswer: We greatly appreciate the reviewer's insightful inquiry into the critical issue of data\nnormalization in small RNA sequencing studies and their reference to the work of Düren et al.\n(2022) on depth normalization in small RNA sequencing. Normalization methods play a crucial\nrole in accurately characterizing gene expression, and we acknowledge the potential impact of\nconfounding factors, normalization biases, and multi-mapped reads in our analysis.\nAt the outset of our study, we recognized the significance of choosing an appropriate normalization\nmethod for our cfRNA data. Given that plasma cfRNA are small RNA fragments, and not subjected\nto manual fragmentation during library preparation, we hypothesized that Reads Per Million (RPM)\nnormalization would be a suitable choice. To validate this hypothesis, we conducted rigorous\nvalidation experiments using the gold-standard technique for RNA quantification. We profiled 181\ndistinct plasma cfRNAs and compared their expression levels normalized by RPM (mapped reads\nas indicator of RPM) to quantitative real-time PCR (qPCR) results. We observed a strong overall\nconcordance (R=0.86) between RPM-normalized expression levels and Ct values measured by\nqPCR (Fig. 2a). This high concordance was consistently observed for various RNA species,\nincluding lncRNA (R=0.83), miRNA (R=0.81), piRNA (R=0.71), mRNA (R=0.87), snRNA\n(R=0.89), rsRNA (R=0.95), snoRNA (R=0.82), ysRNA (R=0.97), and tsRNA (R=0.85)\n(Supplementary Fig5).These findings provide strong evidence for the reliability of RPM\nnormalization in our study.\nTo further address concerns related to potential normalization bias within cohort-based study\nutilizing SLiPiR-seq, we conducted correlation analysis on cfRNA profiles from cancer patients\nand cancer-free individuals within our dataset using RPM normalization method. This analysis\nrevealed a high overall correlation (R=0.987, p<2.2e-16) between the two groups (Fig. 4a), further\nsupporting the robustness and consistency of the normalization method across different sample sets.\nConcerning the issue of multi-mapped reads, we implemented specific measures to mitigate\npotential biases. For RNA species such as rsRNA, ysRNA, and tsRNA, we utilized a custom Python\nscript designed for sequence-based exact matching to minimize any potential bias introduced by\nmultiple-mapping issues (Fig. 3a and Supplementary Fig. 1). For miRNA and piRNA, we employed\nthe bowtie2 mapping algorithm with stringent parameters to reduce the presence of multi-mapped\nreads (Supplementary Fig. 1). Furthermore, for mRNA, lncRNA, snoRNA, and snRNA, we\nimplemented a pre-processing step to generate a distinct Gene Transfer Format (GTF) annotation\nfile for each RNA species, facilitating precise annotation of genomic locations. We then adopted\nstrand-specific (-s, 1) read calling strategies in FeatureCounts, ensuring the accurate assignment of\nreads to their respective RNA species and strands (Supplementary Fig. 1). These efforts have\nsignificantly reduced the incidence of multi-mapped rates.\nIn conclusion, we consider that the observed increase in cumulative expression levels in early stage\nLC patients is primarily attributed to the upregulated expression of specific cfRNAs in cancer\npatients. Despite the challenges associated with data normalization, our comprehensive validation\nexperiments and correlation analyses support the reliability and accuracy of our normalization\nprocess. This allows us to interpret the observed differences in cfRNA expression between groups\nas indicative of biological variations associated with lung cancer.\nQuestion 4: In line 185-261, the authors attempted to use SLiPiR-seq data to develop a cfRNA\nsignature for early detection of cancer. However, given by the relatively small sample size and\ndiverse cancer types, it is unclear if the study has enough statistical power to substantiate the\nclaim in line 260. “Collectively, these results demonstrate that cfRNA signatures identified by\nSLiPiR-seq can precisely discriminate different cancer types”.\nAnswer: We updated the conclusion of this section into “Collectively, these results demonstrate that\ncfRNA signatures identified by SLiPiR-seq are promising in the classification of different cancer\ntypes, motivating us to validate the clinical utility of SLiPiR-seq in larger sample cohorts in future\ninvestigations.”\nQuestion 5: In Fig. 2a, it looks like the dots corresponding to miRNAs are all below the trend\nline. Could it mean that miRNA expression is most likely underestimated by SLiPiR-seq?\nAnswer: We appreciate the reviewer for bringing up this important observation. Upon careful\nreevaluation of our data, we have identified an error in placing the trendline in the original\npresentation of Fig. 2a. We have since corrected the trendline to accurately represent the data.\nFurthermore, to address the concern regarding the potential underestimation of miRNA expression\nby SLiPiR-seq, we conducted a thorough analysis using the data provided in Supplementary Data\n1. We specifically focused on the relationship between Ct values and log2RPM for miRNA\ndetection, as well as the comparison of miRNA expression with that of other RNA species. Our\nanalysis revealed that when considering only miRNAs, there are 31 data points located above the\ntrendline and 41 below it. In contrast, when analyzing the combined expression of the other 8 RNA\nspecies alongside miRNAs, there are 22 data points above and 50 below the trendline. These\nfindings indicate that relative to the other RNA species, miRNA expression might be slightly\nunderestimated by SLiPiR-seq.\nIt is important to emphasize that while there might be a slight underestimation of miRNA\nexpression relative to other RNA species, this observation does not compromise the integrity of our\ndownstream differential analysis. Specifically, our analysis methodology does not involve the direct\ncomparison of the abundance of each RNA species within the cfRNA repertoire detected by\nSLiPiR-seq. Instead, we focus on the relative expression of individual miRNAs, which we have\nfound to be highly correlated with qPCR results.\nQuestion 6: In line 140, the authors classified tsRNAs, rsRNAs and ysRNAs species based on\nthe first 15 nucleotides. I would expect more comprehensive comparison of RNA species using\nestablished databases such as Rfam and RNAcentral.\nAnswer: The other reviewer also pointed out a question about the classification of RNA species on\nthe feedback: “It is interesting, read counts of all subspecies was summed to represent the counts of\none RNA species. But, sometimes, the length of subspecies fragments varies between 16bp and\n50bp in one RNA species (Supplementary Data 3), this issue is worth thinking about deeply.”\nTo address this issue, we removed the analytical step that classifies fragments with identical\nsequences for the first 15 nucleotides from the 5' end as one RNA species. To avoid too many RNA\nspecies, we also tweaked the codes for bowite2 alignment to tolerate zero mismatch. As a result, we\ndetected a total of 45397 rsRNA and 2664 ysRNA unique sequences.\nHowever, the revision on the classification of RNA species inevitably altered other results\ndemonstrated in the paper, resulting in extensive modifications. For example, the number of\ndetected RNA species for tyr-sRNA in fig. 1c and fig. 2d were updated. Fig. 3a was also updated.\nAfter revision, a total of 152858 cfRNA species were used for the differential expression analysis,\nwhich is 1.5 times more than before. This will affect the normalization of the cfRNA read count\nmatrix and the calculation of adjusted p values during the DESeq2 analysis. Therefore, all results\nfrom fig. 4 to fig. 6 and their supplementary figures were inevitably changed. After revision, rsRNA\naccounts for 73.4% of all differentially expressed cfRNA and overshadows the performance of\nother RNA types. To compensate for this issue, we adjusted the inclusion criteria of candidate\ncfRNA selection from log2 fold change>1 to >0.8, which allows more cfRNAs of other types to be\napplied for the subsequent feature selection test.\nOverall, although most results were updated, the inclusion of more tyr-sRNA species does not alter\nthe findings and conclusions of our study. Interestingly, LR and SVM models training with 50\nrsRNA (SVM AUC= 0.819 [0.786-0.843]) or 25 ysRNA (SVM AUC= 0.793 [0.747-0.829])\nfeatures showed high AUC in LC detection. However, models trained by the RF algorithm shared\nno agreement with these results (Supplementary Fig. 11c). Therefore, rsRNA and ysRNA were not\nused in the combination test. Again, the five RNA types of interest (mRNAs, miRNAs, snRNAs,\nsnoRNAs, and tsRNAs) were used in the combination test and the cancer classification test.\nBack to this question, Rfam and RNAcentral are well-known and established databases, and they\nhave comprehensive information about tRNAs, rRNAs and yRNAs. However, they do not contain\ninformation about tsRNAs, rsRNAs and ysRNAs. Therefore, we are not able to conduct comparison\nof these small RNAs using Rfam and RNAcentral.\nQuestion 7: In Fig. 4e, the legends are unclear. The colour blocks below the dendrogram, on\nthe left of the heatmap, at the bottom of the heatmap should be clearly labelled. Description\nof blocks in grey is not defined.\nAnswer: We added ‘left’, ‘top’ and ‘bottom’ on the figure legend of fig. 4e to make the color blocks\nclear. We also removed the gray color from cancer-free controls for the bottom color block.\nQuestion 8: In line 168, should absolute log2 fold change < 1 be used instead?\nAnswer: It depend very much on what we are looking for in the study. If we are screening for\ninternal controls for qPCR quantification, absolute log2 fold change < 1 should be used. On the\nother hand, if we aim to identify cancer biomarkers, we believe that log2 fold change > 1 should\nbe considered. This is because the detection of cfRNA signatures with a larger fold change between\ncases and controls will be easier in a cancer detection test.\nQuestion 9: The surprisingly high sensitivity of 96.15% and specificity of 100% may indicate\noverfitting under high dimensionality of the features. While it may not be feasible to have\nindependent training, testing and validation samples, the use of regularisation techniques and\ncross-validation may alleviate overfitting of features.\nAnswer: We agree with the reviewer’s concern about the overfitting of machine learning models in\nthe validation cohort. In the original submission, while we applied regularization and cross-\nvalidation in our models, we selected the cfRNAs with the highest AUC in the validation cohort\namong three different feature selection methods for the combination of multiple types of cfRNAs,\nwhich may result in overfitting of the model in the validation cohort and compromising the purpose\nof validation. In our revision, we used cfRNAs with the highest AUC in the test set of the discovery\ncohort among three different feature selection methods for the combination test. The validation\nprocess of the combination test is now unbiased. After revision, the best combination\n(“m+sn+sno+ts” panel) shows 100% specificity and 99.28% sensitivity in the discovery cohort and\n95.24% specificity and 76.92% sensitivity in the validation cohort.\nOther revisions we made:\n1. Renamed all ‘mRNA’ to ‘mRNA’ and all ‘lncRNA ’ to ‘lncRNA’ to avoid confusion.\n2. Added ‘Supplementary Data 6- Summary of different cfRNA combinations’ to show AUC,\nsensitivities, specificities of different combinations.\n3. Replaced fig. 5c. To avoid confusion between sensitivity at 100% specificity and the\nsensitivity/ sensitivity calculated from risk scores, we determined to use AUC as indicator\nof model accuracy throughout the paper instead of using sensitivity at 100% specificity.\n4. Reordered supplementary fig. 12-14. Replaced the cumulative read count boxplot in\nsupplementary fig. 14 to risk score boxplot of individuals. The presentation of the\ncumulative read count boxplot is less meaningful than the read count boxplot of individual\ncfRNA in supplementary fig. 12, and their interpretations are duplicative."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I am satisfied with the Author’s response as well as the revision of this manuscript, and have no\nmore comments.", "reviewer_2": "The authors have addressed all my concerns. This work disclosed the necessity of exploring\nmultiple species of RNAs in plasma as biomarkers for early-detection of cancer. I recommend\npublication of the corrected manuscript in Nature Communications.", "reviewer_3": "Thanks for addressing the comments in my first review. Please refer to my further comments\nbelow:\nQs1: No further comments\nQs2: Median and median absolute deviation are not good parameters to substantiate if a\ndistribution is smoother and broader. Please consider using a Goodness-of-Fit test to check the\nsmoothness of a distribution. Meanwhile, kurtosis is a common measure of the breadth\n(tailedness) of a distribution.\nQs3: The authors showed that the RPM values correlates with qPCR Ct values, which supports the\nreliability of RPM normalisation. Ideally, artificial spike-ins that corresponds to the cell number or\nDNA quantity can be used to measure if a global difference in cumulative expression level exists in\nfuture studies.\nThe authors further mentioned the use of sequence-based exact matching and stringent mapping\nparameters to reduce the chance of multi-mapping. Since some RNA species have very similar\nsequences (e.g. miRNA), these filters may affect the sensitivity of detection. Is there any figure or\ntable that summarises the impact of these filtering strategies on the number of detected RNA\nspecies?\nQs4: No further comments\nQs5: What are the potential reasons that might lead to the underestimation of miRNA expression\nby SLiPiR-seq? In relation to my comments in Qs3, could the stringent filtering strategies affect\nthe sensitivity of counting miRNAs?\nQs6: No further comments\nQs7: No further comments\nQs8: No further comments\nQs9: Please show the number of samples in the discovery, testing, and validation cohort. Besides,\nplease indicate the cohort splitting strategy in detail, such as whether random sampling or\nstratified sampling were used."}, "author_response": "RESPONSE TO REVIEWERS’ COMMENTS\nReviewer #3\nThanks for addressing the comments in my first review. Please refer to my further comments below:\nQs1: No further comments\nQs2: Median and median absolute deviation are not good parameters to substantiate if a distribution\nis smoother and broader. Please consider using a Goodness-of-Fit test to check the smoothness of\na distribution. Meanwhile, kurtosis is a common measure of the breadth (tailedness) of a distribution.\nAnswer: We appreciate the insightful suggestion from the reviewer regarding the statistical analysis\nof fragment size distribution. In response to this, we conducted a comprehensive analysis of the\nmRNA and lncRNA fragment size distributions, fitting them to the Gamma distribution and\ncalculating the associated α and β parameters. The obtained results contribute valuable insights into\nthe smoothness and breadth of these distributions.\nThe fragment size distributions for mRNAs in both NEBNext (α=3.32, β=0.10) and SLiPiR-seq\n(α=10.22, β=0.24) were scrutinized, revealing a difference in smoothness, with SLiPiR-seq\nexhibiting a smoother distribution. A similar trend was observed in lncRNA fragment size\ndistributions, where SLiPiR-seq (α=12.07, β=0.35) demonstrated a smoother distribution compared\nto NEBNext (α=6.29, β=0.21).\nAdditionally, we utilized kurtosis as a measure of the breadth or tailedness of the fragment size\ndistribution. The calculated kurtosis values for mRNA fragments indicated a broader distribution\nin SLiPiR-seq (11.25) compared to NEBNext (13.22). Similarly, for lncRNA fragments, SLiPiR-\nseq (14.92) exhibited a broader distribution than NEBNext (30.58).\nIn conclusion, the assessment of the fragment size distribution, incorporating both α and β\nparameters from the Gamma distribution fitting and kurtosis values, indicates that SLiPiR-seq\nproduces smoother and broader mRNA and lncRNA fragment distributions. We have incorporated\nthese refined analyses into the main manuscript as per the reviewer's recommendations to provide\na more accurate comparison of fragment size distributions.\nQs3: The authors showed that the RPM values correlates with qPCR Ct values, which supports the\nreliability of RPM normalisation. Ideally, artificial spike-ins that corresponds to the cell number or\nDNA quantity can be used to measure if a global difference in cumulative expression level exists\nin future studies.\nThe authors further mentioned the use of sequence-based exact matching and stringent mapping\nparameters to reduce the chance of multi-mapping. Since some RNA species have very similar\nsequences (e.g. miRNA), these filters may affect the sensitivity of detection. Is there any figure or\ntable that summarises the impact of these filtering strategies on the number of detected RNA\nspecies?\nAnswer: We appreciated the valuable suggestion made by the reviewer regarding the utilization of\nartificial spike-ins for more accurate quantification of SLiPiR-seq results. This recommendation\nhas been duly noted, and we plan to incorporate this experimental approach in future studies. The\nincorporation of artificial spike-ins is anticipated to provide enhanced precision in evaluating\ndetection accuracy and further enable the assessment of global differences in cumulative expression\nlevels.\nFor the mapping issues, we have implemented sequence-based exact matching strategies for rsRNA,\nysRNA, and tsRNA due to their closely related reference sequences, as detailed in table1 shown\nbelow. Notably, the reference rsRNA and ysRNA sequences were derived from the SAM file, which\nwas generated through mapping the clean reads to human rRNA or RNY genomic sequence.\nTherefore, the reference sequences were comprehensive, and thus minimize the risk of missed or\ninaccurately detected RNA species. Similarly, tsRNA data were obtained from MINTbase, known\nfor its exhaustive coverage of all tsRNA sequences, further mitigating the possibility of\nmisdetection.\nTable 1 Reference sequence of representative rsRNA, ysRNA, and tsRNA.\nID Sequence\nrsRNA：\nHomo-18S-8780 TAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATT\nHomo-18S-5733 GTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATT\nHomo-18S-4933 CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATT\nHomo-18S-4893 TCGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATT\nHomo-18S-16704 GTCGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATT\nHomo-18S-6367 AGTCGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATT\nysRNA:\nRNY4-21 CCCACTGCTAAATTTGACTGGCTT\nRNY4-9 CCCCACTGCTAAATTTGACTGGCTT\nRNY4-16 CCCCCACTGCTAAATTTGACTGGCTT\nRNY4-33 CCCCCCACTGCTAAATTTGACTGGCTT\nRNY4-177 ACCCCCCACTGCTAAATTTGACTGGCTT\nRNY4-433 AACCCCCCACTGCTAAATTTGACTGGCTT\ntsRNA:\ntRF-17-BDHPIWL AAAACATCAGATTGTGA\ntRF-18-BDHPIWD6 AAAACATCAGATTGTGAA\ntRF-19-BDHPIWH7 AAAACATCAGATTGTGAAT\ntRF-20-BDHPIWOM AAAACATCAGATTGTGAATC\ntRF-21-BDHPIWOME AAAACATCAGATTGTGAATCT\ntRF-22-BDHPIWOMP AAAACATCAGATTGTGAATCTG\nIn contrast, miRNA detection is rely on the reference sequences from miRbase, which cannot be\nused to do exact sequence-based matching due to the inherent limitation of not covering all possible\nsequences generated through sequencing. To address this, a two-step filtering process was\nimplemented in this study. The first step excluded fragments shorter than 19 nucleotides, followed\nby mapping using bowtie2 with default parameters. Subsequently, only miRNA reads exhibiting an\nexact match without insertions, deletions, or soft/hard clipping were included, ensuring accurate\nmiRNA counting. This stringent approach aimed to distinguish miRNA from tsRNA, given their\npotential sequence similarities. Specific examples, such as the similarity between has-miR-1260b\nand tRF-17-HR0VX6J, hsa-miR-7977 and tRF-17-YR66EFJ, and hsa-miR-4286 and tRF-17-\n0RER9LJ, were highlighted to underscore the necessity of the filtering step in preventing potential\nmisannotation (Table 2 shown below). The subsequent figure illustrates the impact of these filtering\nsteps on the decrease in detected miRNA species (Fig. 1a). Notably, the majority of the omitted\nmiRNA species exhibited either sequence similarity to tsRNA or lower expression levels.\nTable 2 Comparison of representative miRNA and tsRNA with similar sequences.\nID Sequence\ntRF-17-HR0VX6 ATCCCACCGCTGCCACC\nhas-miR-1260b ATCCCACCACTGCCACCAT\ntRF-17-0RER9LJ ACCCCACTCCTGGTACC\nhas-miR-4286 ACCCCACTCCTGGTACC\ntRF-17-YR66EFJ TTCCCGGCCAACGCACC\nhas-miR-7977 TTCCCAGCCAACGCACCA\nQs4: No further comments\nQs5: What are the potential reasons that might lead to the underestimation of miRNA expression\nby SLiPiR-seq? In relation to my comments in Qs3, could the stringent filtering strategies affect\nthe sensitivity of counting miRNAs?\nAnswer: We appreciate the insightful question raised by the reviewer concerning the potential\nreasons for the observed underestimation of miRNA expression in SLiPiR-seq. The\nunderestimation may indeed be attributed to the applied filtering and read count strategies, as\nexpounded in response to Qs3. To elucidate this, a comparison of the rpm values for the Top 5\ndetected miRNA, both with and without the filtering steps, was conducted (Fig. 1b). Notably, a\nslight decrease in rpm values was observed for miRNAs detected through the filtering steps (Fig.\n1b). This phenomenon may contribute to the observed underestimation of miRNA in the correlation\nanalysis of qPCR with SLiPiR-seq results, as depicted in Fig2a within the manuscript.\nA comprehensive correlation analysis was further conducted, several miRNAs with low expression\nor overlapping sequences with tsRNA were exclusively detected by the non-filtering process (Fig\n1c). These results implies that while the filtering step could enhance the precision of miRNA\nannotation, it may concurrently lead to diminished detection sensitivity, as pointed out by the\nreviewer. Despite these observations, a robust concordance persisted between results obtained\nthrough the filtering and non-filtering strategies (R=0.96, Fig. 1c). This high degree of correlation\nsuggests that, while the filtering steps might impact the sensitivity of miRNA detection, they do not\nalter the relative expression levels of individual miRNAs in plasma. Importantly, this observation\nunderscores that the integrity of our downstream differential analysis remains uncompromised by\nthe filtering processes.\nFig. 1 The impact of filtering procedure on plasma miRNA detection.\nQs6: No further comments\nQs7: No further comments\nQs8: No further comments\nQs9: Please show the number of samples in the discovery, testing, and validation cohort. Besides,\nplease indicate the cohort splitting strategy in detail, such as whether random sampling or\nstratified sampling were used.\nAnswer: There are 245 samples in the discovery cohort (196 samples for training and 49 samples\nfor testing) in the lung cancer detection study. There are 382 samples in the discovery cohort (306\nsamples for training and 76 samples for testing) in the pan-cancer classification study. There are 53\nsamples in the validation cohort for both studies. The discovery cohort was randomly partitioned\ninto training and test sets in an 80-20% manner using the createDataPartition function in the R\npackage caret. The cohort splitting strategy was based on random sampling. The sample\npartitioning and model train-test processes were repeated 100 times for all machine learning\nanalyses to avoid sampling bias.\nWe have added the number of samples in each cohort in the “Method - Sample grouping and\npartitioning” section of our manuscript: “In the lung cancer study, samples from NOR_SZBA and\nLC_SZDE were grouped as a discovery cohort (N=245). In the pan-cancer study, samples from\nNOR_SZBA, LC_SZDE, BRC_NXYK, CRC_NXYK, GC_NXYK and HCC_NXYK were\ngrouped as a discovery cohort (N=382). Samples from NOR_SZDW and LC_SZBU were grouped\nas a validation cohort (N=53) to validate models established in both the lung cancer detection study\nand the pan-cancer classification study.”"}]

Decision letter