NSCLC patients exhibited a unique miRNA profile compared to healthy population
In this retrospective study, 37 plasma samples were subjected for exosome purification and miRNA-sEq. Among them, 23 samples were from 16 phase IV patients with bone metastasis (BM+) and 14 samples from 14 phase IV patients without bone metastasis (BM-). Raw reads of miRNA-seq from plasma samples were normalized to counts per million (CPM) and 1287 miRNAs were retained as expressed miRNAs in plasma exosomes (Supplementary Table 1). Using exact test, 91 miRNAs were identified as differentially expressed miRNAs (DE miRNAs) between H and C including 71 up-regulated miRNAs and 20 down-regulated miRNAs in patient group (Supplementary Table 2). Based on the DE exosomal miRNA profile, samples of healthy individuals (H) could be separated from samples of NSCLC patients (C) using supervised hierarchy clustering (Supplementary Figure S1). All these data proved that NSCLC patients exhibited a unique miRNA profile compared to healthy population.
Detection of co-expression clusters in exosomal miRNAs
To characterize the correlation pattern and predict the function of miRNAs, we applied weighed gene co-expression network analysis (WGCNA) to the 1287 miRNAs detected in the plasma exosomes. WGCNA is widely used in genomic data analysis which can detect clusters of highly correlated genes based on pairwise correlations [23]. As shown in Supplementary Figure S2, we identified three clusters of co-expressed miRNAs which were represented by different color codes (Brown: cluster A; Turquoise: cluster B; Blue: cluster C). Cluster B had 144 miRNAs which is the largest cluster among three. 49 miRNAs were assorted into Cluster A while 95 miRNAs were in Cluster C (Supplementary Table 3).
Furthermore, BM + group showed significant up-regulation in cluster B eigengene value compared to healthy population and BM- group which suggested miRNAs in cluster B might related to the initial of bone metastasis (Fig. 1B). The 144 miRNAs from cluster B basically differentiated BM + group from BM- group using unsupervised clustering which also shed light on the function of the miRNA cluster in bone metastasis (Fig. 2A). We performed pathway analysis of the 144 miRNAs in cluster B using miRNA enrichment analysis and annotation (MiEAA)[24] and top 20 miRNA pathways were shown in Table 1. Interestingly, cluster B was enriched in metabolism processes such as pyruvate metabolism (p = 0.010575), glycolysis and gluconeogenesis (p = 0.014983), purine metabolism (p = 0.014983), propanote metabolism (p = 0.014983) and pyrimidine metabolism (p = 0.027086) (Supplementary Table 4).
In cluster A, BM + group exhibited a trend of increase in eigengene value compared to healthy population and BM- group (Fig. 1A). miRNA enrichment analysis was done with 49 miRNAs in cluster A, however, not many related pathway was found. Nevertheless, we discovered that 43 out of 49 cluster A miRNAs were expressed by chromosome 14 (Supplementary Table 5).
In cluster C, the eigengene values were comparable within healthy population group, BM + and BM- group (Fig. 1C) which indicated that this cluster was not related to bone metastasis. The cluster A miRNAs failed to differ BM + from BM- with unsupervised clustering in Fig. 2B also suggested that cluster C miRNAs were not involve in the bone metastasis. Related pathways of cluster C were shown in Supplementary Table 6. miRNAs in cluster C seemed to related to all kinds of signaling pathways including cell cycle (p = 0.000322), proteasome and lysosome (p = 0.000322), p53 signaling pathway (p = 0.00037), insulin signaling pathway (p = 0.000322) and Ras pathway(p = 0.00037).
Identified differentially expressed miRNAs between BM + and BM- group as potential biomarker for bone metastasis
Using exact test, we were able to identify differentially expressed miRNAs between BM + and BM- group. With cut off at logCPM > 4, p < 0.05, FDR ≤ 0.1, majority of miRNAs were excluded except hsa-miR-574-5p, hsa-miR-328-3p and hsa-miR-423-3p (Table 2). The CPM of three miRNAs in each group was shown in Fig. 3. Hsa-miR-574-5p was significantly down-regulated in BM+ (Fig. 3A). Studies on the function of hsa-miR-574-5p pointed out that hsa-miR-574-5p was a suppressor in Wnt/β-catenin signaling pathway and highly related to development and metastasis of different cancer including thyroid cancer[25], colorectal cancer[26, 27] and breast cancer[28]. Hsa-miR-328-3p and hsa-miR-423-3p were significantly up-regulated in BM + compared to BM- group (Fig. 3B, C). Both miRNAs were reported as activators in Wnt/β-catenin signaling pathway and promoted cancer cell invasion and metastasis in advanced non-small cell lung cancer[29] and colorectal cancer[30]. More importantly, the three DE miRNAs we identified, hsa-miR-574-5p, hsa-miR-328-3p and hsa-miR-423-3p, all belonged to the cluster B which is related with bone metastasis (Supplementary Table 3, Turquoise). These DE miRNAs might regulate cancer metastasis through Wnt/β-catenin signaling pathway and might be potential biomarker for bone metastasis in NSCLC.