RNA expression and risk of venous thromboembolism in lung cancer.

BACKGROUND: The propensity to develop venous thromboembolism (VTE) on the basis of individual tumor biological features remains unknown.
OBJECTIVES: We conducted a whole transcriptome RNA sequencing strategy, focusing on a single cancer type (lung cancer), to identify biomarkers of cancer-associated VTE.
METHODS: Twelve propensity-matched patients, 6 each with or without VTE, were identified from a prospective institutional review board-approved registry at the Cleveland Clinic with available tissue from surgical excision of a primary lung mass between 2010 and 2015. Patients were propensity matched based on age, sex, race, history of prior cancer, date of cancer diagnosis, stage, histology, number of lines of chemotherapy, and length of follow-up. RNA sequencing was performed on tumor tissue, and gene set enrichment analysis (GSEA) was performed on differentially expressed genes.
RESULTS: We identified 1037 genes with differential expression. In patients with VTE, 869 genes were overexpressed and 168 were underexpressed compared to patients without VTE. Of these, 276 overexpressed and 35 underexpressed were significantly different (Q < 0.05). GSEA revealed upregulation of genes in complement, inflammation, and KRAS signaling pathways in tumors from patients with VTE.
CONCLUSIONS: These differentially expressed genes and associated pathways provide biologic insights into cancer-associated VTE and may provide insignts to develop new risk stratification schemes, prevention, or treatment strategies.

Better biomarkers are needed to predict cancer‐associated thrombosis.

We identified RNA expression profiles in patients with lung cancer with and without thrombosis.

Identified genes and pathways involved in inflammation, KRAS, and innate immunity associated with thrombosis in patients with lung cancer.

These differentially expressed genes may promote development of new risk stratification and therapeutic strategies.

INTRODUCTION

The risk of venous thromboembolism (VTE) is increased in patients with cancer, by at least 4‐ to 7‐fold.1, 2 VTE significantly contributes to morbidity and a worse overall prognosis in patients with cancer.3 VTE incidences ranging from 2% to 14% have been reported in different populations of patients with cancer, with the highest VTE risk identified in patients with pancreatic (~11%), lung (~8%), and stomach cancer (~8%).1, 2, 3, 4 Patients with lung cancer have an elevated risk of VTE compared to other malignancies, with incidence rates ranging from 7% to 15%.4, 5, 6 Various clinical predictors for increased VTE risk have been identified; however, further understanding of thrombogenic biomarkers is needed.

The pathophysiology of cancer‐associated VTE involves complex interactions among the tumor (oncogenes and proteins), intracellular signaling pathways, coagulation system, and anticancer treatment.7, 8, 9 A validated risk model for VTE has been developed based on data from a large, prospective US national cohort study of patients initiating a new chemotherapy regimen for solid tumors and lymphoma.10 This risk model has been validated by independent investigators across a range of solid tumors11, 12, 13 and has been incorporated into clinical guidelines for VTE risk assessment.14, 15, 16

Although validated and utilized among a range of solid tumors, this risk model has not been consistent in stratifying lung cancer patients for VTE risk,17, 18 especially in the era of targeted therapies and improved outcomes. Studies assessing epidermal growth factor receptor mutational status have demonstrated that positive status confers an almost 2‐fold increased risk of VTE.17 However, whether this association is related to the mutational status or the active intervention is unknown. Most recently, a global, prospective study (CANTARISK) showed that the risk model was not significantly associated with VTE in patients with non–small cell lung cancer.19

Genomic predictors are widely used in oncology to inform prognosis and predict outcomes.20, 21, 22 For example, Oncotype DX uses a 21‐gene recurrence score to quantify distant recurrence risk and chemotherapy benefit in patients with breast cancer.22 Although various clinical predictors of VTE risk have been identified, the identification of biomarkers, such as genomic predictors with higher sensitivity and specificity are needed to better discriminate patients at increased risk of VTE. Herein, we use an RNA sequencing strategy of lung cancers to nominate biomarkers of cancer‐associated VTE.

METHODS

Identification of cohort and data collection

We evaluated patients with lung cancer from a prospective institutional review board–approved registry at the Cleveland Clinic with available tissue from surgical excision of a primary lung mass between 2010 and 2015. We propensity‐matched patients on the basis of whether they developed VTE. VTE was defined as pulmonary embolism (PE), deep vein thrombosis (DVT), or catheter‐associated thrombus from 3 months prior to diagnosis to 24 months after diagnosis. Nine patients with VTE were propensity‐matched with 9 patients without VTE, with a score difference of <5%. Patients were matched based on the following criteria: age at cancer diagnosis, sex, race, history of prior cancer, date of lung cancer diagnosis, stage, histology (classified as adenocarcinoma, squamous cell carcinoma, mixed small cell cancer/large cell neuroendocrine, or neuroendocrine tumor), number of lines of chemotherapy, and length of follow‐up. A logistic regression model including these variables was used to calculate the probability of VTE. This probability, or propensity score, was used to identify patients with VTE and without VTE with similar propensity scores, which resulted in the two groups being similar with respect to the matching characteristics. Patients were staged surgically if treated with initial surgery, or otherwise staged clinically per the American Joint Committee on Cancer 7th edition.23

Sample preparation for RNA sequencing

Tumors were isolated from surgical specimens. Representative histological sections demonstrated >90% tumor purity in each sample tested. Nevertheless, contributions from stromal and immune cells likely contributed to the sequencing reads. RNA was collected from fresh frozen tumor samples using Allprep kit (Qiagen, Hilden, Germany) and subjected to quality control. Samples were checked using BioA, and libraries were prepared using TruSeq Total Stranded RNA‐RiboZero kit (Illumina, San Diego, CA, USA). Each library was quantified using Qubit and then pooled together. The library pool was quantified using quantitative polymerase chain reaction. The pool was then diluted to 4 nM and loaded on HiSeq at 12 pM. The run was performed using a full flow cell (2 lanes) PE 100; 5% PhiX was spiked in to check run quality (Q > 30 = 96%). Six tumor samples (3 with VTE and 3 without VTE) had insufficient high‐quality RNA to proceed with sequencing analysis, due to tissue viability.

RNA sequencing

RNA was sequenced by HiSeq 2500 (Illumina, San Diego, CA, USA), resulting in paired 100‐nt reads. RNA reads were aligned to the hg19 genome assembly using MapSplice. Gene expression was quantified for transcript models corresponding to TCGA GAF2.13, using RSEM and normalized within sample to a fixed upper quartile. For gene‐level analyses, expression values of 0 were set to the overall minimum value, and all data were log2 transformed. The log2 fold change and adjusted P and Q values (using Benjamini‐Hochberg procedure) were calculated using linear models in combination with moderated t statistic.

Single‐sample GSEA and the information‐based association metric

From 23 710 genes, a total of 1037 genes had differential gene expression. GSEA was performed on these differentially expressed genes. Association between GSEA profiles for each gene set and binary categorical variables were determined using an information‐based similarity metric (RNMI). This quantity is obtained by estimating the differential mutual information24 between the VTE profile and each of the gene sets’ GSEA profiles and then normalizing and rescaling it so that it is defined from 0 (no association) to 1 (perfect association). We estimate the differential mutual information using a kernel‐based method25 which places a Gaussian density centered at each data point and a width determined by a biased cross‐validation estimate.25, 26 The mutual information is then normalized27 using the joint entropy to better account for the intrinsic differences in entropy associated with each single‐sample GSEA profile. Nominal P values for the information‐based association metric scores between the gene sets/pathways and response scores were estimated using an empirical permutation test.

RESULTS AND DISCUSSION

The characteristics of the study population are presented in Table 1. The final study population included 12 patients with lung tissue RNA samples meeting criteria for quality control (6 each with or without VTE) with a median age of 67 years (range, 48‐80), predominantly white (n = 10) males (n = 8) with stage IV disease (n = 7). Histology included 7 of 12 patients with adenocarcinoma, 3 squamous and 2 with mixed/other histology. Eight patients received first‐line chemotherapy with carboplatin doublet therapy. With respect to second‐line treatment, 2 patients received chemotherapy alone, 1 received chemotherapy plus immunotherapy, and 1 received immunotherapy alone. Eight patients (67%) did not receive second‐line treatment, and 9 patients (75%) did not receive third‐line treatment. Median follow‐up was 18 months (range, 1‐69).

Table 1

Characteristics of propensity‐matched patient cohort

Characteristics	All patients, N = 12 (%)	Patients with VTE, n = 6 (%)	Patients without VTE, n = 6 (%)
Age, y
Median (range)	67, (48‐80)	68	66
Gender
Male	8 (66.7)	4 (66.7)	4 (66.7)
Female	4 (33.3)	4 (66.7)	4 (66.7)
Race
White	10 (83.3)	4 (66.7)	6 (100)
Black	2 (16.7)	2 (33.3)
History of coronary artery diseasea	8 (75)	5 (83.3)	3 (66.7)
Prior history of VTE	1 (8.3)	1 (16.7)
Overall stage
IA	1 (8.3)	1 (16.7)
IIA	1 (8.3)		1 (16.7)
IIIA	3 (25.0)	1 (16.7)	2 (33.3)
IV	7 (58.3)	4 (66.7)	3 (66.7)
Tumor stage
T1b	2 (16.7)	2 (33.3)
T2a	3 (25.0)	2 (33.3)	1 (16.7)
T2b	3 (25.0)	1 (16.7)	2 (33.3)
T3	3 (25.0)	1 (16.7)	2 (33.3)
T4	1 (8.3)		1 (16.7)
Nodal stage
N0	3 (25.0)	1 (16.7)	2 (33.3)
N1	2 (16.7)	2 (33.3)
N2	2 (16.7)		2 (33.3)
N3	4 (33.3)	3 (50)	1 (16.7)
Distant metastasis
M0	5 (41.7)	2 (33.3)	3 (50)
M1b	7 (58.3)	4 (66.7)	3 (50)
Tumor status
Primary	8 (66.7)	5 (83.3)	3 (50)
Recurrence	4 (33.3)	1 (16.7)	3 (50)
Histology
Adenocarcinoma	7 (58.3)	3 (50)	4 (66.7)
Squamous cell carcinoma	3 (25.0)	1 (16.7)	2 (33.3)
Mixedb	1 (8.3)	1 (16.7)
Other (neuroendocrine)	1 (8.3)	1 (16.7)
Treatment received (first line)
Carboplatin doublet therapy	8 (66.7)	4 (66.7)	4 (66.7)
Immunotherapy	1 (8.3)		1 (16.7)
None	3 (25.0)	2 (33.3)	1 (16.7)
Treatment received (second line)
Chemotherapy	3 (25.0)	2 (33.3)	1 (16.7)
Carboplatin doublet therapy			1 (16.7)
Pemetrexed/bevacizumab		1 (16.7)
Irinotecan		1 (16.7)
Immunotherapy (Nivolumab)	1 (8.3)		1 (16.7)
None	8 (66.7)	4 (66.7)	4 (66.7)
Treatment received (third line)
Chemotherapy	1 (8.3)		1 (16.7)
Immunotherapy (Nivolumab)	2 (16.7)	1 (16.7)	1 (16.7)
None	9 (75.0)	5 (83.3)	1 (16.7)
Length of follow‐up, mo
Median (range)	18, (1‐69)	19.5	18

Defined as stable angina, acute coronary syndrome, and sudden cardiac death.

Mixed histology signifies combined small cell lung cancer and large cell neuroendocrine.

Of 6 patients without VTE, 3 had a history of coronary artery disease, and no patients had prior history of VTE. Of 6 patients with VTE, 5 (83%) had a history of coronary artery disease, and 1 (17%) had prior history of VTE. Two (33%) had lower extremity DVT, 3 (50%) had PE, and 1 patient (17%) had an extensive catheter‐associated thrombus in the superior vena cava, extending into the right atrium. Two (33%) VTEs were discovered at diagnosis, 1 before treatment, 1 during treatment, and 2 during treatment break. Five patients developed VTE within 1 year of cancer diagnosis, and 1 patient developed VTE about 1.5 years from diagnosis. Specifically, patients developed VTE at days −5, −62, +152, +340, and +566 from respective date of cancer diagnosis.

We identified 1037 genes with differential gene expression. In patients with VTE, 869 genes were overexpressed and 168 were underexpressed compared to patients without VTE (Figure 1). Of these, 276 genes with statistically significant (Q < 0.05) differential expression were identified, which were composed of 241 genes that were overexpressed in patients with lung cancer with VTE and 35 genes that were underexpressed. When gene expression was limited to those that included log (fold change) of +3 or greater and −3 or lower, a total of 146 genes were identified; 130 of which were overexpressed in patients with lung cancer with VTE and 16 genes that were underexpressed. We further limited genes to include a log (fold change) of +4 or greater and −3 or lower, and to include only those with Q < 0.05. This criterion established a total gene population of 27 genes, 15 of which were overexpressed in patients with lung cancer with VTE, and 12 genes that were underexpressed (Table 2).

Figure 1

Scatter plot highlighting differential gene expression (1037 genes) from RNA sequencing for patients with lung cancer with VTE compared to those without VTE. Values in red denote avidly overexpressed genes, up to 4‐fold. Positive log values (869 genes) represent overexpressed genes, while negative values (168 genes) represent underexpressed genes in patients with VTE. VTE, venous thromboembolism

Table 2

Genes with differential expression in lung cancer patients with VTE

Gene	Full gene name	Log (fold change)	Q value	Function
CECR1	Cat eye syndrome chromosome region, candidate 1	4.9	0.000008	Participates in cell proliferation and differentiation
SYT14	Synaptotagmin 14	4.9	0.00001	Mediates membrane trafficking in synaptic transmission
JSRP1	Junctional sarcoplasmic reticulum protein 1	4.8	0.00002	Excitation‐contraction coupling at the sarcoplasmic reticulum, regulates calcium influx and efflux
MIAT	Myocardial infarction associated transcript	4.8	0.000003	Myocardial infarction transcript
SHC4	Squalene‐hopene cyclase adaptor protein 4	4.6	0.00003	Ras activating pathway
LSAMP	Limbic system associated membrane protein	4.4	0.0002	Mediates selective neuronal growth and axon targeting
ASTN1	Astrotactin 1	4.3	0.0004	Neuronal adhesion molecule
IGLL5	Immunoglobulin lambda‐like polypeptide 5	4.3	0.0004	Protein coding
FAM107A	Family with sequence similarity 107 member A	4.3	0.0004	Associated with tumor development
CR1	Complement C3b/C4b receptor 1	4.2	0.0005	Receptor for activated complement
PEX5L	Peroxisomal biogenesis factor 5 like	4.2	0.0005	Hyperpolarization‐activated cyclic nucleotide‐gated channels
TNFSF8	Tumor necrosis factor superfamily member 8	4.1	0.0008	Induces T cell proliferation or death
PROX1	Prospero homeobox 1	4.1	0.0008	Transcription factor involved in developmental processes such as cell fate determination, gene transcriptional regulation, and progenitor cell regulation
NTM	Neurotrimin	4.1	0.0006	Promotes neurite adhesion
ATP8A2	ATPase phospholipid transporting 8A2	4.0	0.001	Involved in lipid flipping
SYT12	Synaptotagmin 12	−3.0	0.02	Mediates calcium‐dependent regulation of membrane trafficking in synaptic transmission
LMX1B	LIM homeobox transcription factor 1 beta	−3.1	0.02	Transcription factor
MIR31HG	MIR31 host gene	−3.1	0.02	Involved in cellular pluripotency, regulates differentiation of myoblasts and other tissue; repress transcription of genes involved in cell senescence
EPHB6	Ephrin type‐B receptor 6	−3.3	0.009	Cell adhesion and migration
KLHL35	Kelch‐like family member 35	−3.3	0.0002	Protein coding
IRX4	Iroquois homeobox 4	−3.4	0.01	Cell differentiation, heart development, multicellular organism development
TGM1	Transglutaminase 1	−3.4	0.004	Crosslinking of proteins
DSG1	Desmoglein 1	−3.5	0.009	Component of desmosome
TRPV6	Transient receptor potential cation channel subfamily V member 6	−3.5	0.007	Involved in calcium channel function
ALOX12B	Arachidonate 12‐lipoxygenase, 12R type	−3.6	0.005	Conversion of arachidonic acid
VSIG10L	V‐set and immunoglobulin domain containing 10 like	−3.6	0.0001	Protein coding
EN1	Engrailed homeobox 1	−4.1	0.0008	Cell development and differentiation

The most highly overexpressed genes were CECR1 and MIAT. CECR1 (4.9‐fold elevation, Q = 0.000008) regulates cell proliferation and differentiation, while MIAT (4.8‐fold elevation, Q = 0.000003) is linked to cardiac thrombosis. SHC4 (4.6, Q = 0.00003) is a gene that activates both Ras‐dependent and Ras‐independent pathways. We also identified genes involved in inflammatory pathways, which regulate the innate immune system, notably, NLRP14 (3.8, Q = 0.0001), CLNK (3.6, Q = 0.005), IL5RA (3.6, Q = 0.006), XCL1 (3.5, Q = 0.008), IL16 (3.4, Q = 0.002), and CCR6 (3.0, Q = 0.02), as well as genes involved in complement activation, CR1 (4.2, Q = 0.0005).

GSEA showed gene sets associated with innate immunity and inflammatory pathways, which included IL2‐STAT5 signaling, IL6‐JAK‐STAT3 signaling, interferon‐ response, complement, and inflammatory response (Table 3). Conversely, genes associated with apoptosis pathways, myogenesis pathways, and epithelial mesenchymal pathways were underexpressed.

Table 3

Gene set enrichment analysis (GSEA)

Associated gene set pathways	Normalized P value	FDR Q value
IL2_STAT5_SIGNALING	0.11	0.44
ALLOGRAFT_REJECTION	0.03	0.25
IL6_JAK_STAT3_SIGNALING	0.50	0.68
INTERFERON_GAMMA_RESPONSE	0.26	0.64
KRAS_SIGNALING_UP	0.42	0.67
COMPLEMENT	0.59	0.67
INFLAMMATORY_RESPONSE	0.40	0.77
APICAL_JUNCTION	0.85	0.84
SPERMATOGENESIS	0.97	0.98
EPITHELIAL_MESENCHYMAL_TRANSITION	0.76	0.89
MYOGENESIS	0.69	1
COAGULATION	0.70	1
KRAS_SIGNALING_DN	0.40	1
TNFA_SIGNALING_VIA_NFKB	0.48	1
APOPTOSIS	0.08	0.73

FDR, false discovery rate.

Several genes identified in our study are consistent with prior genes identified in arterial thrombosis and inflammatory response, for example, MIAT (myocardial infarction‐associated transcript) and NLRP14 (nucleotide‐binding oligomerization domain, leucine rich repeat and pyrin domain containing 14). Recent studies have shown that MIAT is elevated and associated in patients with acute myocardial infarction (P < .001) in the peripheral blood and overexpressed in atherosclerotic plaques in patients undergoing coronary artery bypass surgery.28, 29, 30 Additionally, serum levels of MIAT were positively associated with percentage of lymphocytes and negatively with neutrophils and platelet count, suggesting an inflammatory component.29 Five patients with VTE and 3 patients without VTE had a history of CAD; MIAT was upregulated only in patients with VTE, suggesting that MIAT may be specific for VTE in our cohort, although confirmation is necessary. Additionally, recent studies have identified gene expression profiles in patients with atherosclerosis that do not possess MIAT.31, 32 Similarly, NLRP14 is a member of NLRP3 genes responsible for activation of the inflammasome complex, which may be associated with thrombotic events.33, 34 These data are consistent with the knowledge that inflammatory states such as inflammatory bowel disease, rheumatoid arthritis, and systemic lupus erythematosus are prothrombotic.35, 36, 37

We identified several gene sets that implicated Kras signaling as associated with a higher risk of developing VTE. This finding is consistent with a study demonstrating metastatic colorectal patients (N = 172) with tumor mutant KRAS status (N = 65) were associated with a 14.5% increased risk of VTE.38 If validated, Kras could become a tumor‐based biomarker for assessing VTE risk; the precise means by which activation of Kras can increase the risk of thrombosis remains unclear.

A recent study identified differential gene expression in patients with colorectal cancer (CRC) with VTE only before CRC diagnosis and around the time (±3 months) of CRC diagnosis.39 This study found 30 differentially expressed genes and performed Ingenuity pathway analysis on the top 10 genes in patients with VTE. This study also found inflammatory pathway associated with VTE and several shared genes with our study: SBSPON, DEFA5, PTPRR, SORBS1. Another study identified similar genes and pathways in patients with no cancer with single and recurrent VTE. Both cohorts had differential expression of immune and inflammatory pathways, specifically interferon‐γ and STAT3, and cancer pathway RAS, which is consistent with our study.40 Of the 36 genes presented in this study, overlap of inflammatory genes LTB and SLC7A11 were identified. An additional study identified 559 and 294 differentially expressed genes in patients with no cancer with recurrent and single VTE, respectively. Of these, 202 upregulated and 58 downregulated genes overlapped between the two groups.41 Leukocyte transendothelial migration and JAK‐STAT3 signaling pathways were associated with recurrent VTE, the latter of which is consistent with our study. In both patients with cancer and patients without cancer with VTE, inflammatory and immune pathways were identified across several studies, suggesting that these pathways are involved in thrombosis.

A limitation of this study is our sample size, which limited the power of this analysis. We were limited by tissue quality, as several tissue samples were of insufficient quality for RNA sequencing. Additionally, stromal cells and immune cells may have contributed to the sequencing reads, although representative histological sections demonstrated >90% tumor purity in each sample tested. Bulk sequencing provides insight to the tumor microenvironment and can provide additional knowledge to the function of immune cells in relation to VTE events. The results should be considered provisional, requiring validation through replication in independent study populations and demonstration of biologic plausibility.

In conclusion, to our knowledge, this study is the first to employ whole transcriptome RNA profiling strategy in lung cancer to identify differentially expressed genes and possible pathways that may guide therapeutic targets for cancer‐associated VTE in the future.

RELATIONSHIP DISCLOSURES

MA reports grants from Siemens Healthcare, and grants and personal fees from Bayer AG, outside the submitted work. KM reports personal fees from Pfizer, outside the submitted work. AAK reports personal fees and nonfinancial support from Janssen, Bayer, Sanofi, Halozyme, Pfizer, Seattle Genetics, Pharmacyclics, and AngioDynamics; personal fees from Parexel, Pharmacyte, Pharmacyclics, Leo Pharma, Trisalus, and Medscape; and grants from Merck, Array, Bristol Myers Squibb, and Leap, outside the submitted work. TAS reports nothing to disclose.

AUTHOR CONTRIBUTIONS

TAS collected data, analyzed data, and wrote manuscript; MEA performed research, contributed analytical tools, and analyzed data; KRM interpreted data and reviewed the manuscript; and AAK designed the research, analyzed data, and revised the manuscript.