Expression profile analysis of long noncoding RNA in HER-2-enriched subtype breast cancer by next-generation sequencing and bioinformatics.

BACKGROUND: Human epidermal growth factor receptor 2 (HER-2)-enriched subtype breast cancer is associated with a more aggressive phenotype and shorter survival time. Long non-coding RNAs (LncRNAs) have essential roles in tumorigenesis and occupy a central place in cancer progression. Notably, few studies have focused on the dysregulation of LncRNAs in the HER-2-enriched subtype breast cancer. In this study, we analyzed the expression profile of LncRNAs and mRNAs in this particular subtype of breast cancer.
METHODS: Seven pairs of HER-2-enriched subtype breast cancer and normal tissue were sequenced. We screened out differently expressed genes and measured the correlation of the expression levels of dysregulated LncRNAs and HER-2 by Pearson's correlation coefficient analysis. Gene ontology analysis and pathway analysis were used to understand the biological roles of these differently expressed genes. Pathway act network and coexpression network were constructed.
RESULTS: More than 1,300 LncRNAs and 2,800 mRNAs, which were significantly differently expressed, were identified. Among these LncRNAs, AFAP1-AS1 was the most dysregulated LncRNA, while ORM2 was the most dysregulated mRNA. LOC100288637 had the highest positive correlation coefficient of 0.93 with HER-2, while RPL13P5 had the highest negative correlation coefficient of -0.87. The pathway act network showed that MAPK signaling pathway, PI3K-Akt signaling pathway, metabolic pathways, cell cycle, and regulation of actin cytoskeleton were highly related with HER-2-enriched subtype breast cancer. Coexpression network recognized LINC00636, LINC01405, ADARB2-AS1, ST8SIA6-AS1, LINC00511, and DPP10-AS1 as core genes.
CONCLUSION: These results analyze the functions of LncRNAs and provide useful information for exploring candidate therapeutic targets and new molecular biomarkers for HER-2-enriched subtype breast cancer.

Background

Breast cancer, one of the most common malignant tumors, is a molecularly heterogeneous disease that includes four major tumor subtypes at least.1–3 Human epidermal growth factor receptor 2 (HER-2)-enriched subtype is characterized by HER-2 over-expression and/or amplification with a lack of hormone receptor. HER-2-enriched subtype breast cancer is associated with a more aggressive phenotype and shorter survival time.4,5 Although HER-2-targeted agents such as trastuzumab have dramatically improved prognosis of patients with HER-2 overexpression, resistance to such agents remained a severe problem.6 Hence, the underlying molecular mechanisms of the malignant phenotype of HER-2-enriched subtype breast cancer need to be elucidated, and new therapeutic targets need further research.

Noncoding RNAs (ncRNAs) are RNA transcripts that do not code for proteins. General conventions divide ncRNAs into two main categories: small ncRNAs <200 bp and long noncoding RNAs (LncRNAs) >200 bp.7 Recent discoveries analyzed the functions of LncRNAs in human cancers, supporting the fact that they had essential roles in tumorigenesis and occupied a central place in cancer progression.8,9

Notably, few studies have focused on the dysregulation of LncRNAs in the HER-2-enriched subtype breast cancer. In this study, we sequenced the expression profile of LncRNAs and messenger RNAs (mRNAs) in HER-2-enriched subtype breast cancer samples and adjacent nontumorous tissue. More than 1,300 LncRNAs and 2,800 mRNAs that were significantly differently expressed were identified. Gene ontology (GO) analysis and pathway analysis were used to understand the biological roles of these differently expressed genes. Pathway act network illustrated the pathways occupying a central place in HER-2-enriched subtype breast cancer. Coexpression network revealed several LncRNAs having important regulation and control ability, which might hopefully work as candidate therapeutic targets and new molecular biomarkers for HER-2-enriched subtype breast cancer.

Materials and methods

Patient samples

Written informed consent was obtained from all patients and the study was approved by the Ethics Committee of the First Affiliated Hospital of Wenzhou Medical University, Wenzhou, Zhejiang, People’s Republic of China. Seven breast cancer patients who received modified radical mastectomy were included in this study. All patients were diagnosed with HER-2-enriched subtype breast cancer by pathology and immunohistochemistry (IHC) after surgery. IHC results showing a 2+ HER-2 level were further tested using fluorescence in situ hybridization (FISH). Patients having a 3+ HER-2 level as tested by IHC or having HER-2 gene amplification as evidenced by FISH along with a lack of both estrogen receptor and progesterone receptor were regarded as having HER-2-enriched subtype breast cancer. Primary breast cancer tissues and their adjacent normal breast tissues were snap-freezed in liquid nitrogen immediately after resection and then stored at −80°C before RNA extraction. Detailed information of all cases in the study is summarized in Table S1.

RNA extraction and sequencing

Total RNA was extracted from tissue samples using the TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s protocol. The cDNA libraries for single-end sequencing were prepared using Ion Total RNA-Seq Kit v2.0 (Life Technologies) according to the manufacturer’s instructions. The cDNA libraries were then processed for the Proton (Life Technologies) sequencing process according to the commercially available protocols.

Mapping and identification of differently expressed genes

Before reads mapping, clean reads were obtained from the raw reads by removing the adaptor sequences, reads with >5% ambiguous bases, and low-quality reads. The clean reads were then aligned to the human genome (version: GRCH37) using the MapSplice program (v2.1.6, University of Kentucky, Lexington, KY, USA). We applied EBseq algorithm to screen out the differently expressed genes using the following criteria: 1) fold change (FC) >2 for up- or downregulation and 2) false discovery rate (FDR) <0.05. A volcano plot was drawn by the R based on the differently expressed gene analysis and the color was determined by the filtering criteria. Pearson’s correlation coefficient was calculated to measure the linear correlation of the expression levels of LncRNAs and HER-2.

GO analysis and pathway analysis

GO analysis was carried out to facilitate elucidating the biological implications of unique genes in the significant or representative profiles of the differently expressed genes in the experiment.10 We downloaded the GO annotations from Gene Ontology (http://www.geneontology.org/). Fisher’s exact test was applied to identify the significant GO terms, and FDR was utilized to correct the P-values. Pathway analysis was utilized to find out the significant pathways of the differently expressed genes according to Kyoto Encyclopedia of Genes and Genomes database (KEGG). We turned to the Fisher’s exact test to select the significant pathways, and the threshold of significance was defined by P-value and FDR.11

Construction of pathway act network and coexpression network

We chose genes enriched in significant biological pathways (P<0.05) and used Cytoscape (V3.2.0; Institute of Systems Biology, Seattle, WA, USA) to construct a pathway act network for graphical representations of central pathways.12 A coexpression network based on the correlation analysis between the differently expressed LncRNAs and mRNAs associated with cancer was constructed.13 LncRNAs and mRNAs with Pearson’s correlation coefficients >0.99 were used to construct the network. K-core scoring was used to identify core regulatory genes in networks. A k-core of a given gene indicated its hub or nodal status with connections to other genes in a network.14 Accordingly, a higher k-core score meant a more central location of a gene within a network. Core regulatory factors were determined by the k-core difference (difk-core) between two groups of samples.

Results

LncRNAs and mRNAs expression profiles in HER-2-enriched subtype breast cancer

Volcano plot was used for assessing gene expression variation between the HER-2-enriched subtype breast cancer and adjacent normal tissue (Figure 1A). Genes with FC >2 for both up- or downregulation and FDR <0.05 were identified as significantly differently expressed. Compared to the normal breast tissues, a total of 1,382 LncRNAs displayed differential expression in tumor tissues, including 722 upregulated LncRNAs and 660 downregulated LncRNAs (Figure 1B). We found 2,879 differently expressed mRNAs, of which 722 were upregulated and were 2,157 downregulated (Figure 1C). Hierarchical clustering showed systematic variations in the expression of LncRNAs and mRNAs in the HER-2-enriched subtype breast cancer samples (Figure 1D and E). The most dysregulated LncRNAs and mRNAs are shown in Tables 1 and 2. Among these significantly differently expressed LncR-NAs, AFAP1-AS1 (GeneID 84740) was the most dysregulated LncRNA with an FC of 9.79, while ORM2 (GeneID 5005) was the most dysregulated mRNA with an FC of 9.85.

Figure 1

Summary of RNA sequencing results.

Notes: (A) Volcano plots. The negative log of FDR (base 10) is plotted on the Y-axis, and the log of the FC (base 2) is plotted on the X-axis. The red points on this graph represent LncRNAs and mRNAs that are significantly differently expressed in HER-2-enriched subtype breast cancer (FC >2 and FDR <0.05), the blue points represent LncRNAs and mRNAs with FC >2 and FDR >0.05, the green points represent LncRNAs and mRNAs with FC <2 and FDR <0.05, and the gray points represent LncRNAs and mRNAs with FC <2 and FDR >0.05. (B, C) Pie charts show the distribution of dysregulated LncRNAs and mRNAs according to RNA sequencing in seven pairs of HER-2-enriched subtype breast cancer. (D, E) Differentially expressed LncRNAs and mRNAs in tumor and adjacent normal tissues are analyzed using hierarchical clustering. Each row represents a single LncRNA or mRNA and each column represents one tissue sample. Red indicates high relative expression and green indicates low relative expression.

Abbreviations: FDR, false discovery rate; FC, fold change; HER-2, human epidermal receptor-2; LncRNA, long noncoding RNA; mRNA, messenger RNA.

Table 1

The most dysregulated LncRNAs

Gene symbol	Type of gene	Log2FC	FDR	Status
AFAP1-AS1	ncRNA	9.794139497	4.22995E–14	Up
LOC101926960	ncRNA	−8.701495579	0	Down
LOC101928141	ncRNA	−8.174040861	1.33556E–11	Down
LOC101928958	ncRNA	7.613666463	8.80878E–08	Up
LOC101929722	ncRNA	−7.517970448	8.60645E–13	Down
LOC101927630	ncRNA	7.422258015	8.49117E–08	Up
LOC100507651	ncRNA	7.394177076	1.30833E–09	Up
LOC150622	ncRNA	7.36082617	0	Up
LOC152225	ncRNA	7.205420187	5.39482E–11	Up
DPP10-AS1	ncRNA	7.033724274	7.91208E–08	Up
ADARB2-AS1	ncRNA	6.985229437	4.832E–08	Up
LOC101927344	ncRNA	6.803228275	5.25509E–11	Up
LOC101929440	ncRNA	6.502054052	2.37481E–07	Up
TRHDE-AS1	ncRNA	−6.149875004	6.66134E–16	Down
LINC00922	ncRNA	6.14652324	3.42103E–11	Up
LOC100129931	ncRNA	6.111626135	1.13635E–08	Up
LOC101928219	ncRNA	−6.019679259	7.43652E–06	Down
LINC00470	ncRNA	5.885406732	1.78791E–05	Up
LOC100996635	ncRNA	−5.793089791	4.33931E–12	Down
CASC9	ncRNA	5.790007128	3.34011E–07	Up
ST8SIA6-AS1	ncRNA	5.762996171	1.47072E–10	Up
LOC100507600	ncRNA	5.65612642	8.95442E–06	Up
PGM5-AS1	ncRNA	−5.580468913	1.1483E–11	Down
LINC01187	ncRNA	5.520366066	6.08022E–05	Up
LOC101928845	ncRNA	5.516592549	5.64646E–07	Up
LINC00032	ncRNA	−5.515257945	0.000226728	Down
WT1-AS	ncRNA	5.515183303	3.82263E–07	Up
LOC286442	ncRNA	5.413074265	8.82323E–07	Up
FEZF1-AS1	ncRNA	5.387253346	1.88E–07	Up
LOC101929691	ncRNA	5.350563789	2.04E–05	Up

Abbreviations: ncRNA, noncoding RNA; LncRNAs, long noncoding RNAs; FC, fold change; FDR, false discovery rate.

Table 2

The most dysregulated mRNAs

Gene symbol	Type of gene	Log2FC	FDR	Status
ORM2	Protein-coding gene	9.852721187	0	Up
MMP1	Protein-coding gene	9.809536292	0	Up
TMPRSS4	Protein-coding gene	9.687249493	0	Up
GCNT3	Protein-coding gene	9.629566048	4.44E–16	Up
UGT1A7	Protein-coding gene	9.450805442	1.87E–11	Up
CSAG2	Protein-coding gene	9.432817101	3.14E–12	Up
MYOC	Protein-coding gene	−9.320905139	0	Down
S100A7	Protein-coding gene	9.240768658	0	Up
IVL	Protein-coding gene	9.099686214	0	Up
UGT1A8	Protein-coding gene	9.090895458	1.14E–10	Up
CXorf61	Protein-coding gene	9.063258798	2.96E–10	Up
PAX7	Protein-coding gene	8.973405465	3.33E–16	Up
LOC101929578	Protein-coding gene	8.970975966	2.22E–16	Up
PGC	Protein-coding gene	8.939546792	1.55E–15	Up
INA	Protein-coding gene	8.932432782	3.44E–12	Up
MYH6	Protein-coding gene	8.846381491	2.66E–15	Up
SLC5A5	Protein-coding gene	8.843356207	1.05E–13	Up
TBX10	Protein-coding gene	8.754465826	9.77E–14	Up
CEACAM6	Protein-coding gene	8.67620907	6.66E–16	Up
PGA3	Protein-coding gene	−8.579918705	1.36E–12	Down
CST1	Protein-coding gene	8.497095625	0	Up
ORM1	Protein-coding gene	8.471380655	1.81E–14	Up
PGA5	Protein-coding gene	−8.370889976	6.91E–11	Down
LHFPL5	Protein-coding gene	8.294381419	4.07E–13	Up
FGG	Protein-coding gene	8.204193271	9.77E–15	Up
MMP13	Protein-coding gene	8.201153001	0	Up
CSAG3	Protein-coding gene	8.174130048	4.66E–09	Up
S100P	Protein-coding gene	7.983896233	1.07E–14	Up
ARHGAP36	Protein-coding gene	−7.962901754	1.14E–14	Down
DCD	Protein-coding gene	7.939785332	1.56E–13	Up

Abbreviations: FDR, false discovery rate; mRNAs, messenger RNAs; FC, fold change.

Expression correlation between LncRNAs and HER-2

We performed Pearson’s correlation coefficient analysis to measure the correlation of the expression levels between the 1,382 dysregulated LncRNAs and HER-2 (Table S2). LOC100288637 (GeneID 100288637) had the highest positive correlation coefficient of 0.93 with HER-2, while RPL13P5 (GeneID 283345) had the highest negative correlation coefficient of −0.87.

Function analysis of differently expressed genes

GO and KEGG pathway analysis of differently expressed mRNAs provided a measure of the critical function. We included all differently expressed mRNAs for GO analysis and found that the most enriched GO was correlation with cell adhesion in the GO biological process analysis (Figure 2A). Meanwhile, the majority of the genes were proven to be related to the extracellular region in the cellular component analysis and calcium ion binding in the molecular function (Figure 2B and C). In the KEGG pathway analysis, the dysregulated mRNAs were found to be enriched in 41 pathways (Table S3). The most enriched pathway included cytokine–cytokine receptor interaction, steroid hormone biosynthesis, and protein digestion and absorption (Figure 2D).

Figure 2

(A–C) The top 15 GO terms associated with biological process, cellular component, and molecular function of differently expressed mRNAs are shown. (D) KEGG pathway analysis for differently expressed mRNAs.

Note: The negative log of P-value (base 2) is plotted on the X-axis.

Abbreviations: GO, Gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes database; mRNAs, messenger RNAs.

Pathway act network

A pathway act network was constructed using 41 significantly enriched pathways to illustrate the key pathways in the process of HER-2-enriched subtype breast cancer (Figure 3). Our finding indicated that the MAPK signaling pathway, PI3K-Akt signaling pathway, metabolic pathways, and cell cycle and regulation of actin cytoskeleton were considered to be the most important nodes in the net because the component exchanges with other pathways were strongly dependent on their existence.

Figure 3

Pathway act network.

Notes: A node represents a signaling pathway. The node color is correlated with pathway expression pattern. Red indicates that the signaling pathway is activated, while green indicates that the signaling pathway is suppressed. Yellow indicates that the genes included in the corresponding signaling pathway are both upregulated and downregulated. Lines represent interactive relationship between signaling pathways.

Coexpression network

Coexpression network was constructed for HER-2 enriched subtype breast cancer group and adjacent normal tissue group, using expressed LncRNAs and mRNAs that were significantly different. (Figure 4A and B). The coexpression network in the tumor group comprised 649 network nodes and 2,169 connections, including 27 that were negatively connected and 2,142 that were positively connected. The network in the normal group contains 729 network nodes and 2,341 connections, including 1,785 pairs that presented as positive, and 556 pairs that presented as negative. This result demonstrated that there were obviously different coexpression patterns between the tumor group and the normal group. LINC00636 (GeneID 285205) and LINC01405 (GeneID 100131138) had the highest k-core score in the tumor group. ADARB2-AS1 (GeneID 642394), ST8SIA6-AS1 (GeneID 100128098), LINC00511 (GeneID 400619), and DPP10-AS1 (GeneID 389023) had the highest k-core score in the normal group. Coexpressed genes with higher difk-core scoring were considered to have important regulation and control ability. Our data indicated that the aforementioned six LncRNAs, having highest difk-core scoring, were the central genes within the coexpression network (Figure 4C).

Figure 4

Coexpression networks.

Notes: (A, B) Coexpression networks in the tumor tissue group and adjacent normal tissue group. The lines indicate a correlative relationship. The circles represent mRNAs and arrows represent LncRNAs. Different colors represent the corresponding k-core scoring, and we highlight the highest k-core-scoring hub. (C) Core LncRNAs in HER-2-enriched subtype breast cancer are listed by the difference k-core scoring between tumor group and normal group.

Abbreviations: mRNAs, messenger RNAs; LncRNAs, long noncoding RNAs.

Discussion

With the emergence of studies focusing on the functional attributes of LncRNAs, it has been revealed that LncRNAs may contribute significantly to physiological processes as well as pathological conditions. Some LncRNAs may act as tumor suppressor genes,15–17 whereas others may be defined as oncogenes.18–20 However, LncRNAs have just begun to be understood, and the majority of them have yet to be researched. Xu et al21 analyzed LncRNAs’ expression profile in human breast cancers by using microarrays. Then, Shen et al22 researched LncRNAs’ expression in triple-negative breast cancer and found that a lot of LncRNAs were differently expressed. Breast cancer is a molecularly heterogeneous disease, and so LncRNAs’ expression pattern in other subtypes needs to be identified. To the best of our knowledge, this study is the first comprehensive description of differentially expressed LncRNAs in HER-2-enriched subtype breast cancer. We sequenced seven pairs of HER-2-enriched subtype breast cancer and normal tissue and found significantly differently expressed genes. Furthermore, functions of these genes were analyzed by bioinformatics.

Pathway act network showed that MAPK signaling pathway, PI3K-Akt signaling pathway, metabolic pathways, and cell cycle and regulation of actin cytoskeleton were highly related with HER-2-enriched subtype breast cancer. It has been reported that MAPK signaling pathway could promote proliferation and angiogenesis. HER-2-targeted drugs, trastuzumab and lapatinib, may exert their action through MAPK signaling pathway.23 Agents targeting the PI3K-Akt pathway could restore sensitivity to standard breast cancer therapies, including endocrine, HER-2-targeted agents, and chemotherapy.24 Guerram et al25 found that targeting inhibited apoptotic and metabolic signaling pathways could lead to cancer cell death, especially in HER-2-overexpressed breast cancer. It was also pointed out that suppressed actin cytoskeleton pathway inhibited cell motility in breast cancer.26 These reports are consistent with our finding that these pathways are worth being further studied for potential therapeutic value.

Coexpression network showed six LncRNAs as core genes in HER-2-enriched subtype breast cancer. Except for LINC00511 (also known as onco-LncRNA-12), there is a dearth of research on the other LncRNAs. LINC00511 was reported to be upregulated in breast cancer and led to lung cancer cell growth decrease when LINC00511 was knocked down.27 Hence, LINC00511 was defined as an oncogene. This supports our findings that these LncRNAs are crucial genes in tumorigenesis. The remaining five LncRNAs, which might hopefully become candidate therapeutic targets and new molecular biomarkers, are still worth further research.

Limitations

Our study still has several limitations. For one thing, breast cancer is a highly heterogeneous disease while our sample size of sequenced tissues is limited. Therefore, our results might not represent robust LncRNAs expression signature in HER-2-enriched subtype breast cancer. Another limitation is that the functions of these core LncRNAs have yet to be determined. Further in vitro and in vivo experiments are currently being conducted by our group to investigate the biological effects of the aforementioned LncRNAs.

Conclusion

The dysregulated LncRNAs and mRNAs expression profiles were sequenced and analyze in HER-2 enriched subtype breast cancer in this study. These results analyze the functions of LncRNAs and provide useful information for exploring candidate therapeutic targets and new molecular biomarkers for HER-2 enriched subtype breast cancer.