Biological Validation of RNA Sequencing Data from Formalin-Fixed Paraffin-Embedded Primary Melanomas

Abstract

Purpose: Initiatives such as The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) have generated high-quality, multi-platform molecular data from thousands of frozen tumor samples. While these initiatives have provided invaluable insight into cancer biology, a tremendous potential resource remains largely untapped in formalin-fixed, paraffin-embedded (FFPE) samples that are more readily available, but which can present technical challenges due to crosslinking of fragile molecules such as RNA.

Materials and methods: We extracted RNA from FFPE primary melanomas and assessed two gene expression platforms -- genome-wide RNA sequencing (RNA-seq) and targeted NanoString -- for their ability to generate coherent biological signals. To do so, we generated an improved approach to quantifying gene expression pathways, in which we refine pathway scores through correlation-guided gene subsetting. We also make comparisons to the TCGA and other publicly available melanoma datasets.

Results: Comparison of the gene expression patterns to each other, to established biological modules, and to clinical and immunohistochemical data confirmed the fidelity of biological signals from both platforms using FFPE samples to known biology. Moreover, correlations with patient outcome data were consistent with previous frozen-tissue-based studies.

Conclusion: FFPE samples from previously difficult-to-access cancer types - such as small primary melanomas - represents a valuable and previously unexploited source of analyte for RNA-seq and NanoString platforms. This work provides an important step towards the use of such platforms to unlock novel molecular underpinnings and inform future biologically-driven clinical decisions.

Fig 1.

Iterative subsetting of gene sets to identify top cocorrelating genes. (A) Schematic of the algorithm. The Cancer Genome Atlas (TCGA) immune gene set is shown as a representative. The algorithm stops when the final iteration’s median equals that of the preceding iteration, having thus reached the fixed point. This final median becomes the pathway score. The orange box denotes the final “subsetted” genes that define the score. (B) Median iteration values for the TCGA immune gene set. Iteration 3 equals iteration 2 and so is not shown. (C) Median iteration values for the TCGA melanocyte gene set. (D) Final immune gene set, TCGA data. (E) Final melanocyte gene set, TCGA data.

Fig 2.

RNA sequencing (RNA-seq) and NanoString compare favorably on the same samples. (A) Schematic of the sample processing workflow. (B) Spearman correlation values for all 1,362 genes shared by the RNA-seq and NanoString datasets, as well as stratified by expression and dynamic range. Probes with low expression were assigned to the “low expression” category even if they had low dynamic range. Student t test was performed for nonoverlapping categories. (*)P < .001. (C) Topmost highly and poorly correlating genes between the platforms, presented as a heatmap of the RNA-seq data. (D) Median absolute deviation divided by the median (MAD/M), a measurement of dynamic range, for the RNA-seq data. Each dot is a gene corresponding to one in (C). Avg, average; Corr., correlation.

Fig 3.

Pathway-level correlations within and between platforms. (A) Heatmap of immune genes derived from Fig 1D for the RNA-seq and NanoString datasets, with the immune score on top. (B) Intraplatform correlation to the immune score for individual immune genes, plotted as RNA-seq versus NanoString. Several key immune genes are highlighted. (C-E) Correlation of the (C) immune, (D) cell cycle, and (E) skin scores between platforms. (F, G) Correlation of the EMT and pigment scores within (F) the TCGA primary melanomas and (G) the current FFPE RNA-seq dataset. EMT, epithelial-mesenchymal transition; RNA-seq, RNA sequencing; TCGA, The Cancer Genome Atlas.

Fig 4.

Correlation between pathway scores and histopathology. (A, B) Correlation of the RNA-sequencing cell cycle score to (A) pHH3 immunohistochemistry (IHC) and (B) mitotic rate. (C-F) Correlation of the RNA-sequencing CD8 or CD3 gene expression data to their respective IHC at either the (C, D) center or (E, F) periphery of the tumor.

Fig 5.

Correlation between pathway scores and patient data. (A, B) Correlation of the RNA-seq skin score to (A) Breslow thickness or (B) tumor surface area. (C) Anticorrelation of the RNA-seq pigment score to the chronic sun-damage score. (D) Gene-set enrichment analysis output showing mitotic gene sets as enriched in samples from patients in whom disease eventually recurred. (E) Correlation of RNA-seq or NanoString cell cycle scores to samples stratified by patient-level recurrence status. RNA-seq, RNA sequencing.

Fig. A1.

Example of preserved gene–gene correlation across datasets. Correlation of TOP2A and BUB1 in The Cancer Genome Atlas, GSE19234, GSE22155, and current formalin-fixed paraffin-embedded tissue datasets, all assessing patients’ melanoma samples.

Fig. A2.

Additional validation of the immune gene set. (A) Correlation of our The Cancer Genome Atlas (TCGA)-derived immune score and the official TCGA immune score. (B) Venn diagram illustrating the overlap in gene membership between our TCGA-derived immune gene set and the official TCGA immune gene set. (C) Top: Heatmap of the 38 formalin-fixed paraffin-embedded (FFPE) RNA-seq dataset for all 407 TCGA-derived immune gene-set genes. Bottom: The median absolute deviation divided by the median (MAD/M) for all 407 immune genes for each sample, a measure of correlative discordance within the immune signature. (D) Comparison of the immune MAD/M for the TCGA and FFPE datasets.

Fig. A3.

Heatmaps showing the final fixed-point gene sets for (A) epithelial-mesenchymal transition (EMT), (B) cell cycle, and (C) skin for RNA sequencing (RNA-seq) and NanoString. For each gene set, the RNA-seq and NanoString samples have the same order, from highest to lowest RNA-seq score. The melanocyte score genes are not available for NanoString because no melanocyte genes were present on the panel. (D) Formalin-fixed paraffin-embedded (FFPE) NanoString versus FFPE RNA-seq correlations, comparing the final gene members of the five pathways with all other genes shared between the two platforms.

Fig. A4.

Intragene-set correlations for the five pathways. Each gene-set member is correlated with its respective overall gene-set score, for primary melanomas only in the (A) The Cancer Genome Atlas (TCGA) and current formalin-fixed paraffin-embedded (FFPE) RNA-seq and (B) GSE7553 and GSE15605 datasets. The percentages at the top indicate that of genes in FFPE tissue with ρ > 0.5. (C) Hierarchical clustering of the top 697 most-variable pan-cancer genes for 21 cancer types plus our FFPE melanoma dataset. Shaded triangles highlight clusters. The yellow box indicates melanocyte-specific genes enriched across melanoma samples, and includes such genes as TYR, TYRP1, and MLANA. ACC, adrenocortical carcinoma; BLCA, bladder carcinoma; BRCA, breast carcinoma; COADREAD, colorectal adenocarcinoma; EMT, epithelial-mesenchymal transition; FFP, this study; GBM: glioblastoma multiforme; HNSC, head and neck squamous carcinoma; LIHC, hepatocellular carcinoma; KIRC, renal clear-cell carcinoma; KIRP, renal papillary carcinoma; LAML, acute myeloid leukemia; LGG, low-grade glioma; LUAD, lung adenocarcinoma; LUSC, lung squamous carcinoma; Mets, metastases; MESO, mesothelioma; PAAD, pancreatic adenocarcinoma; PRAD, prostate adenocarcinoma; OV, ovarian; SKCM met, metastatic cutaneous melanoma; SKCM pri, primary cutaneous melanoma; THCA, thyroid carcinoma; UCEC, uterine carcinoma; UVM, uveal melanoma.

Fig. A5.

(A) Hierarchical clustering of the top 5% most variable genes. Enriched gene classes are shown for three clusters on the basis of Metacore analyses. The final cluster did not have a clear enrichment. We note that the MITF/melanocyte pathway was not recovered, because of the smaller dataset (n = 38 samples here v n = 333 in TCGA) and the lower variability in melanocyte gene expression, resulting in those genes not being among the top 5%. (B) Average reads per kilobase million (RPKMs) for all genes within the five gene sets used in this study and for three gene sets not expected to be highly expressed in melanoma. Percentages on top indicate the percentage of genes within each gene set above the 1.0 average RPKM threshold (dashed line). The three control gene sets were derived from the GTEx normal tissue database. Error bars refer to median plus interquartile range. The y-axis scale is log₁₀.

Fig. A6.

Correlation between pathway scores and histopathology for NanoString data. (A, B) Correlation of the NanoString cell cycle score with (A) pHH3 immunohistochemistry (IHC) and (B) mitotic rate. (C-F) Correlation of the NanoString CD8 or CD3 gene expression data with their respective IHC at either the (C, D) center or (E, F) periphery of the tumor.