Swarm Intelligence-Enhanced Detection of Non-Small-Cell Lung Cancer Using Tumor-Educated Platelets

Abstract

Blood-based liquid biopsies, including tumor-educated blood platelets (TEPs), have emerged as promising biomarker sources for non-invasive detection of cancer. Here we demonstrate that particle-swarm optimization (PSO)-enhanced algorithms enable efficient selection of RNA biomarker panels from platelet RNA-sequencing libraries (n = 779). This resulted in accurate TEP-based detection of early- and late-stage non-small-cell lung cancer (n = 518 late-stage validation cohort, accuracy, 88%; AUC, 0.94; 95% CI, 0.92-0.96; p < 0.001; n = 106 early-stage validation cohort, accuracy, 81%; AUC, 0.89; 95% CI, 0.83-0.95; p < 0.001), independent of age of the individuals, smoking habits, whole-blood storage time, and various inflammatory conditions. PSO enabled selection of gene panels to diagnose cancer from TEPs, suggesting that swarm intelligence may also benefit the optimization of diagnostics readout of other liquid biopsy biosources.

Keywords: NSCLC; RNA; blood platelets; cancer diagnostics; liquid biopsies; particle-swarm optimization; self-learning algorithms; splicing; swarm intelligence; tumor-educated platelets.

Graphical abstract

Figure 1

Technical Performance of thromboSeq (A) Platelet total RNA yield in ng/μL isolated from 6 mL whole blood in EDTA-coated Vacutainers tubes. p Value calculated by independent Student's t test. (B) Correlation plot (Pearson's correlation) of estimated RNA input to the output SMARTer cDNA yield. Each dot represents a sample, color-coded by clinical group. (C) Correlation plot (Pearson's correlation) of SMARTer cDNA yield to the Truseq cDNA library yield. Each dot represents a sample, color-coded by clinical group. (D) Bioanalyzer profiles of samples with spiked, smooth, and intermediate spiked/smooth traces for total RNA as measured by RNA 6000 Picochip (left column), SMARTer amplified cDNA as measured by DNA High Sensitivity chip (middle column), and Truseq cDNA libraries as measured by DNA 7500 chip (right column). (E) Shown are for spiked, smooth, and intermediate spiked/smooth SMARTer cDNA trace profiles for both non-cancer and NSCLC are the relative cDNA yield in nmol following SMARTer amplification (top), average cDNA length per 1,000 bp following SMARTer amplification (middle), and the number of intron-spanning spliced RNA reads (bottom). (F) Selection of intron-spanning spliced RNA reads for thromboSeq analysis. Stackplot indicates the distribution of reads for each sample, subspecified from intron-spanning (orange), exonic (yellow), intronic (green), intergenic (blue), and mitochondrial (purple) regions. (G) Selection of samples with >750 genes detected for thromboSeq analysis. (H) Leave-one-sample-out cross-correlation filtering step for which counts of each sample were correlated to the median counts of all other samples. (I) Number of genes detected with confidence in the platelet RNA samples using shallow thromboSeq. (J) Deep thromboSeq versus shallow thromboSeq for matched platelet samples. Median total read counts (min-max) in millions: 59.7 (43.2–96.2) for deep thromboSeq and 12.9 (11.6–20.0) for shallow thromboSeq. The three genes with highest expression in deep thromboSeq are highlighted. The boxes of the boxplots indicate the interquartile range (IQR), the horizontal black line indicates the median values, and the whiskers range 1.5× the IQR. See also Figure S1 and Table S1.

Figure 2

Analysis of the Spliced RNA Repertoire of TEPs from NSCLC Patients (A) Schematic figure represents the read distribution analyses procedure. A total of 100 bp reads were mapped to the human genome and reads mapping to four distinct regions were quantified. mtDNA, mitochondrial genome. (B) Boxplots indicate for non-cancer (green, n = 104) and NSCLC (blue, n = 159) the median and spread of reads mapping to mitochondrial (mtDNA), exonic, intronic, or intergenic regions, and the median and spread of intron-spanning and exon boundary reads (normalized to one million total genomic reads). The boxes of the boxplots indicate the IQR, the horizontal black line indicates the median values, and the whiskers range 1.5× the IQR. (C) Unsupervised hierarchical clustering of differentially spliced RNAs between non-cancer (n = 104) and NSCLC (n = 159) individuals, with FDR < 0.01. Columns indicate samples, rows indicate genes, and color intensity represents the Z score-transformed RNA expression values. Clustering of samples showed non-random partitioning (p < 0.0001, n = 263, Fisher's exact test). (D) PAGODA GO analysis. Most significant results by adjusted Z score, indicating high statistical significance, were clustered and visualized. Color code indicates a green-to-orange (low-to-high) score per sample per gene cluster. See also Table S2.

Figure 3

Alternative Splicing and Exon Skipping in TEPs of Patients with NSCLC (A) Schematic figure represents the development of an isoform count matrix that contains, per sample for each expressed RNA isoform, the number of reads supporting that particular isoform. RNA isoforms were inferred from the RNA-seq data using MISO. The isoform count matrix is employed for non-cancer versus NSCLC differential splicing analysis. The left pie chart indicates the gene loci (in total n = 571) with one or more differentially spliced isoforms (encoded in the color bars). Two gene loci with 10 isoforms each were not indicated. The right pie chart indicates the number of gene loci with multiple differentially spliced isoforms that show concordant up (blue box), concordant down (red box), or both directions (alternative isoforms; green box). (B) Schematic figure represents the MISO algorithm mapping exon-skipping events, thereby calculating the percent spliced in (PSI) value. Subsequently, ANOVA statistics were applied to individual exons with sufficient coverage in the dataset (n = 230) comparing read levels in non-cancer and NSCLC individuals, and assessed for overall in- or exclusion by calculation of the ΔPSI value. Histogram shows the direction of the PSI value for 27 exons with FDR < 0.01 (positive PSI values favors inclusion in NSCLC, whereas negative PSI values favors exclusion in NSCLC). Individuals gene names associated with the 27 exons are listed in the box.

Figure 4

P-Selectin Signature (A) Correlation plot (Pearson's correlation) of proportion of reads mapping to exonic coordinates (x axis) versus the log₂-transformed, RUV-corrected, and counts-per-million (logCPM) of P-selectin. Each dot represents a sample, colored by clinical group. (B) Distribution of Pearson's correlation coefficients of the correlation between logCPM levels of 4,722 genes and the logCPM of P-selectin. A subset of the genes show a strong positive correlation with P-selectin (r approximates 1), whereas others do not (r approximates 0). (C) Venn diagram overlay of genes upregulated in the NSCLC TEP signature and genes with a significant positive Pearson's correlation (FDR < 0.01) toward P-selectin. Number of overlapping genes is indicated in the Venn diagram.

Figure 5

RNA-Binding Protein Analysis of TEP-Derived RNA Signatures (A) Schematic biological model highlighting the difference between nucleated cells and anucleated platelets in the context of regulation of translation. Nucleated cells (left) are able to regulate and maintain the transcriptome by transcription factor (TF)-mediated DNA transcription, in contrast to platelets. (B) Schematic representation of the RNA-binding protein (RBP)-thromboSearch engine algorithm. Reference genome sequences from 4,722 platelet genes were matched with 547 motif sequences deconvoluted from 102 RBPs. (C) UTR-read coverage filter. Blue dots represent mean counts across all samples, and gray shading indicates the respective standard deviations. (D) Boxplots indicate the average RBP RNA expression levels (n = 43 RBPs identified) in log₂-transformed counts-per-million (logCPM) sorted by median normalized intron-spanning read level. The boxes of the boxplots indicate the IQR, the horizontal black line indicates the median values, and the whiskers range 1.5× the IQR. (E) Heatmap with unsupervised clustering of RNA levels from 43 RBPs in the age, smoking status, and blood storage time-matched cohort (n = 263 samples; n = 104 non-cancer, n = 159 NSCLC). The RBP RNA levels enable non-random clustering of the samples (p < 0.0001, n = 263, Fisher's exact test). (F) Enrichment of identified RBP binding sites per UTR region. The x and y axes represent the mean binding sites for the 5′ and 3′ UTR per RBP (dots, n = 80). Several RBPs are specifically enriched in the 3′ UTR, whereas others are enriched in the 5′ UTR. (G and H) Heatmap of all RBPs (rows) and all 5′ UTR (G) and 3′ UTR (H) regions detected with sufficient coverage in platelets (column). Number of binding sites is indicated in the color coding. (I) Spearman's rank correlation analysis between n binding sites of an RBP and the logarithmic fold-change (logFC) of genes in the NSCLC/non-cancer differential splicing analysis. Positive correlations indicate an enrichment in binding sites with an increase of the logFC, whereas negative correlations indicate the opposite. Plots indicate the relation between the Spearman's correlation coefficient (x axis) and the concomitant differential splicing FDR value.

Figure 6

PSO-Enhanced thromboSeq for NSCLC Diagnostics (A) Schematic representation of the particle-swarm optimization (PSO) approach. Yellow-to-red colored dots represent AUC values of a cohort classified using a thromboSeq classification algorithm, with use of 100 randomly compiled biomarker gene panels (left) or 100 biomarker gene panels proposed by swarm intelligence (right). Dots were mirrored twice for visualization purposes. Most optimal AUC value reached by PSO-enhanced thromboSeq is indicated in both plots with an asterisk. (B) ROC analysis of PSO-enhanced thromboSeq classifications using only the matched non-cancer and NSCLC cohorts. Gray dashed line indicates ROC evaluation of the training cohort assessed by LOOCV, red line indicates ROC evaluation of the evaluation cohort (n = 40), blue line indicates ROC evaluation of the validation cohort (n = 130). Indicated are cohort size, most optimal accuracy, and AUC value. (C) ROC analysis of PSO-enhanced thromboSeq classifications using matched training and evaluation of non-cancer and NSCLC cohorts and validation in the remaining matched and unmatched samples. Gray dashed line indicates ROC evaluation of the training cohort assessed by LOOCV, red line indicates ROC evaluation of the evaluation cohort (n = 88). ROC evaluation of the late-stage NSCLC (Validation Late-st.; n = 518) and locally advanced NSCLC (Validation Loc.-adv.; n = 106) is plotted in the blue and green lines, respectively. Indicated are cohort size, most optimal accuracy, and AUC value. Acc., accuracy. See also Figures S3 and S4.