Integrative radiogenomic analysis for genomic signatures in glioblastomas presenting leptomeningeal dissemination

Abstract

Despite therapeutic advances, the prognosis for glioblastoma (GBM) remains poor. In particular, leptomeningeal dissemination (LMD) has a dismal prognosis. The aim of this study was to identify tumor molecular phenotype, which has a great propensity to develop LMD. Between May 2004 and December 2012, a total of 145 GBM tumor samples were obtained from data registry. A total of 20 of the 145 patients with GBM were found to develop LMD. A specialized radiologist confirmed the diagnosis of LMD on magnetic resonance imaging. To clarify the genomic signatures in GBM with LMD, we performed integrative analysis of whole transcriptome sequencing and copy number alteration in the radiological features indicating LMD phenotypes in GBM. Eleven newly diagnosed patients with GBM with LMD had worse prognosis than those without LMD (median 5.55 vs. 12.94 months, P < 0.0001). Integrating analysis using gene expression based on the change of copy number revealed that SPOCK1, EHD2, SLC2A3, and ANXA11 were highly expressed with the gain of copy number, compared with the gene expression in the non-LMD group. In addition, it was demonstrated that NME2, TMEM100, and SIVA1 were downregulated with the loss of copy number. We also found that mesenchymal subtype accounted for 50% in LMD group, whereas mesenchymal subtype consisted of 29% in non-LMD group, even though there was no statistical significance (P = 0.06). Through this radiogenomic analysis, we suggested the possibility of finding candidate genes associated with LMD and highlighted the significance of integrating approach to clarify the molecular characteristics in LMD.

1 Introduction

Despite advances in surgical techniques, radiotherapy, and chemotherapy, glioblastoma (GBM) has still a poor prognosis. In particular, leptomeningeal dissemination (LMD) in the progression of GBM carries a worse prognosis. This condition has been considered a rare and serious condition, in which survival ranges from 12 to 20 weeks.^[1–3] To date, only a few studies have reported an accurate incidence of LMD in GBM.^[1,4] The overall incidence of LMD in GBM remains low. There have been potential risk components of LMD, specifically related to geographical location of the tumor or the proximity to the ventricular system. However, LMD can be sometimes found primarily in initial presentation as well as in the progressive stage of GBM regardless of communication with ventricle or subarachnoid cistern. Recent progress in the treatment of GBM has contributed to the prolongation of survival in GBM. This prolonged survival has led to the speculation that LMD has been reported with increasing frequency in recent years, occurring in up to 20% to 25% of newly diagnosed or recurrent GBM.^[2,5–7] In addition, recent development of new molecular-targeting drugs is attributable to convert the tumor biology to a more invasive phenotype.^[7–9] Mandel et al reported that increase in tumor invasiveness and predilection to disseminate into the cerebrospinal fluid (CSF), causing LMD, might be primarily related to molecular characteristics of the tumor.^[7] To date, molecular signature associated with the development of LMD has not been identified. To address the hypothesis that GBM with LMD has worse prognosis and reflects underlying intertumoral molecular profiles, we performed integrative analysis of whole transcriptome sequencing (RNA-Seq) gene expression patterns, copy number alteration (CNA), and whole exome sequencing (WES) in the radiological features indicating LMD phenotypes in GBM.

2 Materials and methods

2.1 Patient sample preparation

Between May 2004 and December 2012, a total of 195 GBM tumor samples (145 patients) with available clinical and pathology reports were obtained from data registry in Samsung Medical Center (Seoul, Korea) because this group had the most extensive and complete follow-up data. All tissue samples had a pathologically confirmed diagnosis of GBM according to WHO criteria and collected with written informed consent under a protocol approved by the Institutional Review Board of the Samsung Medical Center (2010-04-004). All samples were from adult persons. Exclusion criteria included patients without histological confirmation of grade IV GBM and without agreement of this study. We identified 20 patients with GBM who developed LMD from our institutional database. The diagnosis of LMD was determined by radiographic findings on magnetic resonance (MR) imaging. Diagnostic lumbar puncture for CSF cytology was not performed, because the definition of LMD was not limited to the spinal arachnoid space, but to the intracranial dissemination. Data collected included demographic characteristics (age and sex), tumor characteristics (newly diagnosed or secondary GBM), survival from GBM diagnosis, time from GBM diagnosis to LMD diagnosis, and survival from the time of LMD diagnosis.^[7] There were 120 newly diagnosed GBM and 25 recurrent setting of GBM. Patients consisted of 83 males and 62 females and the median age of the patients was 53 years (range, 16–80 years) (Table 1).

Table 1: Characteristics of glioblastoma samples by status of LMD.

2.2 Definition of leptomeningeal dissemination on MR finding

The LMD was defined as positive findings in imaging evaluation of brain. All MRI exams were performed on a 1.5- or 3.0-T Signa Echospeed scanner (GE Medical Systems, Milwaukee, USA). To diagnose LMD, we set the radiographic criteria that were contrast enhancement of the leptomeninges around the outlines of the gyri and sulci or multiple nodular deposits in the subarachnoid space on T1 and/or T2 fluid attenuated inversion recovery gadolinium-enhanced MR images (Fig. 1). Imaging characteristics were established through a consensus of specialized neuroradiologists. All images were evaluated by consensus in a blinded fashion by 2 board-certified radiologists (S.T. Kim and J.H. Cha). Both readers were blinded to the genomic signatures and other clinical details at the time of image interpretation.

Figure 1: Schematic outline of this study.

2.3 Genomic DNA and RNA isolation

DNA was isolated using the DNeasy kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol. Total RNA was extracted from brain tumor tissue using the RNeasy kit (Qiagen, Hilden, Germany).

2.4 Survival analysis

Imaging and clinical follow-up were available for all patients. Patients were binary classified as presenting LMD or not on MR imaging as described earlier. The mainstay of survival analysis was the Cox proportional hazards model using R 3.0.1 (Vienna, Austria; http://www.R-project.org/) and P value <0.05 was deemed statistically significant. In order to assess the effects of a given predictor, we generally used the log-rank test.

2.5 Whole exome sequencing and copy number alteration

Capturing exonic DNA fragments was performed using either the Illumina TruSeq Exome-capture kit (for Case S780) or the Agilent SureSelect kit (for 100 other cases). Captured exonic DNA fragments were sequenced on the Illumina HiSeq2000 to generate 2 × 101 bp paired-end reads. BAM file was created by aligning paired-end WES reads using Burrows-Wheeler Aligner (version 0.6.2) to the human reference genome (hg 19). Sequences were sorted using SAMtools (version 0.1.18), and the duplicate reads were removed using Picard (version 1.73, http://picard.sourceforge.net). The predicted somatic mutations were annotated using Variant Effect Predictor (version 73).^[10] When the matched normal WES data were not available, we used the mpileup command in SAMtools^[11] with the preprocessed BAM files for tumor samples. The mpileup results were processed further to extract bases with a base quality of at least 15. To assemble a list of high-value mutations for initial validation, we selected all mutations meeting the following criteria: mutant allele read count ≥2, mutant allele fraction ≥0.01, and chromosomal coordinate of a mutation recorded as “confirmed somatic mutation” in COSMIC database (version 63).^[12] Insertions/deletions (indels) were not called for tumors without the matched normal WES data. We used the patient-matched normal WES data to estimate the fold change in copy number in tumors. When the cases are not available to the patient-matched normal data, we made a “pseudo-normal” profile according to averaging a pool of 20 randomly chosen normal WES data, generated using the same sequencing platform and analysis pipeline as the tumor data. To do downstream analyses, including segmentation and calculation of the copy number for each gene, we used the ngCGH python package (version 0.4.4, http://github.com/seandavi/ngCGH), which is a tool for producing aCGH-like of next-generation sequencing data.

2.6 Whole transcriptome sequencing (RNA-Seq)

For each samples, we used the illumine TruSeq RNA Sample preparation kit to prepare RNA-Seq libraries. For analysis of mRNA level, we obtained the trimmed reads, which include 30 nucleotides from the 5′ end of each read. The trimmed reads were aligned to the human reference genome (hg19) using GSNAP^[13] (version 2012-12-20). We used R package DEGseq^[14] to calculate gene expression as reads per kilobase per million mapped reads with RefSeq gene annotation (“refFlat” table downloaded from UCSC genome browser, last accessed on 2012/08/06).

2.7 Integration of copy number alteration and genome-wide expression analyses

RNA deep sequencing datasets were available in 144 of 145 patients and CNA datasets in 100 of 145 patients. Data from CNA and genome-wide expression profiles were analyzed individually. To identify the impact of CNA on gene expression, we utilize the following statistical approach. WES data were log-transformed. For each gene, the WES data were represented by a vector that was labeled 1 for amplification (log value >0) and 0 for no amplification. To identify the significant genes that exhibited CNA and gene expression changes, we modified an R package (CNAmet) that integrates high-throughput copy number data and expression data.^[15] CNAmet algorithm followed next step. First, the signal-to-noise ratio statistic was used to link gene expression values to CNAs^[16]:

where m₀ and σ₀, and m₁ and σ₁ represent the sample means and sample standard deviations for the expression level for nonamplified and amplified samples, respectively. Second, the weight values were combined to a score that indicated genes whose expression alterations were attributed to changes in copy number levels. If the difference of means of the groups is bigger and standard deviations within the groups are small, signal-to-noise statistics results in a large weight. Finally, we calculated corrected α values with a permutation test. We used permutation test in order to attain statistical significance. We conducted 10,000 permutations and got an α value. A low α value represents a close association between gene expression and gene amplification.

3 Results

3.1 GBM with LMD has worse clinical outcome than GBM without LMD

A total of 20 patients out of the 145 patients were identified to have developed LMD, yielding a rate of 13.8% of patients for this group. Eleven patients showed LMD in newly diagnosed GBM and 9 other patients showed LMD in the recurrent GBM (Table 1). Newly diagnosed tumor location was predominantly supratentorial (17 of 20 cases, 85%) and the remainder (3 cases) were infratentorial lesions. We analyzed the possible clinical factors to assess whether or not they affected LMD. The median follow-up period was 356 days. During the follow-up period, 5 of 20 patients with LMD survived, and 52 of 125 patients without LMD survived. In a total of 145 patients, median survival from the diagnosis as GBM was 12.3 months. In 20 patients with LMD, the interval between diagnosis of GBM and the occurrence of LMD ranged from 0 to 1442 days (mean, 214.85 days). Excluding the recurrent GBM in LMD group, from the diagnosis of GBM, overall survival analysis revealed that newly diagnosed patients with GBM with LMD (11 patients) had worse prognosis than those without LMD (median 5.55 vs. 12.94 months, P < 0.0001) (Fig. 2). Other clinical parameters such as age, performance status, and extent of resection were not significantly different between the groups (Table 1) (S1 Fig., https://links.lww.com/MD/B86).

Figure 2: Kaplan–Meier analysis of overall survival between newly diagnosed patients with glioblastoma with and without LMD. P values are according to the log-rank test. LMD = leptomeningeal dissemination.

3.2 Integrative approach for identification of differentially expressed genes

Compared to single type–based approach (gene expression or CNA), this proposed integrative approach (CNAmet) yielded a considerably smaller number and fewer portions of such synergistically impacted genes. As a result, it provided prevention from overestimation of candidate genes and a large number of genes fell into the class impacted definitely by only 1 genomic type under our integrative framework. As a result of integrating analysis using gene expression based on the change of copy number, we found that SPOCK1, SLC6A6, ODZ3, EHD2, EFNA5, SLC2A3, KIF3C, ANXA11, STK10, and C10orf90 were highly expressed with the gain of copy number, compared with the gene expression in the non-LMD group. In addition, it was demonstrated that GUK1, FAM195B, NUBP2, ROMO1, NME2, ANAPC11, TMEM100, RPL39L, NUDFB7, and SIVA1 were downregulated with the loss of copy number (Fig. 3A). The gene list can be found in Table 2. We performed GO and KEGG pathway analysis using Expression2Kinase software (http://www.maayanlab.net/X2K).^[17] By importing the complete list of upregulated genes in LMD or non-LMD into the KEGG pathway and GO categories, we found that highly expressed genes in the LMD group were related to membrane part such as “anchored to plasma membrane (GO:0046658)” and “HSA04360_AXON_GUIDANCE” (Supplemental Table 1A, https://links.lww.com/MD/B86). Also we found that highly expressed genes in the non-LMD group were related to ubiquitination such as “HSA04120_UBIQUITIN_MEDIATED_PROTEOLYSIS” and “NADH dehydrogenase (ubiquinone) activity (GO:0008137)” (Supplemental Table 1B, https://links.lww.com/MD/B86).

Figure 3: Gene expression and CNA across LMD status. (A) Heat map demonstrating supervised hierarchical clustering of gene expression based on the change of CNA between patients with glioblastoma with and without LMD. Gene list defined using fold change cutoff (≥1.2 or ≤0.7) and P-value cutoff (<0.05). Top 10 genes and bottom 10 genes are represented with fold change, z-score. Red to blue color scale indicates the range from the highest positive to the highest negative correlation. (B) The distribution of CNA delineates amplification/deletion status in LMD and non-LMD. Chromosomes 1 to 22 are represented from left to right. Y-axis shows the relative number of patients according to status. Amplification means CNA ≥2.8 copy and deletion means CNA ≤1.4 copy. CNA = copy number alteration, LMD = leptomeningeal dissemination.

Table 2: Description of 104 differentially expressed genes list integrating the change of copy number variation.

Table 2: (Continued) Description of 104 differentially expressed genes list integrating the change of copy number variation.

3.3 Molecular subtypes in GBM with leptomeningeal dissemination

Recently, GBM populations have been clustered into 5 molecular subtypes (proneural, neural, classical, mesenchymal, and unknown) based on gene expression profiles.^[18,19] We found that mesenchymal subtype accounted for 50% in LMD group, whereas mesenchymal subtype consisted of 29% in non-LMD group (Table 1). There was an increasing trend toward significance in the mesenchymal subtypes in LMD group, although there was no statistical significance (P = 0.06). We mapped the above-mentioned reference genes onto the 5 molecular subgroups defined by the TCGA network in non-LMD group (S2 Fig., https://links.lww.com/MD/B86).

3.4 Somatic mutation and copy number alteration in LMD group

Approach for GBM-related somatic mutation and CNA was to identify the characteristics of LMD group. There was no significant difference between patients with GBM with LMD and those without LMD in CNA (Fig. 3B). In detail, we found PTEN deletion (50%), EGFR amplification (35%), CDKN2A/B deletion (40%), CDK6 amplification (30%), MET amplification (30%), and RB1 deletion (15%) in 20 LMD samples (Fig. 4). Sixteen of 20 samples had available indel information, whereas 4 LMD samples were not available due to the matched normal pairs. PTEN, TP53, ATRX, and NF1 indel were very rare in LMD group. No IDH1 mutation (R132H) was found in the LMD group, compared with the incidence of 6 of 125 in the non-LMD group. We also found that only a few samples harbored EGFR mutation (A289V) (1 LMD and 4 non-LMD samples). Accordingly, there were no significant somatic mutations in the LMD group.

Figure 4: Key gene driver alterations in patients with glioblastoma with LMD. The 16 cases have WES data available for both the tumor and the paired normal tissues. WES has enabled detection of point mutations and indels (MAF ≥0.2 and read depth ≥20) and copy number changes (≥0.5 or ≤−0.5, log 2 scale). CN = copy number, CN (corr) = corrected copy number, Expr = expression (mRNA), indel = insertion/deletion, LMD = leptomeningeal dissemination, MAF = mutant allele fraction, RPKM = read per kilobase per million, SNV = single nucleotide variant, WES = whole exome sequencing.

4 Discussion

Understanding the mechanism of leptomeningeal spreading and identification of genomic data leading to LMD is very critical to prolong the survival in patients with GBM. Distant metastasis of intracranial GBM such as LMD or spinal cord metastasis has been described with increasing frequency in recent years because of higher survival rates and prolonged survival times.^[4] To date, the risk factors associated with LMD are not clearly understood.^[2,20] Although some studies have demonstrated that communication with ventricular system is a critical factor for LMD in patients with malignant gliomas,^[20–22] other studies reported that the proximity of tumor margin to ventricular system was not a risk factor for CSF dissemination.^[4,23] This study revealed no significant association between ventricular opening and LMD, because adjacent cortical spread of LMD as well as distant leptomeningeal seeding around the brainstem or ventricular system was included in this study. The outcome of LMD in this study was compatible with the poor prognosis of patients with GBM with LMD as previously described.^[2] Similarly, GBM with LMD at the initial presentation had a worse prognosis to newly diagnosed GBM without LMD in this study (Fig. 2). It still remains controversial whether the prognosis in patients with GBM with primarily presented LMD is poor or not.^[24] In an attempt to clarify the risk factors for LMD, we investigated the molecular features in GBM with LMD in this study.

In the recent literature, there has been a little information about molecular characteristics in the LMD. Korshunov et al showed that gains of the 1p36 locus were closely related to symptomatic LMD of supratentorial GBM.^[24] Another study suggested that mutations of the PTEN gene and a high MIB-1 labeling index were correlated with LMD in GBM.^[25] In this study, we found that some genes were highly expressed in the LMD group based on the CNA change. As representative genes among them, SPOCK1 can activate PI3K/Akt signaling to block apoptosis and promote proliferation and metastasis in vitro and in vivo.^[26] It is also known as a novel metastasis-related biomarker in lung cancer^[27] and promotes hepatocellular carcinoma regulated by CHD1L.^[28] EHD2, EH-domain containing 2, has a role of migration and invasion in the breast cancer cells.^[29] As another candidate gene associated with LMD, ANXA11 and SLC2A3 are known to play a role of cellular invasion and metastasis in laryngeal cancer and hepatic carcinoma.^[30,31] In particular, the mutation and deregulation of ANXA11 was reported to enhance cancer metastasis, invasion, and drug resistance through the platelet-derived growth factor receptor pathway and the mitogen-activated protein kinase/P53 pathway.^[32] Among the downregulated genes in LMD group, NME2, TMEM100, and SIVA1 can be candidate genes promoting non-LMD. NME2 is known to reduce proliferation, migration, and invasion of gastric cancer cells to limit metastasis.^[33]

We acknowledge that this study has some limitations such as relatively small sample size to draw a definite conclusion. In addition, the diagnosis of LMD in this study did not include cytological confirmation. In this study, the diagnosis of LMD was entirely based on the recent advances of radiological MR techniques, because CSF study has low sensitivity of the diagnosis of intracranial LMD, despite a higher specificity. Despite these limitations, we presented a trend of transcriptomic characterization between LMD and non-LMD cohorts based on genomic profiles for suggestive candidate genes associated with LMD prognosis.

Through this radiogenomic analysis integrating gene expression, CNA, and mutation analysis, we suggested the possibility of finding candidate genes associated with LMD and highlighted the significance of integrating approach to clarify the molecular characteristics in LMD. To prove several candidate genes associated with LMD, we believe that further molecular experimental studies should be investigated in the future.