A comparative analysis of RNA sequencing methods with ribosome RNA depletion for degraded and low-input total RNA from formalin-fixed and paraffin-embedded samples

Background Formalin-fixed and paraffin-embedded (FFPE) blocks held in clinical laboratories are an invaluable resource for clinical research, especially in the era of personalized medicine. It is important to accurately quantitate gene expression with degraded and small amounts of total RNA from FFPE materials. Results High concordance in transcript quantifications were shown between FF and FFPE samples using the same kit. The gene expression using the TaKaRa kit showed a difference with other kits, which may be due to the different principle of rRNA depletion or the amount of input total RNA. For seriously degraded RNA from FFPE samples, libraries could be constructed with as low as 50 ng of total RNA, although there was residual rRNA in the libraries. Data analysis with HISAT demonstrated that the unique mapping ratio, percentage of exons in unique mapping reads and number of detected genes decreased along with the decreasing quality of input RNA. Conclusions The method of RNA library construction with rRNA depletion can be used for clinical FFPE samples. For degraded and low-input RNA samples, it is still possible to obtain repeatable RNA expression profiling but with a low unique mapping ratio and high residual rRNA.


Background
With the development of massive parallel sequencing, RNA-Seq has become an useful tool for transcriptome analysis, as well as for the identification of novel transcripts, SNPs, gene fusion and alternative splicing events [1]. Formalin-fixed and paraffin-embedded (FFPE) blocks held in clinical laboratories are an invaluable resource for clinical research, especially in the era of personalized medicine. FFPE samples are easy to store, preserve tissue morphology for clinical and pathological observation, and preserve nucleic acids for molecular biology research [2]. Currently, many clinical tests are based on the expression of certain genes, such as the Mamma-Print test, to assess recurrence risk in early-stage breast cancer [3] and the tissue of origin (TOO) test to find the site of the primary tumor. In addition, RNA expression profiles have become an important source of new biomarkers with potential values in cancer metastasis and disease prognosis [4,5]. The discovery and development of these diagnostic and prognostic biomarkers will rely heavily on retrospective studies on historical FFPE samples [6]. Therefore, it is important to accurately quantitate the gene expression with total RNA from FFPE materials.
RNA-seq requires the enrichment of mature mRNAs, or the depletion of highly abundant ribosomal RNAs (rRNAs) from total RNA before sequencing. RNAs from FFPE materials are usually degraded to small sizes without the 3′poly (A) tail; moreover, recent studies suggest that certain functionally important mRNAs are non-poly (A) RNAs [7]. Therefore, capturing the 3′poly (A) tail is not a compatible method, especially when the starting materials are from FFPE samples. Another method for RNA-seq of FFPE samples is cDNA hybrid capture using a whole exome DNA probe to hybridize to the total RNA library. The yield of on-exon data was increased significantly due to the cDNA-capture, while the accuracy of quantitated gene expression was decreased [8,9]. The signals of low gene expression might be missed by decreased uniformity of the exome probe.
For RNA-seq of FFPE samples, rRNA depletion from total RNA is the optimal solution. Nucleic acids extracted from FFPE blocks are fragmented and chemically modified, making them controversial to use in molecular diagnosis. rRNA depletion protocols could keep as much information as possible from the total RNA. There are several rRNA depletion protocols. The first method that is commonly used hybridizes the rRNA to a DNA probe and degrades the rRNA: DNA hybrids using RNase H. In the second method, rRNA is captured by complementary DNAs, which are coupled to paramagnetic beads, and the mixture is removed from the reaction [10]. Several studies have shown that FFPE RNA-seq data produced high concordance with RNA-seq results from matched frozen fresh samples [11,12]. Previous studies have confirmed that for low-quality RNA, especially for degraded FFPE RNA, the RNase H method performed best [13]. The third method, which is suitable for low-input and low-quality samples, first transcribes total RNA to cDNA, and then the ZapR enzyme digests all rRNA: DNA hybrids. With an increasing number of commercially available RNA library preparation kits based on the principle of rRNA removal, we can make the best use of clinical FFPE samples. All those kits utilizing these methods are available, but the effect of the efficiency of rRNA removal on RNA-seq data is still unclear.
In this study, we compared four FFPE RNA library preparation kits (KAPA, TaKaRa, QIAGEN and Vazyme) based on two principles of rRNA depletion, with degraded RNA from FFPE samples and paired FF samples as starting materials (Fig. 1). Takara Kit only requires input of 5 to 50 ng total RNA with chemical modifications, such as those extracted from FFPE tissue and input of 250 pg to 10 ng total RNA for FF samples. After total RNA was fragmented or denatured, cDNA was synthesized, including cDNA from rRNA. In the next step, the synthesis of cDNA was added full-length Illumina adapters by a first round of PCR amplification (PCR1), including barcodes. And then, originating from rRNA of the ribosomal cDNA was cleaved by ZapR in the presence of the R-Probes. Finally, untouched and originating from non-rRNA molecules were enriched by a second round of PCR amplification (PCR2), and purified the final library. KAPA kit has been validated for library construction from 25 ng to 1 μg of total RNA. This kit using Oligo Hybridization and rRNA Depletion eliminated the effect of ribosomal RNA on library. The rRNA duplexed to DNA oligos was digested by RNase H treatment. Before the cDNA synthesis, hybridization oligos were removed from the sample by DNase I digestion. The rRNAdepleted RNA is eluted and fragmented to the desired size using high temperature in the presence of Mg 2+ . And then, 1st strand and 2nd strand cDNA was synthesized successively, of which 2nd strand cDNA was marked by dUTP. The dAMP was then added to the 3′-end of dscDNA fragments, and 3′-dTMP adapters are ligated to 3′-dAMP library fragments. After fragment separation, PCR amplification was performed on the final library.
Vazyme kit is mainly applicable to the total RNA of human, mouse and rat with a starting value of 0.1-1 μg, and also applicable to the construction of the library for the degradation of RNA samples of the above species. QIAGEN Kit need 1-100 ng enriched, poly(A) + RNA. So we used the first few steps of Vaths™ Total RNA-seq (H/ M/R) Library Prep Kit protocol to get the poly(A) + RNA. The removal of ribosomal RNA from both Vazyme and QIAGEN kits was similar to KAPA kit.
In addition, we evaluated the effect of bioanalysis tools on the total mapping rate, unique mapping rate, exon percentage and number of detected genes using FF samples and FFPE samples. HISAT (hierarchical indexing for spliced alignment of transcripts) allows scientists to align reads to a genome, assemble transcripts, compute the abundance of these transcripts in each sample and compare experiments to identify differentially expressed genes and transcripts [14]. STAR (Spliced Transcripts Alignment to a Reference) can discover noncanonical splices and chimeric (fusion) transcripts and is also capable of mapping full-length RNA sequences [15]. STAR generates output files that can be used for many downstream analyses, such as transcript/ gene expression quantification, differential gene expression, novel isoform reconstruction, signal visualization, and so forth [16]. Both tools are free, open-source methods for comprehensive analysis of RNA-seq experiments.
In the last part of this study we evaluated the performance of two kits allowing for lower input of total RNA because many clinical studies need to use RNA, even though a low quality and a very low input of RNA can be extracted from clinical FFPE samples. We also validated the reproducibility of low-quality and low-quantity samples.

Performance of four RNA-seq preparation kits for FF and FFPE samples
To evaluate the performance of four RNA-seq preparation kits, we collected total RNA from GM12878 FF and FFPE samples. The quality of the two RNA samples is shown in Additional file 1: Figure S1. We constructed RNA-seq libraries following the recommended protocols respectively. After sequencing, the raw data of all eight libraries were down sampled to 18 G and analytical comparisons were focused on several fields including the yield of libraries, GC content, rRNA depletion efficiency, genome alignment profiles, transcriptome coverage, transcript quantification, etc. (Table 1).
The recommended input is even lower for the TaKaRa kit than the other three kits, so we input 10 ng of total RNA for preparing the library, while the input of the other kits was 100 ng. The library yields and exon percent in the unique mapping data of the FFPE sample with the TaKaRa kit was the highest (Table 1 and Figure 2), which indicated that the TaKaRa kit is intended for low-input starting material. The performance of the other three kits showed a similar tendency of the library yields and exon percentage in the unique mapping data of the FFPE samples being much lower than that of the FF samples. Residual rRNA in the TaKaRa library was also the highest and had the least clean data, which was due to the removal of ribosomal cDNA (cDNA fragments originating from rRNA molecules) after cDNA synthesis using probes specific to mammalian rRNA.
As shown in Figure 3, the total number of genes detected from the FFPE samples was similar among the four libraries. The number of genes detected in the TaKaRa library of the FF sample was more than twice as much as detected in the other libraries, even with using less input total RNA. We also used sample 13, sample 14 and sample 15 which were from native external quality assessment samples to test the four RNA-seq library preparation kits. As shown in Additional file 1: Table S1, we got the similar results to FFPE sample of GM12878.
RNA-seq is an established platform for quantifying gene expression using high-quality RNA. To evaluate the gene expression performance of the FF and FFPE samples across the four kits, we compared the consistency of transcript quantification from matched pairs of FF and FFPE samples using the same kit ( Figure 4). The results showed high concordance in transcript quantifications between FF and FFPE samples using the same kit (R (FF vs FFPE) = 0.96 for the TaKaRa kit, R (FF vs FFPE) = 0.97 for the Vazyme and QIAGEN kits, R (FF vs FFPE) = 0.98 for the KAPA kit). In addition, we compared the consistency of FF (or FFPE) samples between different kits. The consistency among the KAPA, Vazyme and QIA-GEN kits was higher than that of the four kits. Among the four kits, KAPA and QIAGEN showed the highest consistency, not only for    Table S2. To clarify the difference between the TaKaRa kit and any one of the other three kits in FFPE samples and FF samples, we chose the differential transcripts, which had more than a 50-fold difference. There were a total of 37 differential transcripts in the FF sample and 58 differential transcripts in the FFPE sample (Additional file 1: Table S3). There were 16 differential transcripts found both in the FF sample and in the FFPE sample. Most of these differential transcripts were mitochondrially encoded RNA, small nucleolar RNA, and 5S ribosomal pseudogene, all of which were noncoding RNA. Only a few transcripts were from coding RNA, such as the PET117 homolog, Karyopherin subunit alpha 7, and BolA family member 2B. The FPKMs of these transcripts in TaKaRa libraries were higher than those in other libraries, but not more than 10. These results indicate that the main difference between the TaKaRa libraries and the other three libraries was caused by noncoding residual RNA, and for the quantification of transcripts from coding RNA, there was no significant difference among the four RNA-seq libraries.

Comparison of two bioanalysis methods with FF and FFPE samples
We evaluated the effect of bioanalysis tools on the total mapping rate, unique mapping rate, exon percentage and number of detected genes using FF samples and FFPE samples. For all the samples, there was almost no differences between HISAT and STAR on the quality data (Additional file 1: Table S4) regardless of RNA-seq preparation kits. Due to time and computer space, we used the HISAT analysis method to analyze data in our assay.

RNA-seq library kit for degraded and lower input of total RNA from FFPE samples
Many clinical studies, such as fusion detection, gene expression profiling, identification of novel transcripts and detection of alternative spicing events, want to use RNA, even though a low quality and a very low input of RNA can be extracted from clinical FFPE samples. To meet this need, we tested two kits allowing for a lower input of total RNA. The detailed results are shown in Table 2. We used the recommended cycles for each kit and obtained significantly higher library yields from the TaKaRa kit than from the KAPA kit. When raw data were down-sampled to 20 G, fewer clean data were left in the TaKaRa library because there were more reads from rRNA in its library. Although the total mapping rate in the TaKaRa library was also lower than it was in the KAPA library, exon % in the TaKaRa library was higher. A similar number of genes were detected by both kits. The correlations of transcript quantification between the two inputs and two kits are shown in Figure 5. This result demonstrated that the performance of the TaKaRa kit may be sufficient when the total RNA input is as low as 10 ng, which may be more compatible for use with RNA coming from valuable FFPE samples while reducing the depletion of samples.

Performance of two kits with different quality of input total RNA
Another serious problem for use of clinical FFPE samples is low quality. The Agilent RNA Integrity Number (RIN) of most FFPE samples was so poor that it was not sensitive enough to evaluate the quality of RNA from degraded FFPE samples. Here, we used the reference of DV200%, the percentage of RNA fragments > 200 nucleotides, to assess FFPE RNA quality. We tested the two kits with 15 different qualities of FFPE RNA samples (Additional file 1: Figure S2). The total RNA input was 50 ng for all the samples, and the recommended PCR cycles were used for each kit. As shown in Table 3, the KAPA kit failed to construct a library for some poor quality RNA samples, or the library was insufficient to obtain more data, while all the TaKaRa libraries were successfully constructed and sequenced. Moreover, more transcripts were detected from the TaKaRa libraries than from the KAPA libraries. Similar to previous results, for all the samples when the raw data were down-sampled, fewer data were left in the TaKaRa library because residual rRNA in the TaKaRa library was much more than that of the KAPA library. The worse the quality of RNA is, the lower the percentage of exons in unique mapping reads.
To test the reproducibility of the TaKaRa kit with low quality samples, we repeated the RNA library of five FFPE samples (sample 22 to 27 except sample 26 due to insufficient total RNA). The reproducibility performance of five low-quality clinical samples was shown in Table 3. As shown in Figure 6, the results showed high concordance (R > 0. 8          two batches. The reproducibility may be related to low quality of input total RNA.

Discussion
RNA-seq of clinical FFPE samples could provide more important and reliable information for discovery and validation of biomarkers. Previous research [17,18] and the results of this study also showed that FFPE RNA-seq provided reliable gene expression data, comparable to that obtained from fresh frozen tissue with the method of rRNA depletion. Standard practice in tissue fixing and paraffin-embedding has little impact on the expression analysis of RNA samples, which makes archived FFPE samples valuable and retrospective studies feasible. However, there is a limitation in the study in that we used a freshly cultured cell line and newly prepared cell FFPE blocks, and we did not compare FFPE samples with longer storage. The difference in the principle of rRNA depletion could still result in a difference in library yield, residual rRNA, the percent of exons in unique mapping reads and transcript quantification. The KAPA kit, QIAGEN kit and Vazyme kit used the same principle of rRNA depletion, so a high concordance of transcript quantification was shown among the three kits. These three kits removed rRNA from total RNA before cDNA synthesis, using a rRNA probe to combine rRNA and then digesting rRNA by RNase H and removing the rRNA probe by DNase I. The library of the TaKaRa kit showed a difference with higher yield, residual rRNA, exon percentage and the number of detected genes, using a lower RNA input compared to the other three kits, which might result from the unique method of rRNA depletion. The workflow used in the TaKaRa kit takes advantage of a novel technology allowing removal of ribosomal cDNA (cDNA fragments originating from rRNA molecules) after cDNA synthesis using probes specific to mammalian rRNA. The specificity and number of probes could have an effect on the rRNA depletion efficiency, especially in low-quality RNA of FFPE samples, so there was higher residual rRNA in the TaKaRa library. In our results, we found that there were 9 transcripts (ENSG00000201998.  1000) in the TaKaRa libraries, while their FPKM were very low in the other kits. These transcripts included mitochondrially encoded 16S RNA, mitochondrially encoded 12S RNA, small nucleolar RNA, the 5S ribosomal pseudogene, and 5S ribosomal 9, all of which were noncoding RNAs. These transcripts detected in the TaKaRa libraries were due to low efficiency of rRNA depletion. For the same reason, we only obtained fewer clean data with the TaKaRa kit compared to the KAPA kit when using poorquality RNA to construct libraries. The shortcoming of high residual rRNA would waste sequencing reads and increase the cost of RNA-seq.
On the other hand, the strategy of rRNA depletion after cDNA synthesis made the TaKaRa kit especially well-suited for working with very small quantities of total RNA. A similar strategy was adopted not only by TaKaRa kit but also by the Nugen kit [19]. We tried to decrease the input of total RNA, and both the TaKaRa kit and the KAPA kit showed good concordance with higher RNA input. For the TaKaRa kit, the concordance of two libraries (50 ng vs. 10 ng) was 0.99. For seriously degraded FFPE RNA, a RNA-seq library was still successfully constructed and repeated by the TaKaRa kit, but not the KAPA kit, which indicated that initial rRNA depletion from total RNA was not very effective and often leaves an insufficient amount of material for preparation of high-quality libraries.

Conclusions
The concordance between FF and FFPE samples is excellent for any of four RNA-seq library kits. Therefore, FFPE could be used for the RNA-seq profiling with the methods of rRNA removal from total RNA. The difference between TaKaRa and the other three kits for FF and FFPE samples might be due to the different principle of rRNA removal or different input of total RNA. Both the KAPA and TaKaRa kits allowed low total RNA input and consistent transcript quantification was obtained between the lowest input and a higher input. When the quality of input total RNA was high, in which the DV 200% was more than 30%, both the KAPA and TaKaRa kits performed well. When the DV 200% of degraded RNA was less than 30%, lower quality data with lower unique mapping and lower exon percentage or failure of library construction will be evident. The TaKaRa RNA-seq library kits could be used for RNA-seq library construction of low-quality and low-quantity FFPE samples. Although the rRNA residual is a little higher, it could detect more transcripts and showed good reproducibility with low-quality and low-quantity FFPE samples.

GM12878 fresh cell and the preparation of cell FFPE
GM12878 cell line, which was originated from human B-lymphoblastoid cells and now often used as control sample in NGS, was obtained from Cobioer Biological Technology (Nanjing, China) and cultured using the recommended culture conditions. Briefly, GM12878 was incubated at 37°C in an incubator (Haier, China) with 5% CO 2 in air atmosphere, and culture media RPMI 1640 (Gibco, Cat#12633-020) contained 10% fetal bovine serum (FBS) (Gibco, Cat#10437028). About 1 × 10 7 cells were treated with 10% neutral formalin for 1 h and then the fixed cells were collected by centrifugation of 3000 rpm for 10 min (Eppendorf, 5810R). The cells were suspended with 0.5 ml 1xPBS and then was added into 2% agarose gel. The mixture was cooled and solidified into block. The GM12878 cell block was processed according to standard formalin-fixed and paraffin-embedded (FFPE) methods [20].

Clinical samples and ethics
There were 5 FF clinical samples (sample 1 to 5) and 22 FFPE samples (sample 6 to 27). FFPE sample 13, 14 and 15 were samples from native external quality assessment. All the 24 clinical samples (sample 1 to 27 except sample 13,14 and 15) were collected from the Cancer Hospital Chinese Academy of Medical Sciences & Peking Union Medical College. FFPE tissue slides were examined by expert pathologists including a minimum of 20% cancer cells. All these cancer samples from patients who signed the informed consent forms and were allowed to be used in other researches (Additional file 1: Table S5).

RNA isolation and assessment of quality
GM12878 fresh cell and 5 FF samples were isolated total RNA with column purification of AllPrep DNA/RNA Mini Kit (QIAGEN, Cat#80204), and three 5um-sections of the FFPE samples were used to extract total RNA using AllPrep DNA/RNA FFPE Kit (QIAGEN, Cat#80234), both according to the manufacturer's recommendations. We could get enough RNAs from above samples just following the detailed user guide of the kits.
The final RNA concentration was typically measured by Qubit RNA HS assay kit (Thermo Fisher Scientific, Cat#Q32855). The integrity of RNA was determined by the RNA integrity number (RIN) and DV200 (percentage of RNA fragments greater than 200 nt) with Eukaryote total RNA pico 6000 Assay of the 2100 Bioanalyzer (Agilent). We constructed libraries strictly according to the user guide of each kit and libraries prepared with the above four kits were sequenced using 150-bp paired-end runs on Illumina NGS systems (HiSeq® 2500 and Xten) after quantification by the Qubit dsDNA Assay Kit and determination of fragment length by the Agilent 2100 Bioanalyzer with the DNA 1000 Kit.

Bioinformatic analysis
Raw data (raw reads) of fastq format were firstly processed through in-house scripts. In this step, raw reads, Q20, Q30 and GC content were calculated. Then, Trimmomatic v0.36 [21] was used to trim reads containing adapter and low quality reads from raw data. Clean data (clean reads) were obtained by removing reads mapping to rRNA reference genome using bowtie2. All the downstream analyses were based on the clean data with high quality. GENCODE GRCh37 (version 19) reference genome and gene model annotation files were downloaded from genome website directly. Then HISAT (hierarchical indexing for spliced alignment of transcripts) and STAR (Spliced Transcripts Alignment to a Reference) were both used to analyze the data of FF samples and FFPE samples. Index of the reference genome was built using HISAT2 v2.1.0 [22] or STAR v2.6.0c [15] and pairedend clean reads were aligned to the reference genome using HISAT2 v2.1.0 or STAR v2.6.0c. RSeQC v2.6.4 [23] was used to calculate how mapped reads were distributed over genome feature (exon, intron and intergenic), and nucleotide composition for each position of read. FeatureCounts v1.6.1 [24] or HTSeq v0.10.0 [25] was used to count the reads numbers mapped to each gene. FPKM, expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced, considers the effect of sequencing depth and gene length for the reads count at the same time, and is currently the most commonly used method for estimating gene expression levels. And then FPKM of each gene was calculated based on the length of the gene and reads count mapped to this gene. FPKM was used to determine concordance between each kits by Pearson correlation coefficient (R).