Plasma extracellular RNA profiles in healthy and cancer patients

Abstract

Extracellular vesicles are selectively enriched in RNA that has potential as disease biomarkers. To systemically characterize circulating extracellular RNA (exRNA) profiles, we performed RNA sequencing analysis on plasma extracellular vesicles derived from 50 healthy individuals and 142 cancer patients. Of ~12.6 million raw reads for each individual, the number of mappable reads aligned to RNA references was ~5.4 million including miRNAs (~40.4%), piwiRNAs (~40.0%), pseudo-genes (~3.7%), lncRNAs (~2.4%), tRNAs (~2.1%), and mRNAs (~2.1%). By expression stability testing, we identified a set of miRNAs showing relatively consistent expression, which may serve as reference control for exRNA quantification. By performing multivariate analysis of covariance, we identified significant associations of these exRNAs with age, sex and different types of cancers. In particular, down-regulation of miR-125a-5p and miR-1343-3p showed an association with all cancer types tested (false discovery rate <0.05). We developed multivariate statistical models to predict cancer status with an area under the curve from 0.68 to 0.92 depending cancer type and staging. This is the largest RNA-seq study to date for profiling exRNA species, which has not only provided a baseline reference profile for circulating exRNA, but also revealed a set of RNA candidates for reference controls and disease biomarkers.

Figure 1. Quality control of sequencing libraries.

(A) Number of raw reads vs. number of mapped reads. The horizontal and vertical lines were the median levels of raw reads and mapped reads, respectively. (B) Size distribution of sequencing library inserts. The theoretical gel-size selection region is shown in grey. (C) Mature miRNA saturation analysis using raw read counts. (D) piwiRNA saturation analysis using raw read counts. The solid red lines in (C,D) were modeled by non-linear regression.

Figure 2. Technical effects on RNA expression profiles.

(A) Principal component analysis showing effect of RNA isolation dates on RNA profile changes. (B) Principal component analysis showing effect of gel size selection dates on RNA profile changes. (C) Effect of gel size selection on the percentage distribution of miRNAs (precursor and mature miRNAs) and piwiRNAs. Batches of 40 polyacrylamide gels are shown. The upper and lower horizontal lines are the median levels of miRNAs and piwiRNAs, respectively.

Figure 3. Statistical summary of RNA species detected by RNA-seq across 192 libraries.

(A) Percentage of each RNA species in all mapped RNAs. (B) Number of mapped unique RNA references in different RNA species. (C) Percentage (sorted from high to low) of each detected mature miRNAs in all mapped miRNA reads. The top 32 miRNAs (quantile 0.98) are highlighted in red and are also shown in embedded graph. (D) Percentage (sorted from high to low) of each detected piwiRNAs in all mapped reads. Top 10 piwiRNAs (quantile 0.98) are highlighted in red and also shown in embedded graph. RNA transcripts with >4 RPM were used. (E) The top 32 miRNAs and their isoform proportions. The color bars represent the percentage of each miRNA isoform.

Figure 4. exRNA stability testing.

(A) Ten miRNAs with the highest stability. Smaller ranking values had higher stability. (B) Boxplot showing abundance level of each corresponding miRNA.

Figure 5. Association of exRNAs with biological and clinical characteristics.

(A) Boxplot of sex-associated small RNAs in 50 healthy samples. (B) Stage-dependent analysis of differentially expressed RNA transcripts in colon cancer. Number of colon cancer-associated RNA transcripts increased with disease progression. (C) Association of selected miRNAs with cancer types and clinical stages. Abundance levels of miR-125a-5p and miR-1343-3p were constantly lower in all cancer types and stages. S1-4: stages I–IV. HS: hormone sensitive prostate cancer. CR: castration resistant prostate cancer. Pan: pancreatic cancer. *FDR > 0.05. **FDR < 0.05. (D) Colon cancer statistical model showing AUC values from 0.68 to 0.81 in clinical stages I–IV. (E) Prostate cancer model demonstrating higher AUC values from 0.89 in HSPC to 0.92 in CRPC, likely contributed by more advanced stages in prostate cancer patients.