Abstract

Background

The pancreas is both an exocrine and endocrine gland, playing important roles in the digestive and endocrine systems. There are two kinds of cells in the pancreas: exocrine cells and endocrine cells. When exocrine cells grow out of control, they may form pancreatic exocrine tumors. About 95% of pancreatic cancers can be classified into pancreatic exocrine tumors. One kind of pancreatic exocrine tumor called pancreatic ductal adenocarcinoma (PDAC) is the most common type, making up more than 85% of all pancreatic cancers. PDAC is the fourth leading cause of cancer related death in the world with a 5-year survival rate of only 5% [1]. Surgery is by far the most effective treatment strategy for PDAC, but less than 20% of PDAC patients have resectable tumors at the time of diagnosis [2, 3], with the improving 5-year survival rate after resection to 10–25% [4, 5]. The etiology of PDAC are poorly understood. However, several factors like cigarette smoking [6], family history of pancreatic cancer [7], diabetes [8] and chronic pancreatitis [9] are contributing factors for PDAC.

Like other malignancies, the intra-tumoral heterogeneity makes PDAC not a single disease, but a group of biologically and clinically distinct diseases [10, 11]. Thus, there is a great need to identify homogeneous groups which is an essential step towards personalized treatment of PDAC. Traditional classification of PDAC has been carried out by pathologists based on histologic appearance and phenotypic traits. However, in reality, tumors with similar morphological appearance may have very distinct molecular features and clinical outcomes [12, 13]. Recent advancements in genome wide molecular profiling may change these situations by providing an opportunity to investigate the tumor heterogeneity at the whole genome level. Gene expression profiling, one of the most commonly used molecular profiling approaches, is the measurement of the expression levels of thousands of genes simultaneously. And, microarray and RNA sequencing (RNA-Seq) are the two most used techniques. Gene expression profiling have allowed researchers to classify cancers into homogeneous groups with improved diagnosis [14–16] and correlated better with survival information than traditional classification of cancers [17]. Over the last few years, increasing molecular classification studies have been conducted in PDAC which proved that it can be classified into 2 to 4 subgroups [18–24]. However, these studies used tumor samples ranging from dozens to more than few hundreds as their discovery cohort. They may not fully represent the intra-tumoral heterogeneity and limit the ability to identify rare subtypes of PDAC.

Another concern in dissecting the tumor heterogeneity is the methods used in the identification process. Given a set of gene expression profiles, clustering, a machine learning technique, can be used to group data objects of similar characteristics together into distinct clusters without prior assignment (unsupervised classification). There are three kinds of clustering strategies [25]: first, gene-based clustering, which the genes are treated as the objects, while the samples are the features. Second, sample-based clustering which the samples are the objects and genes are the features. And third, biclustering (or subspace clustering) which capture clusters formed by a subset of genes across a subset of samples. The previous two strategies apply a global model to identify clusters. That is, each sample in a subtype is determined by the activity of all the genes. Similarly, each gene in a given gene cluster is defined using all the samples when performing the clustering analysis [26]. Since subsets of genes are active or silent only under certain experimental conditions, and behave almost independently under other conditions [26], the classification results are relatively poor when using the global model [27].

Only biclustering employ a local model to identify coherent patterns in an expression matrix. Instead of clustering gene and sample separately, biclustering allows simultaneous clustering of genes and samples [26]. Thus, biclustering has become a popular technique and lots of algorithms are proposed, such as distance-based [28, 29], factorization-based [30, 31] and geometric-based biclustering [32, 33]. Most biclustering algorithms [34–38] allow bi-clusters to have partially overlap, and some objects (samples or genes) may not belong to any bi-cluster at all [39, 40]. This character of biclustering, although useful in some instances [26], is not good for interpretation. Non-negative Matrix Factorization (NMF), a dimensionality reduction and factorization-based biclustering algorithm, aims to find groups of linear combination of metagenes representing local patterns in the expression matrix. NMF has been proven useful in many cancer subtyping studies [18, 20, 23, 41, 42] due to its easy interpretation and desired performances.

In our study, we focused on using NMF to extract biclusters from gene expression data, thus to describe and characterize the heterogeneity of PDAC. We overcame the sample shortage by combining different sources of PDAC into a single and large dataset. Specifically, we collected publically available PDAC gene expression profilings from 11 microarrays and 3 RNA-Seq datasets. In total, our study involves more than 1200 PDAC patients, and 796 of them were used as the discovery cohort. This is the largest cohort of PDAC gene expression profiles investigated so far, which greatly increased the statistical power and provided more robust results. We identified six molecular and clinical distinct subtypes, and provided a deep learning-based classification system for PDAC. Compared with previous studies [18–24], our study has several advantages. First, we included more PDAC cases to increase statistical reliability. Second, we selected genes as subtype-specific biomarkers directly from biclusters. Third, we identified six subtypes to provide and describe a more complete picture of the PDAC heterogeneity. Last but not least, we used deep learning to build a classification system for PDAC, which can be used to classify new patients. The classification model will be publicly available upon request.

Methods

Generation of the PDAC classifier and classification

A classifier was built on the most representative samples and most predictive genes for each cluster. Silhouette width [47] was computed to identify the most representative samples using the R package Cluster. Subtype specific genes were determined using the extractFeatures function in the NMF package [46], with the largest row feature scores. Then, SAM (Significance Analysis of Microarrays) [48] analysis was performed to filter out unstable genes between clusters. Figure 1 summarized the classifier building process.

Open in a separate window

Fig. 1

The flowchart of the classifier building process. a Data processing step. Fourteen datasets were collected and separated into training and validation datasets. Four hundred eleven most variable genes were then selected based on the median absolute deviation (MAD >0.5), and were kept for the clustering process. b NMF clustering step. Six-cluster resulted the maximum cophenetic correlation coefficient was chosen as the optimal number of clusters. Then, NMF clustering were performed of 200 times with optimal number of clusters to obtain the consensus matrix. c Classifier building steps. A classifier was built on the most representative samples and most predictive genes for each cluster

We trained a deep learning model as the PDAC classifier using the H2O R package [49]. We split the training dataset into three parts when building the model: 60% for training, 20% for validation and the remaining 20% for testing. The parameters we used were as follows: TanhWithDropout activation, bernoulli distribution, and two hidden layers with 500 neurons each. The other parameters were set as default. The classification performance of the classifier was verified on the training and validation datasets.

Results

NMF identifies six subtypes in PDAC

We applied NMF to the merged training dataset (796 PDAC patients), and obtained 2 to 6 well-defined clusters (Additional file 3: Figure S2). Cophenetic correlation coefficients were calculated to determine the optimal number of clusters, and a peak was found at k=6 (Fig. 2a). The consensus matrix heatmap contains sharp and crisp boundaries, which implies stable and robust clustering for the samples (Fig. ​(Fig.2b).2b). Silhouette width analysis was subsequently performed to select the most representative samples within each cluster (Fig. ​(Fig.2c).2c). The average silhouette width was 0.55 (range, from 0.41 to 0.64), indicating the robustness of the classification. A total number of 781 samples (~98%) with positive silhouette width were retained to build the classifier.

Open in a separate window

Fig. 2

Classification of PDAC into 6 subtypes. a Unsupervised classification of PDAC using NMF. A peak cophenetic correlation was observed for k=6 classes. b Consensus matrix for k=6 is shown. c Silhouette information for k=6 classes. d Patient distribution in the training dataset (n=796). e Patient distribution in the merged validation dataset (n=472)

Next, 160 metagenes identified by NMF were selected as features (Table 2), together with 781 sample’s Z-score normalized data to build a deep learning classifier of PDAC. We used the H2O package to split the merged training dataset into three parts: internal training set (470 PDAC, 60%), internal validation set (152 PDAC, 20%) and internal test set (159 PDAC, 20%). The internal training set was used for building the model, the internal validation set for early stopping, and internal test set for testing the classification error. The classification errors on the internal training set and internal test set were 0.8 and 13%, respectively (Additional file 4: Table S2). The classifier can be used to classify all the 796 PDAC patients in the training dataset into six subtypes: L1 (174 patients, 21.9%), L2 (176 patients, 22.1%), L3 (93 patients, 11.7%), L4 (113 patients, 14.2%), L5 (56 patients, 7.0%) and L6 (184 patients, 23.1%) (Fig. ​(Fig.2d).2d). We also did the classification with the combined validation dataset. Patients in this dataset can also be classified into six subtypes with a similar proportions of patients being distributed among subtypes (Fig. ​(Fig.2e).2e). In addition, we found that there were 65 overlapped samples between our training and combined validation dataset. More specifically, 65 samples were overlapped between ICGCarray set (178 PDAC, microarray platform) and Bailey set (71 PDAC, RNA-Seq platform). We extracted the 65 predicted sample labels from these two cohorts and compared the similarities between them. Result shows that the two lists were similar, except that there were 17 samples with inconsistent classification results, which may be jointly caused by platform differences and the classification error of the classifier.

Table 2

Subtype specific gene lists

Tumor-specific subtypes			Stroma-specific subtypes
L1	L2	L6	L3	L4	L5
AGR2	ADM	AKAP7	ACTG2	ADAMTS1	ABCC8
ALDOB	ANGPTL4	ALB	CDH11	C1orf162	ADAMTSL2
ANXA13	C19orf33	ANPEP	COL10A1	CCL2, 19, 21	C7
AQP5	CCNB2	AQP8	COL12A1	CCR7	CHGA
ARL14	CDH3	CEL	COL8A1	CD2, CD3D, CD6, CD8A, CD36, CD48, CD52, CD69, CD79B, CD163, CD247	COLEC11
C4BPB	CDKN2A	CLPS	COLEC12	CFD	CPE
CA2	COL7A1	CPA1	GAS1	CILP	F2RL2
CDCA7	CRABP2	CPB1	GREM1	CXCL9, 10, 12	FRZB
CTSE	CST6	CTRC	LRRC17	CXCR4	G6PC2
CYP3A5	DCBLD2	CTRL	MFAP5	EVI2B	IAPP
DMBT1	DHRS9	GATM	MYH11	FAM107A	NPTX2
DPCR1	DKK1	KLK1	PDGFRL	FCN1	PAX6
F5	ENO2	LEFTY1	RGS16	FOSB	PTGDS
FAM3D	IFI44L	LGALS2	SCUBE2	FPR1	QPCT
GPX2	IFIT1	MT1G	SFRP2	FYB	RAB26
LGALS4	IGF2BP3	PLA2G1B		GIMAP7	SCG5
MMP1	IRX3	PPY		GZMA	STMN2
NPC1L1	ISG15	PRSS3		GZMB	THBS4
PGC	KRT7	REG1A		HBB	ZBTB16
PIGR	LAMA3	SERPINA5		HLA-DQA1
ST6GALNAC1	LAMB3	SLC30A2		IL1B
TFF1	LAMC2	SLC3A1		IL33
TFF2	MYEOV	TMED6		IL6
TFF3	PHACTR3			IL7R
VILL	PSCA			LTB
VNN1	PTGS2			S100A8
VSIG2	S100A4			SCARA5
	SFN			SFRP1
	SLC2A1			SLIT3
	SPRR3			SPOCK2
	UBE2C			SRGN

Open in a separate window

Functional annotation of PDAC subtypes

There are distinct gene expression patterns between subtypes as observed in the heatmaps from both training and merged validation datasets (Fig. 3a-​-b).b). In the heatmap, columns correspond to PDAC patients, and rows to 160 genes. Gene expression matrices were median centered and expression values were represented by different colors, red means higher expression values, and green, lower ones. We have found that carbohydrate metabolism genes such as ALDOB, CA2, NPC1L1 and PGC are highly expressed in L1. Cell proliferation and epithelium-associated genes, such as CCNB2, CDKN2A, SFN, UBE2C, SPRR3, DHRS9 and CRABP2 are enriched in L2 subtypes. GREM1, MFAP5, COL12A1, COL10A1, COL8A1 and other collagen or ECM-related genes are upregulated in L3. Immune related genes such as CCL, CCR7 and CD gene families are enriched in L4 subtype. Neuroendocrine-associated genes such as PAX6, IAPP, G6PC2, ABCC8 and ZBTB16 are highly expressed in L5. And lastly, CLPS, PLA2G1B, CEL, ALB, CPA1, CPB1, CTRL, SLC3A1, PRSS3 and ANPEP, which are involved in lipid and protein metabolism, are enriched in the L6 subtype (Table ​(Table22 and Additional file 5: Figure S3).

Open in a separate window

Fig. 3

Functional annotation of PDAC subtypes. a Heatmap showing six subtypes of PDAC in training dataset using the 160 subtype specific genes, which reveals distinct gene expression patterns between subtypes. b Heatmap also showing six subtypes of PDAC in merged validation dataset using the 160 subtype specific genes, with similar gene expression patterns (subtype specific genes are highly expressed in the corresponded subtype) as observed in the training dataset. c GSEA analysis reveals distinct enriched gene sets between subtypes. In the heatmap, rows are defined by the selected 60 gene sets, and columns by consensus scores for each subtype. Subtype enriched gene sets are highlighted by different color, L1 (light red), L2 (light brown), L3 (light blue), L4 (light orange), L5 (light purple) and L6 (light green). d-f Heatmaps showing six subtypes of PDAC in three independent validation datasets (GSE71729, Bailey and Collisson) using the 160 subtype specific genes, with similar gene expression patterns (subtype specific genes are highly expressed in the corresponded subtype) as observed in the training dataset

To identify gene sets enriched in each subtype, we then performed GSEA analysis. GSEA is a widely used method to interpret expression data at the level of gene sets, or groups of genes that share a common biological function, or regulation [52]. We subsequently selected altogether 60 most representative gene sets for L1-L6 to build a pathway heatmap, which reveals distinct gene sets enriched in each subtype (Fig. ​(Fig.3c).3c). Based on the biological functions of the selected gene sets, we further grouped the six-subtype into tumor-specific and stroma-specific subtypes. Tumor-specific subtypes include L1, L2 and L6, which are associated with cell proliferation and metabolism-related gene sets. Specifically, L1 has enriched carbohydrate metabolism-related gene sets. L2 is enriched for cell proliferation and epithelium-associated gene sets. And L6 is enriched for lipid and protein metabolism-related gene sets. Stroma-specific subtypes include L3, L4 and L5, which may contain high nonepithelial contents such as collagen, immune and islet cells, respectively. For instance, L3 is enriched for collagen and ECM related gene sets. L4 is enriched for various immune related gene sets. And L5 is enriched for neurotransmitter and insulin secretion related gene sets. Significantly enriched gene sets for each subtype were displayed in Additional file 6: Table S3.

Clinical characterization of PDAC subtypes

About 348 patients (~43.7%) in the training dataset have clinical information. Their subtype labels and associated overall survival information were used to perform survival analysis and clinical characterizations. Kaplan-Meier analysis indicated that L2 has the worst clinical outcomes compared with other five subtypes (Fig. 4a). During the first 24 months after diagnosis, approximately 75% patients in L2 and L3, respectively, were censored. And the death rate in L2 was larger than that in L3, as observed in a steeper slope in the survival curves (Fig. ​(Fig.4a).4a). Although there were no significant survival differences in L1, L3, L4 and L6 during the first 20 months after diagnosis, the survival differences were observed after 20 months, and the death rate of L3 and L6 rapidly increased compared with L1 and L4. L5 always has good clinical outcomes compared with the other 5 subtypes. We also observed a similar overall survival differences between subtypes in the merged validation dataset (Fig. ​(Fig.4b).4b). Lastly, we did the Fisher’s exact test to investigate if the subtype memberships have any associations with other clinical factors, such as age, gender, race, tumor stage and grade. Results shows that only tumor grade have certain associations with subtypes (p-value <0.01). For example, more than 97% patients in L2 and more than 95% patients in L3 have moderately or poorly differentiated tumor cells, whereas about 32% patients in L5 have well differentiated tumor cells (Additional file 2: Table S1). This analysis demonstrates that other clinical factors (such as age, gender, race and tumor stage) cannot predict overall survival, and supports the use of subtypes as a new and reliable prognostic factor in PDAC.

Open in a separate window

Fig. 4

Clinical characterization of PDAC subtypes. a Kaplan-Meier survival curve comparing survival of L1 (red), L2 (brown), L3 (blue), L4 (orange), L5 (purple) and L6 (green) patients in the training dataset. Survival difference was tested using the log-rank test. P-values of less than 0.05 were considered statistically significant. b Kaplan-Meier survival curve in the merged validation dataset. c-e Kaplan-Meier survival curves in GSE71729, Bailey and Collisson datasets

Cross comparison of the identified subtypes with published studies

To compare our classification system with three previously published results [18, 20, 23], we then used our PDAC classifier to classify these three cohorts, separately. Gene expression heatmaps (Fig. ​(Fig.3d3d-​-f)f) and survival curves (Fig. ​(Fig.4c4c-​-e)e) show much similar patterns between validation datasets and the training dataset, which indicate the existence of six subtypes in other cohorts as well. Although some inconsistent results exist, such as the log rank p-value was not significant in GSE71729 dataset, and the survival curves in all three datasets were not followed the exact patterns as observed in the training dataset. We believe such inconsistency were caused by the smaller sample size in the validation datasets (145 PDAC in GSE71729, 71 PDAC in Bailey and 27 PDAC in Collisson set), as compared with a larger cohort size in the training dataset (796 PDAC). The corresponding sample labels in these three datasets were downloaded from the published papers, contingency tables were subsequently built and visualized by heatmaps (Fig. 5a-​-d).d). L1 and L6 were much similar to the GSE71729’s classical subtype. L2 was close to the GSE71729’s basal subtype. L4 was resemble to the GSE71729’s normal subtype. L6, L1 and L2 were similar to the GSE71729’s activated subtype. In the Bailey dataset, L6 was similar to the ADEX subtype. L4 and L1 were close to the immunogenic subtype. L2 was resemble to the squamous subtype, and L3 was similar to the pancreatic progenitor subtype. Lastly, L1 and L3 were similar to the Collison’s classical subtype. L6 was close to the Collison’s exocrine-like subtype. L2 was related to the Collison’s quasi-mesenchymal subtype. All these similarities corresponded well with the molecular and clinical characteristics of the six subtypes identified in our study, which confirmed the correctness of the characteristics we found on these six subtypes.

Open in a separate window

Fig. 5

Cross comparison of identified subtypes with published results. We used our classifier and three published classifiers to classify PDAC patients, respectively, which produced two-dimensional matrices with rows correspond to our classification results and columns correspond to other classification results. a Contingency heatmap of GSE71729 dataset. Numbers in the heatmap represent patient numbers. Row labels: our classifier’s results, and column labels: GSE71729 stroma classifier’s results. b Contingency heatmap of GSE71729 dataset. Numbers in the heatmap represent patient numbers. Row labels: our classifier’s results, and column labels: GSE71729 tumor classifier’s results. c Contingency heatmap of Bailey dataset. Numbers in the heatmap represent patient numbers. Row labels: our classifier’s results, and column labels: Bailey classifier’s results. d Contingency heatmap of Collisson dataset. Numbers in the heatmap represent patient numbers. Row labels: our classifier’s results, and column labels: Collisson classifier’s classification results

Discussion

Heterogeneity makes a cancer not just a single disease and this poses a significant challenge to the treatment of cancer patients. With the advent of genome-wide molecular profiling of cancers, especially the advancements in gene expression profiling technologies, researchers can depict genetic changes to better understand the heterogeneity of cancers. Compared with traditional classification of cancers, gene expression based classification can be used to classify cancers into subgroups with distinct molecular characteristics and clinical implications. Gene expression based classification of cancer was first proposed by Golub et al. [12]. The expression pattern of the 50 most informative genes was measured, and self-organizing maps (SOMs) clustering method was applied [53] to classify 38 leukemia patients into two prognostic groups without previous knowledge of these classes. This demonstrated the fidelity of cancer classification based solely on gene expression patterns [12]. In our study, we applied NMF to perform gene expression based classification of PDAC. We identified six molecular and clinical distinct subtypes, which not only proved that PDAC is a highly heterogeneous disease, but also demonstrated that gene expression based classification of cancer is molecular and clinical significant.

The identification of cancer subtypes can be difficult due to the lack of tumor samples available for study. The majority of PDAC patients (~80%) were first diagnosed with advanced tumor stages and were not suitable for resection. Some studies have overcome this problem by combining different sources of samples into their studies to increase the sample size [18, 54, 55]. Concatenating different datasets into a single dataset can be both significant and challenging. On one hand, integrating samples from various of independent studies can increase the statistical power and robustness. On the other hand, there exist batch effects or called non-biological differences between these datasets. Luckily, methods like Empirical Bayes (EB) [56], Surrogate Variable Analysis (SVA) [57] or Distance Weighted Discrimination (DWD) [58] can be used to remove such batch effects. For example, TCGA’s glioblastoma (GBM) subtyping study [59] integrated gene expression data from 200 GBM assayed on three platforms (Affymetrix HuEx array, Affymetrix U133A array and Agilent 244 K array) into a single dataset. Factor analysis and consensus hierarchical clustering [60] were subsequently performed for feature selection and cluster identification, respectively. The above work also used an independent dataset which contains 260 GBMs from four previously published datasets as validation dataset, and subtypes were predicted using 840 gene expression profiles and ClaNC (a nearest centroid-based classifier) [61]. In a recent publication of diffuse glioma subtyping study from TCGA [62], the authors used ComBat batch effect removal method [45] to combine multi-platform and multi-tumor mRNA expression data.

Using different patient cohorts, gene expression platforms and clustering methods can produce totally different classification results. For example, epithelial ovarian cancer (EOC) has been classified into 4 to 6 subtypes [63–65], colorectal cancer (CRC) 3 to 6 subtypes [66], and PDAC 2 to 4 subtypes were identified by different research groups [18–23]. Thus, integrating multiple patient cohorts to reduce the racial/ethnic and platforms differences, together with a unified clustering method for classification is necessary and important. In our study, in order to build a generalizable classification model for PDAC, we combined multiple PDAC gene expression datasets, and adjusted the introduced batch effect using ComBat. Our study involves more than 1200 PDAC patients, therefore, the statistical power was significantly increased. We have several advantages compared with previous studies [18–23], such as we identified novel subtypes in PDAC, we used NMF biclustering method to extract features which are more subtype specific, and finally we built a deep learning-based classification system for PDAC which can be used to classify new patients.

The expression profiling of the 160 genes identified from our study can stratify PDAC patients into six subtypes. And each subtype is characterized by the expression of a subset of genes which sharing similar biological functions, respectively. For example, L1 and L6 subtypes have enriched with metabolism-related genes; L2 and L3 have enriched with epithelium-associated and ECM-related genes, respectively; immune response genes in L4, and neuroendocrine related genes in L5. These specific expression profiles can be used to predict the clinical outcomes for each subtype, such as epithelium and cell proliferation gene profiles in L2 are related with poor prognosis; metabolism and ECM profiles in L1, L6 and L3 are associated with intermediate survival; and immune and neuroendocrine-associated profiles in L4 and L5 are correlated with relatively good clinical outcomes. The subtyping results from our study can also be interpreted at the cellular level. Low tumor cellularity and the presence of abundant stroma intermixed with normal cells are the common features of PDAC [20, 24]. Although microdissection can be used to enrich tumor cells, non-tumor components still account for a significant proportion in PDAC tissue biopsies. For example, stroma comprises on average 48% of Moffitt et al. [20] primary tumor samples with a standard deviation of 30%; and in the TCGA’s samples [24], the tumor purity ranged from 0 to 53% (median 18%). Current tissue-level expression profiling technologies process thousands of tumor and non-tumor cells at the same time, so differences or heterogeneity between patients may also result from changes in the proportions of cell types in samples. If so, then what machine learning models learned from bulk data is the cell-type proportions among samples, which can be benefit from a large group of patient’s data. Perhaps cell-type proportions are really informative, which have important implications in the treatment strategies for cancer patients. In our study, epithelial cells concentration in tumor-specific subtypes (L1, L2 and L6) may greater than the stroma-specific subtypes (L3, L4 and L5), which may suggest that L3, L4 and L5 should be treated differently from L1, L2 and L6. For instance, malignant epithelial cells in L2 may account for the largest proportion, thus more intensive treatments should be considered for L2 patients. L1 and L6, two metabolism-related subtypes may be treated by some metabolic drugs [67]. Furthermore, collagen-targeted therapies for L2, immunotherapies for L4, and endocrine cell therapies for L5.

Conclusions

In summary, we have identified six biologically informative subtypes of PDAC, which corresponding well with their molecular features and clinical outcomes. The 160 subtype specific biomarkers and the deep learning model have the potential to drive personalized therapies [68] and risk prediction [69] for the PDAC patients.

Additional files

Acknowledgements

Abbreviations

Authors’ contributions

LZ conceived and designed the study; carried out data processing, computational analysis and wrote the manuscript. HYZ helped check the results. HY did language correction, and provided the funding support for the study. All authors read and approved the final manuscript.

Notes

Footnotes

Contributor Information

Lan Zhao, Phone: +852-3442-9902, Email: kh.ude.uytic.ym@c-5oahznal, Email: moc.liamg@10104102oahznal.

Hongya Zhao, Email: kh.ude.uytic@oahzyh.

Hong Yan, Email: [email protected].

References