Abstract
Melanoma is the most serious type of skin cancer that frequently
spreads to other organs of the human body. Especially melanoma
metastases to the brain (intracranial metastases) are hard to treat and
a major cause of death of melanoma patients. Little is known about
molecular alterations and altered mechanisms that distinguish intra-
from extracranial melanoma metastases. So far, almost all existing
studies compared intracranial metastases from one set of patients to
extracranial metastases of an another set of melanoma patients. This
neglects the important facts that each melanoma is highly individual
and that intra- and extracranial melanoma metastases from the same
patient are more similar to each other than to melanoma metastases from
other patients in the same organ. To overcome this, we compared the
gene expression profiles of 16 intracranial metastases to their
corresponding 21 patient-matched extracranial metastases in a
personalized way using a three-state Hidden Markov Model (HMM) to
identify altered genes for each individual metastasis pair. This
enabled three major findings by considering the predicted gene
expression alterations across all patients: (i) most frequently altered
pathways include cytokine-receptor interaction, calcium signaling,
ECM-receptor interaction, cAMP signaling, Jak-STAT and PI3K/Akt
signaling, (ii) immune-relevant signaling pathway genes were
downregulated in intracranial metastases, and (iii) intracranial
metastases were associated with a brain-like phenotype gene expression
program. Further, the integration of all differentially expressed genes
across the patient-matched melanoma metastasis pairs led to a set of
103 genes that were consistently down- or up-regulated in at least 11
of the 16 of the patients. This set of genes contained many genes
involved in the regulation of immune responses, cell growth, cellular
signaling and transport processes. An analysis of these genes in the
TCGA melanoma cohort showed that the expression behavior of 11 genes
was significantly associated with survival. Moreover, a comparison of
the 103 genes to three closely related melanoma metastasis studies
revealed a core set of eight genes that were consistently down- or
upregulated in intra- compared to extracranial metastases in at least
two of the three related studies (down: CILP, DPT, FGF7, LAMP3, MEOX2,
TMEM119; up: GLDN, PMP2) including FGF7 that was also significantly
associated with survival. Our findings contribute to a better
characterization of genes and pathways that distinguish intra- from
extracranial melanoma metastasis and provide important hints for future
experimental studies to identify potential targets for new therapeutic
approaches.
Keywords: Melanoma metastases, Patient-matched intra- and extracranial
melanoma metastasis pairs, Personalized transcriptome analysis
Background
Patients with melanoma brain metastasis (intracranial melanoma
metastasis) still have a very unfavorable prognosis. The mean overall
survival of patients with untreated intracranial melanoma metastases is
as little as 4 months [[35]1], improving only to a median survival of
22.7 months with surgery followed by immunotherapy [[36]2]. There are
several emerging treatment options like targeting MAPK signaling using
BRAF or MEK inhibitors or immunotherapy using anti-CTLA4/anti-PD-1
antibodies that improve patient survival. These treatment options are
effective for patients with extracranial metastases [[37]3]. For
patients with intracranial metastases, the treatment response duration
of BRAF/MEK inhibitors is only limited to a few months and the efficacy
of immune checkpoint inhibitors is substantially reduced in symptomatic
patients [[38]1, [39]4–[40]8]. Results of important clinical
immunotherapy trials for melanoma patients with intracranial metastases
are summarized in [[41]9]. Clinical outcomes of melanoma patients with
intracranial metastases treated with stereotactic radiosurgery and
various systemic therapies have been analyzed in [[42]10]. The overall
12 month survival rates for the combined anti-PD-1-CTLA4 therapy was
68%, 62% for BRAF/MEK inhibitor treatment, 59% for anti-PD-1 therapy,
and 45% for anti-CTLA4 therapy compared to only 21% for BRAF inhibitor
treatment and 15% for conventional chemotherapy. The combined treatment
of patients with anti-PD-1-CTLA4 therapy showed the best median overall
survival of about 21 months. The updated five year data from patients
with intracranial melanoma metastases enrolled on the ABC trial
investigating nivolumab plus ipilimumab or nivolumab alone confirmed
the high anti-tumor activity of nivolumab plus ipilimumab in patients
with asymptomatic intracranial metastases (5-year intracranial
progression-free survival 52%, 5-year overall survival 55%) [[43]11].
The NIBIT-M2 trial showed persistent therapeutic efficacy of ipilimumab
plus nivolumab with a seven year overall survival rate of 42.8% in
asymptomatic patients [[44]12]. Further, in a recently published
retrospective study with 376 patients with intracranial melanoma
metastases treated with nivolumab plus ipilimumab, long-term survival
was seen in treatment-naive, asymptomatic, steroid-free patients as
well as in those patients that received stereotactic radiosurgery in
combination with nivolumab plus ipilimumab [[45]13].
Still, the therapy resistance of intracranial metastases is the leading
cause of death of melanoma patients [[46]14, [47]15]. Unfortunately,
intracranial metastases are also very common affecting about 50% of all
metastatic melanoma patients [[48]1]. Therefore, a detailed molecular
characterization of differences between intra- and extracranial
metastases is needed to better characterize molecular alterations and
mechanisms that distinguish both types of metastases to provide a basis
for future developments of novel therapeutic strategies.
Several studies were done over the last years to identify differences
between intra- and extracranial melanoma metastases. Different genetic
alterations in genes like BRAF, NRAS and CDKN2A were suggested to
contribute to a unique molecular profile of intracranial metastases
[[49]16]. Up-regulation of PI3K/Akt signaling in intracranial
metastases has been reported as a key mechanism for uncontrolled
proliferation of intracranial metastases [[50]17–[51]20].
Epigenetically regulated genes (e.g. CSSP1, GRB10, NMB, PDXK, PRKCZ,
RASL11B, STK10 and WDR24) with altered promoter methylation and
corresponding differential expression in intra- compared to
extracranial metastases were recently identified [[52]21]. Single-cell
sequencing of metastases revealed a neuronal-like cell meta-program of
intracranial melanoma metastases [[53]22]. Another single-cell
sequencing study delineated brain metastases programs into a
proliferative and an inflammatory archetype [[54]23]. All these studies
compared intracranial metastases from several patients to extracranial
metastases from other patients. However, melanoma metastases in
different organs from a patient are more similar to each other than to
metastases from other patients in the same organs. This
patient-specific heterogeneity of molecular data from melanoma
metastases has been observed in other studies before [[55]18, [56]24,
[57]25]. Therefore, such differences between individual patients should
be included in the data analysis to further improve the identification
of molecular differences between intra- and extracranial melanoma
metastases.
To account for the inter-patient heterogeneity, some studies already
started to investigate patient-matched metastasis pairs at different
molecular layers. Chen et al. [[58]18] analyzed gene mutations, DNA
copy number alterations and gene expression profiles identifying the
PI3K/Akt signaling pathway as a potential therapeutic target. Fischer
et al. [[59]26] performed RNA sequencing and multiple immune-relevant
sequencing analyses and identified a significant immunosuppression and
enrichment of oxidative phosphorylation in intracranial melanoma
metastases. We have recently performed a personalized analysis of
genome-wide DNA-methylation profiles and revealed a global decrease of
DNA-methylation intracranially [[60]25]. These DNA-methylation changes
between patient-matched intra- and extracranial melanoma metastases
affected many genes involved in cellular signaling, growth, adhesion
and apoptosis and further supported the presence of a neuronal
phenotype. Further, we were also able to predict potential downstream
targets of genes with altered promoter methylation, which allowed to
group heterogeneous patient-matched melanoma metastasis pairs into
three homogeneous subgroups utilizing a network-based approach
[[61]27]. In addition, mutations in driver genes (most frequently
ARID1A, ARID2 and BRAF), which distinguished intra- from extracranial
melanoma metastases, were identified by targeted next-generation
sequencing [[62]28].
Here, we perform a personalized analysis of patient-matched gene
expression profiles of intra- and extracranial melanoma metastasis
pairs to account for the common developmental origin of patient-matched
metastases. To realize this, each patient-matched metastasis pair was
analyzed by a Hidden-Markov Model (HMM) approach to identify
differentially expressed genes for each pair. This was done by
transferring the recently used HMM-approach for the personalized
analysis of DNA-methylation profiles by [[63]25] to the analysis of
gene expression profiles. The predicted differentially expressed genes
were further used to identify frequently affected cellular pathways and
to derive a set of genes that was altered in the same manner in the
majority of patients. In-depth literature analysis in combination with
comparisons to independent related studies were performed to provide
further hints which genes potentially play an important role to
establish molecular differences between intra- and extracranial
melanoma metastases. Additionally, significant expression associations
of several of these genes with decreased survival of melanoma patients
from a large public patient cohort indicate the relevance of our
findings.
Methods
Gene expression data of melanoma metastases
The considered gene expression data set comprises 37 melanoma
metastasis samples from 16 patients. Each of these patients developed
an intra- and an extracranial melanoma metastasis. The extracranial
metastases included lung, lymph node, liver, skin, small intestine and
soft tissue metastases. A board-certified pathologist had marked
metastasis regions with a high percentage of tumor cells of at least
90% and minimal percentage of normal tissue, necrosis or hemorrhage. In
addition, histologically different regions in the extracranial
metastases of four patients (P04, P08, P18, P42) were included. The
gene expression profiles of these distinct histological regions were
included as separate patient-specific samples, because such distinct
regions in a metastasis may represent different subclones within a
metastasis. An overview of the data set and the patient-specific sample
composition is given in Table [64]1. Additional information about the
patients are provided in Additional file [65]18: Table S13. Gene
expression profiles of all metastases samples were measured by
RNA-sequencing (RNA-seq) and processed as described in [[66]21].
Details to the RNA sequencing protocol and the preprocessing of the
reads are provided in Additional file [67]22: Text S3. The considered
RNA-seq data of seven patients were taken from [[68]21] considering all
patients for which transcriptomes of patient-matched intra-
vs. extracranial melanoma metastasis pairs were available. The RNA-seq
data of the nine other patients were additionally measured for this
study to further increase the patient cohort. The expression read
counts of all metastases samples were normalized with the voom function
of the R-package limma performing a cyclic loess normalization
[[69]29]. To avoid the inclusion of very weakly expressed genes, only
protein-coding genes of the human genome annotation (hg19) with more
than one count per million (CPM > 1) in more than 25% of the metastases
samples were considered. This resulted in a gene expression data set
that comprised the absolute expression levels (log[2]-CPMs) of 14,946
genes in their chromosomal order across the 37 melanoma metastases
samples (Additional file [70]6: Table S1).
Table 1.
Melanoma metastasis patient and sample overview
Patient Intracranial metastasis Extracranial metastasis Data origin
P03 1 brain sample 1 lung sample Westphal et al. [[71]21]
P04 1 brain sample 2 skin samples Westphal et al. [[72]21]
P08 1 brain sample 3 soft tissue samples Westphal et al. [[73]21]
P13 1 brain sample 1 lymph node sample Newly sequenced
P16 1 brain sample 1 lung sample Westphal et al. [[74]21]
P18 1 brain sample 2 lung samples Westphal et al. [[75]21]
P39 1 brain sample 1 lung sample Westphal et al. [[76]21]
P42 1 brain sample 2 lymph node samples Westphal et al. [[77]21]
P74 1 brain sample 1 lymph node sample Newly sequenced
P77 1 brain sample 1 lymph node sample Newly sequenced
P78 1 brain sample 1 small intestine sample Newly sequenced
P101 1 brain sample 1 liver sample Newly sequenced
P106 1 brain sample 1 lymph node sample Newly sequenced
P107 1 brain sample 1 lung sample Newly sequenced
P108 1 brain sample 1 lymph node sample Newly sequenced
P111 1 brain sample 1 lymph node sample Newly sequenced
[78]Open in a new tab
Each melanoma patient developed an intra- and an extracranial
metastasis in the course of its disease. Extracranial metastases
developed either in lung, lymph node, skin, liver, small intestine or
soft tissue. Multiple samples of the same metastasis were taken if the
metastasis contained histologically different regions. The
corresponding metastasis transcriptomes of 7 of 16 patients were taken
from our previous study by Westphal et al. [[79]21] and the metastasis
transcriptomes of 9 of 16 patients were newly sequenced. Additional
information about the age of the patients at the resection of the brain
metastasis, sex, therapies, and mutational states of BRAF and NRAS are
provided in Additional file [80]18: Table S13
Gene expression data of normal tissues
RNA-sequencing data of 28 normal tissue samples were available to
analyze the purity of the melanoma metastases samples. These normal
tissue samples comprised six intracranial and 22 extracranial normal
tissues (8 lymph nodes, 7 lung and 7 skin/soft tissues) from 27
patients. In total, 11 of these samples were taken from [[81]21] and 17
samples were newly sequenced. These normal tissues were jointly
normalized with the melanoma metastasis samples using the same
procedure as described before. The normal tissue samples of three
patients (P03, P04, P16) matched with our melanoma metastases cohort in
Table [82]1. All considered normal tissue samples had a tumor cell
content of 0%. The normalized gene expression data of the normal
tissues are provided in Additional file [83]7: Table S2. An overview of
all normal tissue samples is given in Additional file [84]8: Table S3.
Hierarchical clustering of melanoma metastases and joint clustering with
normal tissue
The gene expression profiles of all melanoma metastases samples
(Additional file [85]6: Table S1) were clustered hierarchically to
identify similarities and differences between the metastases. The
clustering was performed using the R package pvclust [[86]30] with
10,000 bootstrapping repetitions to obtain stability estimates for all
clusters. The Manhattan distance was used as distance measure between
individual samples and Ward’s minimum variance linkage method (ward.D2)
[[87]31] was used as cluster algorithm. The stability of a cluster was
quantified by the approximate unbiased p-value (AU value). The AU value
ranges from 0 to 100, where a larger value represents greater stability
and a value of 100 means that a cluster is perfectly stable.
In addition, to further test and ensure the high tumor content in our
metastasis data, all transcriptomes of metastases and normal samples
were also hierarchically clustered together (Additional file [88]1:
Figure S1). This clustering and its stability analysis were performed
using the same methods as described above.
HMM-based analysis of gene expression profiles of patient-specific metastasis
pairs
Transcriptomes of metastases samples from the same patient were more
similar to each other than to metastases samples of other patients.
Therefore, a personalized analysis of each patient-matched intra- and
extracranial melanoma metastasis pair was done to predict individual
gene expression changes between the intra- and the extracranial
metastasis of each pair. First, pair-specific log[2]-ratio gene
expression profiles were computed for all 21 possible patient-matched
pairs (Additional file [89]9: Table S4). Such a relative expression
profile was obtained for each pair by subtracting for each gene its
absolute expression (log[2]-expression value from Additional file
[90]6: Table S1) in the extracranial metastasis from the absolute
expression of this gene in the intracranial metastasis followed by
sorting of the genes by their chromosomal order. The obtained
log[2]-ratios quantify the expression changes of genes: (i) a
log[2]-ratio clearly less than zero indicates reduced expression, (ii)
a log[2]-ratio about zero indicates unchanged expression, and (iii) a
log[2]-ratio clearly greater than zero indicates increased expression
of a gene in the intra- compared to the extracranial metastasis of a
patient-matched metastasis pair (Additional file [91]2: Figure S2).
Next, an autocorrelation analysis of the log[2]-ratio gene expression
profiles was performed in a chromosome-specific manner. Each
autocorrelation profile of each chromosome was weighted according to
the number of genes that were measured on that chromosome and then
averaged across all chromosomes for the corresponding lags. This
weighted autocorrelation analysis revealed that genes in close
chromosomal proximity tend to correlate more strongly in their
expression with each other than expected by chance (Fig. [92]2A).
Fig. 2.
[93]Fig. 2
[94]Open in a new tab
Autocorrelations of gene expression levels in close chromosomal
proximity and illustration of the utilized Hidden Markov Model for
decoding of gene expression states. A, Autocorrelations in the
chromosomal order of genes (red) are significantly greater than the
autocorrelations of 1000 randomly permuted gene expression profiles
(black). The autocorrelation was calculated chromosome-wise and
weighted according to the number of genes on the chromosome. Decreasing
autocorrelations of log[2]-expression-ratios between intra- and
extracranial metastases of genes in chromosomal order with increasing
lag between the genes show that chromosomal distant genes have less
similar expression than genes that are closer to each other on a
chromosome. The yellow ribbon shows the standard deviations of the
observed autocorrelations. B, Illustration of the three-state Hidden
Markov Model (HMM) with state-specific Gaussian emission densities that
was used to perform the personalized analysis of the patient-matched
metastasis pairs. Genes with unchanged expression in the intra-
compared to the corresponding extracranial metastasis are assigned to
the state ’
[MATH: = :MATH]
’, genes with decreased expression in the intra- compared to the
extracranial metastasis are assigned to the state ’−’, and genes with
increased expression in the intra- compared to the extracranial
metastasis are assigned to the state ’
[MATH: + :MATH]
’. The arrows that connect the states represent possible state
transitions for directly neighboring genes on a chromosome and the
corresponding values represent the learned state transition
probabilities
This observation motivated the usage of a Hidden Markov Model (HMM) for
the personalized analysis of the expression profiles of the individual
metastasis pairs, because an HMM can utilize such local chromosomal
dependencies to improve the predictions of differentially expressed
genes [[95]32–[96]35]. In more detail, the chromosome-specific
log[2]-ratio gene expression profiles of all patient-specific
metastasis pairs were analyzed by a standard first-order three-state
HMM with state-specific Gaussian emission densities specifically
developed for the analysis of individual expression profiles [[97]33].
This HMM has three states to classify each gene according to its most
likely expression state (Fig. [98]2B): (i) state ’−’ represents genes
with decreased expression in the intracranial metastasis compared to
the corresponding extracranial metastasis, (ii) state ’
[MATH: = :MATH]
’ represents genes with unchanged expression, and (iii) state ’
[MATH: + :MATH]
’ represents genes with increased expression in the intracranial
metastasis compared to the corresponding extracranial metastasis of a
patient-matched metastasis pair. Motivated by the distribution of the
gene expression log[2]-ratios of all patient-matched metastasis pairs
(Additional file [99]2: Figure S2), the initial means of the
state-specific Gaussian emission densities were set to -3, 0, 3 and the
corresponding initial standard deviations were set to 0.5, 1, 0.5 for
the states ’−’, ’
[MATH: = :MATH]
’, and ’
[MATH: + :MATH]
’, respectively. In addition, the initial state probabilities were set
to 0.1 for the states ’−’ and ’
[MATH: + :MATH]
’ and to 0.8 for state ’
[MATH: = :MATH]
’ to account for the fact that the vast majority of genes was unchanged
in their expression. The resulting initial HMM was trained based on all
patient-matched log[2]-ratio gene expression profiles using the
Bayesian Baum-Welch algorithm that enables to integrate prior knowledge
about expected gene expression changes into the training [[100]32,
[101]33]. To realize this, a grid search was done to determine
hyperparameter settings for the state-specific Gaussian emission
densities that led to biologically meaningful state representations of
the trained HMM following the basic strategy developed in [[102]25]. In
more detail, this search was used to identify a hyperparameter setting
that minimized the number of wrongly classified genes (genes with
negative log[2]-ratios assigned to state ’
[MATH: + :MATH]
’ or genes with positive log[2]-ratios assigned to state ’−’) while
maximizing the total number of genes assigned to state ’
[MATH: + :MATH]
’ or ’−’ to obtain a clear separation between genes with reduced and
genes with increased expression. A good separation between reduced,
unchanged, and increased expression levels was obtained by setting the
values of the scale parameter (scaleMeans) to 2500, 1000, 2500 and
those of the shape parameter (shapeSds) to 5000, 10, 5000 for the
states ’−’, ’
[MATH: = :MATH]
’, and ’
[MATH: + :MATH]
’, respectively. All other hyperparameters were kept at their
pre-defined standards. The training of the corresponding initial HMM
took 118 iteration steps and was finished in less than one minute on a
standard laptop (i7-8565U CPU: 1.80 GHz, 8 GB RAM) using the HMM
implementation from [[103]33]. State-posterior decoding was used to
determine for each gene in a log[2]-ratio gene expression profile its
most likely underlying expression state (Additional file [104]10:
Table S5).
Pathway enrichment analysis of differentially expressed genes
Enrichment analysis of differentially expressed genes was done for
known cancer- and immune-relevant pathways from KEGG [[105]36]
including ’Pathways in cancer’ (hsa05200). Corresponding pathway genes
were obtained using KEGGREST version 1.26.0 [[106]37] based on the KEGG
release 103.0. The obtained cancer signaling pathway annotations are
listed in Additional file [107]11: Table S6 and the immune pathways are
provided in Additional file [108]12: Table S7. Based on this, each gene
was annotated using these pathway annotations. Next, each pathway was
tested for a statistical significant enrichment of genes with decreased
or increased expression in each metastasis pair using Fisher’s exact
test (R function fisher.test). Each pathway was considered to be
significantly enriched for decreased or increased expression when its
FDR-adjusted p-value was less then 0.05 (R function p.adjust)
[[109]38].
Determination of differentially expressed genes shared across the majority of
patients
To obtain genes that were altered in their expression in the majority
of patients, a ranking of genes across all patients was done according
to the number of patients in which they showed the same differential
expression state. To realize this, multiple metastasis pairs from the
same patient had to be summarized to count each of the four affected
patients only once (Table [110]1). For these four patients, a majority
vote for the predicted HMM expression states ’
[MATH: + :MATH]
’ or ’−’ was performed for each gene across the multiple metastasis
pairs of the same patient. In that process, a small neglectable
proportion of 0.52% of all genes had an equal number of ’
[MATH: + :MATH]
’ and ’−’ predictions in metastasis pairs of the patients with multiple
pairs and were therefore not further considered in the analysis.
Considering all patients, the genes were ranked according to their
decreasing frequency of decreased expression (number of patients for
which a gene was predicted to have the HMM expression state ’−’) and
another ranking was made for all genes according to their decreasing
frequency of increased expression (number of patients with HMM
expression state ’
[MATH: + :MATH]
’). Additional file [111]13: Table S8 contains both rankings.
Cellular functions of differentially expressed genes shared across the
majority of patients
To characterize major cellular functions, a gene ontology (GO) analysis
was performed considering all differentially expressed genes whose
expression was altered in the same direction in at least 8 of 16
patients (Additional file [112]13: Table S8). This was done using
clusterProfiler [[113]39] with a targeted focus on altered biological
processes. The results of this GO analysis are provided in Additional
file [114]14: Table S9. In addition to the automatic GO enrichment
analysis, a manual but in-depth hand-curated gene function and
literature analysis was performed for the more stringent set of the top
103 genes that either showed increased or decreased expression in
intra- compared to corresponding extracranial metastases in at least 11
of 16 patients (Additional file [115]13: Table S8) using UniProt
[[116]40], GeneCards [[117]41] and Pubmed
([118]https://pubmed.ncbi.nlm.nih.gov/).
Independent validation of differentially expressed genes shared across the
majority of patients considering three related studies
The expression behavior of the top 103 genes was compared to their
expression behavior in three related studies that performed
transcriptome analyses of intra- and extracranial melanoma metastases
[[119]18, [120]22, [121]26]. Therefore, the supplementary table with
the corresponding top differential gene expression candidate gene set
was downloaded for each of this three studies. Each of our 103
candidate genes was then checked for the presence and the direction of
its differential expression in intra- compared to extracranial
metastases in each of these candidate gene sets.
Association of candidate gene expression with patient survival
Expression levels of candidate genes were tested for associations with
patient survival using the melanoma cohort from The Cancer Genome Atlas
(TCGA) [[122]42] that comprises primary and metastatic melanomas. Only
patient samples with available survival information and a known tumor
content of at least 80% were considered. Following the TCGA gene
expression data pre-processing approach in [[123]27], the corresponding
raw gene expression counts of the patients were normalized by cyclic
loess normalization [[124]29]. Next, only genes with more than one
count per million (CPM) reads in at least 50% of the patients were kept
to exclude lowly expressed genes. For each of the 38 candidate genes
from our study that were also measured in the TCGA cohort, the TCGA
patients were divided into two groups: (i) a high expression group
including patients with expression of the specific candidate gene in
the fourth quartile of all expression levels of this gene, and (ii) a
low expression group including patients with expression levels of the
specific candidate gene in the first quartile of all expression levels
of the specific gene across the TCGA patients. Based on this, for each
of the 38 candidate genes differences in survival between both groups
were analyzed by a Kaplan-Meier analysis and tested for statistical
significance by performing a log-rank test using the R package
survminer version 0.4.1. The two corresponding patient groups for each
gene along with the survival information are provided in Additional
file [125]15: Table S10 and the corresponding normalized gene
expression levels are contained in Additional file [126]19: Table S14.
Each gene was considered to be significantly associated with survival
if its FDR-adjusted p-value was less then 0.1 (R function p.adjust)
[[127]38].
Results
Hierarchical clustering of melanoma metastasis expression profiles suggests
the need for a personalized analysis of patient-matched metastasis pairs
Our cohort comprises 16 melanoma patients that all developed an intra-
and an extracranial metastasis during the course of their disease
(Table [128]1). The corresponding metastasis transcriptomes of seven
patients were taken from [[129]21] including all patients for which
patient-matched transcriptomes of intra- and extracranial melanoma
metastasis samples were available. In addition, the transcriptomes of
metastasis pairs of nine additional patients were newly sequenced for
this study. The extracranial metastases of four patients (P04, P08,
P18, P42, Table [130]1) contained histologically distinct regions. The
transcriptomes of these distinct regions from the same metastasis
sample were included as separate patient-specific extracranial samples
to enable the analysis of molecular differences between the potential
subclones within a metastasis and in relation to their patient-matched
intracranial metastasis. Further, the majority of patients (9 of 16)
were not treated before metastasis resection, five of the seven
pre-treated patients were neither treated with BRAF/MEK nor with immune
checkpoint inhibitors, and two patients received targeted therapies
(Additional file [131]18: Table S13). Additional file [132]18:
Table S13 also contains the mutation status of BRAF and NRAS of the
individual metastases, the sex of the patients, and the age of the
patients at the time when the intracranial metastasis was resected.
First, it was reassured that the analyzed melanoma metastases had high
tumor content. This was done by performing a hierarchical clustering of
all initially profiled metastasis samples in a joint analysis together
with transcriptomes of normal samples of different tissues comprising
brain, lymph node, lung, soft tissue and skin (Additional file [133]1:
Figure S1). The metastases and the normal tissues formed separate
disjoint clusters, except for one lymph node metastasis sample of
patient P06 that co-clustered together with the normal tissues.
Therefore, this patient was excluded from the study.
Next, a hierarchical clustering of the transcriptomes of all remaining
metastases was performed in combination with a cluster stability
analysis to characterize global similarities and differences between
all individual metastases (Fig. [134]1). This hierarchical clustering
showed three main subclusters, which were completely stable (AU value =
100). The left subcluster (pink bar) only comprises intracranial and
lymph node metastases. The middle cluster (brown bar) includes all
samples from P08 and lymph node metastasis samples from P42. Further,
it is important to note that P08 was the only patient that had a soft
tissue metastasis. There were no associations with available patient
meta-information that explained the observed co-clustering of
metastasis samples of P08 and P42. The right subcluster (yellow bar)
comprises metastases from six of seven extracranial tissue types (all
except soft tissue).
Fig. 1.
[135]Fig. 1
[136]Open in a new tab
Hierarchical clustering of transcriptomes of all patient-matched
melanoma metastases. Each patient developed an intracranial metastasis
(B, grey) and an extracranial metastasis in the course of its disease.
Extracranial metastases appeared in either lung (Lun, blue), lymph node
(Lym, green), skin (Ski, yellow), liver (Liv, pink), small intestine
(Smi, light green) or soft tissue (Sof, purple). Multiple samples of
the same metastasis were taken if the metastasis showed histologically
different regions. Metastases of the same patient mostly co-cluster
together. The few exceptions from this observation are marked with an
asterisk ’*’. Affiliation with one of the three main clusters is marked
in the lowest color bar (left cluster: pink, middle cluster: brown,
right cluster: yellow). Stability of individual clusters is quantified
by the red AU value, where 100 means that a cluster was perfectly
stable
In addition, the clustering also showed that intra- and extracranial
metastases from the same patient were more similar to each other than
to metastases in the same tissue from other patients. All observed
patient-specific subclusters of intra- and extracranial metastasis
pairs were fully stable (Fig. [137]1, AU value = 100). In more detail,
for 13 of 16 patients, the corresponding intra- and extracranial
metastases all formed completely stable patient-specific subclusters
that represented patient-matched metastases co-clustered together
(Fig. [138]1). For the other three patients, their intra- and
extracranial metastases did not directly cluster together: P13 (brain
and lymph node metastases), P108 (brain and lymph node metastases), and
P42 (brain and lymph node metastases). Despite the fact that these
three patients all developed a lymph node metastases, this observation
cannot be generalized for all lymph node metastases in the cohort. The
other four patients that developed a lymph node metastasis formed
completely stable patient-specific subclusters that contained the
corresponding intracranial metastasis co-clustered together with the
patient-matched lymph node metastasis (P106, P111, P77, P74;
Fig. [139]1).
Generally, such a patient-specific co-clustering of melanoma metastases
has also been observed in closely related studies with gene expression
or DNA-methylation data of patient-matched intra- and extracranial
melanoma metastasis pairs [[140]25, [141]26]. This observation
motivates the necessity for a personalized analysis of the
patient-matched metastasis pairs, which is performed subsequently.
Personalized analysis of the expression behavior of patient-matched
metastasis pairs using a Hidden Markov Model
Due to the observed patient-specific co-clustering of metastases, the
transcriptomes of the patient-matched melanoma metastases were analyzed
in a personalized way to identify for each metastasis pair genes with
increased or reduced expression in the intra- compared to the
extracranial metastasis. To realize this, pair-specific log[2]-ratio
gene expression profiles were considered (Additional file [142]9:
Table S4). Such a pair-specific gene expression log[2]-ratio represents
one of three possible expression states of a gene: (i) a value clearly
smaller than zero suggests decreased expression of the specific gene in
the intra- compared to the extracranial metastasis, (ii) a value of
about zero indicates unchanged expression of a gene between both
metastases, and (iii) a value clearly greater than zero suggests
increased expression of the specific gene in the intra- compared to the
extracranial metastasis.
The considered chromosomal log[2]-ratio gene expression profiles of the
patient-matched melanoma metastases pairs show a clear positive
correlation for the expression levels of genes in close chromosomal
proximity (Fig. [143]2A). This positive autocorrelation of gene
expression changes of neighboring genes motivates the usage of a Hidden
Markov Model (HMM) for a personalized analysis of the patient-matched
pair-specific gene expression profiles. An HMM can utilize such
dependencies between neighboring gene expression levels to make
reliable predictions of the underlying gene expression states. This has
already been demonstrated successfully in different studies for tumor
gene expression profiles (e.g. [[144]32, [145]33]). For this study, a
three-state HMM with state-specific Gaussian emission densities was
trained to enable a personalized analysis of the individual
patient-matched melanoma metastasis pairs (Fig. [146]2B, see Methods
for details). Connected to the interpretation of the defined
gene-specific log[2]-ratios above, the obtained HMM was used to assign
each gene in a patient-matched metastasis pair to its most likely
underlying expression state: (i) decreased expression ’−’, (ii)
unchanged expression ’
[MATH: = :MATH]
’, and (iii) increased expression ’
[MATH: + :MATH]
’ in the intra- compared to extracranial metastasis.
The numbers of predicted genes with decreased or increased expression
in intra- compared to extracranial metastases are shown in Fig. [147]3A
for each patient-matched metastasis pair. Overall, on average 1234
genes with decreased expression (standard deviation 429) and on average
1003 genes with increased expression (standard deviation 435) in
intracranial metastases were found across all metastasis pairs. There
was no general trend towards decreased or increased expression of genes
in intra- compared to extracranial metastases across all
patient-matched metastasis pairs. In more detail, the numbers of
predicted differentially expressed genes varied from pair to pair
(Fig. [148]3A: max/min of 2192/504 genes with increased expression for
P42_BLym-2/P77_BLym; max/min of 2047/693 genes with decreased
expression for P13_BLym/P77_BLym).
Fig. 3.
[149]Fig. 3
[150]Open in a new tab
Overview of differentially expressed genes and top-ranked altered
pathways across all patient-matched metastasis pairs. The x-axis shows
each individual metastasis pair grouped by and color-coded according to
the tissue in which the extracranial metastasis occurred (blue: lymph
node, green: lung, yellow: soft tissue, purple: skin, pink: liver,
cyan: small intestine). A, Absolute number of genes with increased
expression in each intracranial metastasis compared to the
corresponding extracranial metastasis (red bars) and the absolute
number of genes with decreased expression in each intracranial
metastasis compared to the extracranial metastasis (blue bars). B,
Percentage of genes associated with the top six cancer-relevant
signaling pathways that were either increased or decreased expressed in
the intracranial metastasis compared to the corresponding extracranial
metastasis of a patient-matched metastasis pair. C, Percentage of genes
associated with the top four immune-relevant pathways that were either
increased or decreased expressed in the intracranial metastasis
compared to the corresponding extracranial metastasis of a
patient-matched metastasis pair. Significant enrichment of decreased or
increased expressed genes in pathways of subpanels B and C are marked
with ’x’ (FDR-adjusted p-value < 0.05)
Characteristic expression alterations of cancer-relevant signaling pathways
of individual patient-matched metastasis pairs
The HMM was used to predict differentially expressed genes for each
patient-matched melanoma metastasis pair (Fig. [151]3A, Additional file
[152]10: Table S5). These predictions formed the basis for an
individual metastasis pair-specific gene enrichment analysis to analyze
which known cancer-relevant signaling pathways were significantly
affected by gene expression alterations in each metastasis pair
(Additional file [153]3: Figure S3A). Every metastasis pair was
enriched for differential expression of at least one cancer-relevant
signaling pathway. Overall, there were six known cancer-relevant
signaling pathways (cytokine-receptor interaction, calcium signaling,
ECM-receptor interaction, cAMP signaling, Jak-STAT and PI3K/Akt
signaling) that were frequently affected by differential expression in
at least ten metastasis pairs (Fig. [154]3B). More details to specific
pathway overrepresentations of individual metastasis pairs are provided
in Additional file [155]20: Text S1.
Characteristic expression alterations of immune-relevant pathways in
individual patient-matched metastasis pairs
Other interesting pathways for an analysis of differentially expressed
genes between patient-matched intra- and extracranial melanoma
metastases are immune signaling pathways. Therefore, a similar
enrichment analysis was done for the pair-specific differentially
expressed genes predicted by the HMM with a focus on known
cancer-relevant immune signaling pathways (Additional file [156]3:
Figure S3B, Fig. [157]3C). Significant enrichments of differential
expression in immune pathways were almost always only observed for
genes with decreased expression in the intra- compared to the
corresponding extracranial metastasis of the patient-matched metastasis
pairs, except for patient P42, who showed an enrichment of immune
signaling genes with increased expression in the intracranial
metastasis. These enrichments of genes with decreased expression were
not restricted to a specific tissue type in which the extracranial
metastases occurred. Further, no significant enrichments of immune
signaling pathways were found for eight metastasis pairs of seven
patients (P106_BLym, P111_BLym, P107_BLun, P16_BLun, P18_BLun-1,
P18_BLun-2, P111_BLiv, P78_BSmi). Generally, four pathways were most
frequently significantly enriched for differentially expressed genes
across all patient-matched metastasis pairs (Fig. [158]3C: Antigen
processing and presentation, Natural killer cell mediated cytotoxicity,
Th1 and Th2, Th17 pathway). All these pathways showed an enrichment of
genes with decreased expression in the intracranial metastases for the
majority of patients, except for patient P42 which showed an enrichment
of genes with increased expression. Details to specific immune pathway
overrepresentations of individual metastasis pairs are summarized in
Additional file [159]21: Text S2. Patient P42, the only patient with a
significant overexpression of immune signaling, had a NRAS mutation
present in both metastases. Roughly a quarter of the patients that did
not show significant enrichments of immune signaling had a BRAF
mutation. A bit more than a third of the patients that showed
significant downregulations of immune signaling had a BRAF and/or NRAS
mutation (Additional file [160]18: Table S13).
Multiple patient-matched metastasis pairs of histologically different regions
are highly similar to each other
The analyzed melanoma metastases cohort contains four patients with
multiple melanoma metastasis pairs that assign the intracranial
metastasis of a specific patient to histologically different regions in
the corresponding extracranial metastasis of this patient
(Table [161]1: P04, P08, P18, P42). These histologically different
regions within the extracranial metastases were marked by an
experienced board-certified pathologist to enable separate analyses. In
Fig. [162]3, it is clearly noticeable that the metastasis pairs from
the same patient show very similar alteration patterns for enrichment
of differentially expressed genes in signaling and immune pathways.
Therefore, similarities of genome-wide gene expression alterations from
patients with multiple metastasis pairs were analyzed in more detail.
Hierarchical clustering of the patient-matched pair-specific gene
expression log[2]-ratio profiles comparing the intra- to the
corresponding extracranial metastasis showed that metastasis pairs of
patients with multiple marked extracranial regions from histologically
different regions always co-clustered together in a patient-specific
manner (Fig. [163]4A). Further, the multiple metastasis pairs from the
same patients showed a significantly greater overlap of predicted gene
expression states than non-patient-matched intra- versus extracranial
metastasis pairs (Fig. [164]4B, two-sided Mann–Whitney-U-test: p <
0.003, median 86.3% vs. 78.7% overlap, Additional file [165]16:
Table S11). Thus, despite histologically different regions in some
extracranial metastases, the corresponding metastasis pairs of patients
with multiple extracranial samples still showed highly patient-specific
expression profiles. These expression profiles are more similar to each
other than to expression profiles of patient-matched pairs or
non-patient-matched pairs from other patients.
Fig. 4.
[166]Fig. 4
[167]Open in a new tab
Similarities of gene expression alterations of metastasis pairs from
the same and different patients. A, Hierarchical cluster dendrogram
representing the similarities of genome-wide log[2]-ratio gene
expression profiles comparing the intra- to the corresponding
patient-matched extracranial metastasis sample. The color bars
represent the tissue type in which the extracranial metastases occurred
(upper bar, blue: lung; green: lymph node, yellow: skin, purple: soft
tissue, pink: liver, light green: small intestine) and highlight the
patient number by a color code (lower bar). Multiple metastasis pairs
from the same patient are labeled with a number at the end of the
corresponding identifier of the metastasis pair (see Table [168]1).
Multiple metastasis pairs of each patient formed distinct clusters.
Stability of individual clusters is quantified by the red AU value,
where 100 means that a cluster was perfectly stable. B, Overlaps of the
HMM-based predictions of differentially expressed genes for all
pairwise comparisons of either two metastasis pairs from different
patients (left boxplot) or for two metastasis pairs from the same
patient (right boxplot). Overlap of HMM predictions are significantly
different between pairs from different and same patients
(Mann–Whitney-U-test: p < 0.003)
Gene ontology analysis identifies biological processes frequently affected by
decreased and increased expression of genes in patient-matched metastasis
pairs
Next, genes that shared the same differential expression across
multiple patients were determined. Since metastasis pairs of the same
patient showed very similar expression behavior (Figs. [169]3 and
[170]4), HMM-based gene expression state predictions from the
metastasis pairs of the same patient were summarized to one expression
state per gene (see Methods). This allowed to include each patient only
once in the ranking of genes to avoid patient-specific biases. The
genes were separately ranked according to their number of increased
expression state predictions (state ’
[MATH: + :MATH]
’ predicted by HMM) across the intracranial metastases of the patients
and according to their number of decreased expression state predictions
(state ’−’ predicted by HMM). Both rankings are provided in Additional
file [171]13: Table S8. In general, more genes with decreased
expression than genes with increased expression were observed for all
patient cutoffs (Fig. [172]5A). Only one gene, CCL19, showed decreased
expression in intra- compared to extracranial metastases in all 16
patients. CCL19 is a cytokine involved in immunoregulatory and
inflammatory processes indicating that these processes are potentially
downregulated in intracranial metastases. None of the genes shared
increased expression in the intra- compared to the corresponding
extracranial metastasis in 15 or more patients, but there were two
genes, ITIH2 and GAP43, which showed increased expression in 14
patients. ITIH2 contributes to stability of the extracellular matrix
and GAP43 is a growth factor highly expressed in neuronal growth.
Fig. 5.
[173]Fig. 5
[174]Open in a new tab
Overlap of differentially expressed genes between melanoma patients and
gene ontology (GO) enrichment analysis of top differentially expressed
genes. A, Number of candidate genes that were shared across multiple
patients. The x-axis shows the number of patients that commonly share
the corresponding differentially expressed genes. The y-axis shows the
total number of shared differentially expressed genes (shared genes
with decreased (blue)/increased (red) expression in patient-matched
intra- compared to extracranial metastasis pairs). B, GO analysis of
biological processes of the 242 candidate genes with increased
expression in intracranial metastases compared to their corresponding
extracranial metastases that were shared across at least eight or more
patients. Colors of the bars mark biological process groups (red:
synaptic processes, orange: brain-specific cell development, dark
brown: neuronal processes, light brown: no specific group). C, GO
analysis of biological processes of the 459 candidate genes with
decreased expression in intracranial metastases compared to their
corresponding extracranial metastases that were shared across at least
eight or more patients. Colors of the bars mark biological process
groups (purple: immune response, green: activation of immune system,
darkgreen: chemokine-related processes, blue: immune-cell development,
light blue: other immune-related processes)
To further analyze biological processes of the top differentially
expressed genes that shared the same expression state in at least 50%
(8 of 16) or more of all patients, the corresponding 242 genes with
increased expression and the corresponding 459 genes with decreased
expression genes were considered for a gene ontology enrichment
analysis [[175]39] (Additional file [176]14: Table S9).
Interestingly, the 242 genes that showed increased expression in
patient-matched intra- compared to extracranial melanoma metastasis
were frequently involved in synaptic processes (e.g. regulation of
trans-synaptic signaling, synapse organization), brain-specific cell
development (e.g. forebrain development, gliogenesis, astrocyte
differentiation) and neuronal processes (e.g. sensory perception,
calcium ion-regulated exocytosis of neurotransmitter, adult behavior)
(Fig. [177]5B). The observation of such a brain-like phenotype for
intracranial melanoma metastases is supported by other related studies
including single-cell RNA sequencing [[178]22, [179]43] and was also
found in our corresponding previous genome-wide DNA-methylation
analysis of the patient-matched metastasis pairs [[180]25]. An
additional analysis of the gene expression levels of the genes
associated with the brain-like phenotype by directly comparing our
intracranial melanoma metastases to our normal brain tissues showed
significantly lower expression for 63.5% and significantly increased
expression for 6% of the genes in the intracranial metastases, whereas
30.5% of the genes were expressed at the same level like in normal
brain tissues (Additional file [181]4: Figure S4). This indicates that
the observed brain-like phenotype is potentially jointly driven by
normal cells in the metastases microenvironment and tumor cells of the
intracranial metastases.
Further, the 459 genes that showed decreased expression in
patient-matched intra- compared to extracranial metastasis were
frequently involved in immune responses (e.g. humoral immune response,
adaptive immune response, positive regulation of immune response),
immune-cell development (e.g. leukocyte differentiation) and
chemokine-related processes (e.g. chemokine-mediated signaling,
response to chemokine) (Fig. [182]5C). This is also supported by a
related study [[183]26], which identified significant immunosuppression
in intracranial melanoma metastases.
Overall, this suggests a down-regulation of immune-related processes
and an up-regulation of processes involved in the establishment of a
brain-like phenotype in patient-matched intra- compared to extracranial
metastasis pairs.
In-depth analysis of most frequently differentially expressed genes in
patient-matched intra- compared to extracranial metastasis pairs
The identification of genes that are altered in the same manner in
their expression between intra- and extracranial melanoma metastases is
an important step to identify potential key genes that contribute to
the establishment of molecular differences between both metastasis
types. Therefore, a more stringent top candidate gene set was created
that consisted of genes that were differentially expressed in intra-
compared to extracranial metastasis in 11 or more patients. This
resulted in 103 genes comprising 75 genes with decreased and 28 genes
with increased expression in intracranial metastases. Most of these 103
candidate genes showed a homogeneous expression behavior across all or
the majority of metastasis pairs and were further manually assigned to
more specific cancer-relevant functional categories by an in-depth gene
function and literature analysis (Fig. [184]6). The expression levels
of the candidate genes in intracranial metastases were further compared
to normal brain tissues to provide hints how tumor cells of the
intracranial metastases and normal cells of the microenvironment may
express individual genes (Additional file [185]5: Figure S5).
Fig. 6.
[186]Fig. 6
[187]Open in a new tab
Heatmap of gene expression differences of most frequently altered genes
between patient-matched intra- and extracranial melanoma metastasis
pairs. The heatmap shows the expression behavior of the 28 genes with
increased and the 75 genes with decreased expression state predictions
by the HMM that were altered in at least 11 or more patients. A color
gradient from blue over white to red is used to visualize the
log[2]-ratios of the gene expression levels of the intra- compared to
the corresponding extracranial metastasis of each patient-matched
metastasis pair. The columns of the heatmap represent the individual
metastasis pairs with labels color-coded by the tissue in which their
extracranial metastasis occurred (blue: lung, green: lymph node,
purple: soft tissue, yellow: skin, pink: liver, cyan: small intestine).
The rows of the heatmap represent the individual candidate genes. The
annotation block on the left side of the heatmap shows general gene
function groups of individual genes and the expression behavior of each
gene in three related studies (Fischer et al. [[188]26], Chen et al.
[[189]18] and Biermann et al. [[190]22])
Considering the 75 genes with decreased expression in intracranial
metastasis in relation to their functional annotations (Fig. [191]6,
lower blue block), several of these genes (26/75) are involved in
immune response (e.g. CCL19, CLEC10A, CD8B, CD79A). Ten genes are
involved in cell growth (FGF7, RSPO3, TNFSF11, EGFL6, SULF1, ADRA1A,
GREM1, FCRL1, IGF1, ANGPTL1), and another nine genes are involved in
signal transduction (GPR68, NPY1R, LOXHD1, P2RY10, RGS13, CILP, STAC,
P2RY14, GPBAR1).
Considering the 28 genes with increased expression in intracranial
metastasis in relation to their functional annotations (Fig. [192]6,
upper red block), seven of them are involved in cell growth (GAP43,
HEPACAM, MAPK4, PPBP, LGI1, TUBB1, SCG3). Six genes are involved in
cellular transport processes (PMP2, SLC38A11, SLC52A3) frequently
including ion transport (GFAP, SLC24A, ATP1A2). Three genes are
involved in cell adhesion (MOG, CNTN2, GLDN), and three other genes
(11%) are involved in metabolism (ITIH2, CYP4F11, CYP4F3).
Comparison of the most frequently differentially expressed genes to three
related studies
In addition, we compared our candidate gene set of 103 genes to the
candidate gene sets of three closely related melanoma metastasis
expression studies that compared intra- to extracranial metastases (see
Methods for details, Fig. [193]6, Additional file [194]17: Table S12).
There was a significant overlap of our candidate gene set with the
candidate gene sets of the three related studies (Fischer et al.
[[195]26]: 38% overlap, p < 0.001; Chen et al. [[196]18]: 3.8% overlap,
p
[MATH: = :MATH]
0.001; Biermann et al. [[197]22]: 47% overlap, p < 0.001; Fisher’s
exact test). Further, the direction of expression changes between
intra- and extracranial metastases were the same for all overlapping
genes with the related bulk gene expression studies from Fischer et al.
[[198]26] and Chen et al. [[199]18], but differed for some genes for
the single-cell gene expression study by Biermann et al. [[200]22].
The genes that were consistently found in multiple independent studies
are potentially the most relevant candidate genes for future
experimental studies. Therefore, the eight genes (MEOX2, FGF7, DPT,
LAMP3, CILP, TMEM119, PMP2, GLDN) that were predicted as differentially
expressed between intra- and extracranial metastases in our and at
least two of the three related studies (Fig. [201]6) were summarized in
Table [202]2 to provide an overview of their expression behavior,
biological functions, and associations with cancer.
Table 2.
Summary of in-depth literature analysis of gene candidates that were
also found to be differentially expressed in at least two of the
closely related melanoma studies
Gene ’
[MATH: + :MATH]
’ ’−’ Support Biological function Cancer association References