Impact Factor 3.258 | CiteScore 2.7
More on impact ›


Front. Genet., 30 October 2019 |

Optimization Techniques to Deeply Mine the Transcriptomic Profile of the Sub-Genomes in Hybrid Fish Lineage

Zhong Wan1*, Jiayi Tang1, Li Ren2, Yamei Xiao2 and Shaojun Liu2*
  • 1School of Mathematics and Statistics, Central South University, Changsha, China
  • 2State Key Laboratory of Developmental Biology of Freshwater Fish, Hunan Normal University, Changsha, China

It has been shown that reciprocal cross allodiploid lineage with sub-genomes derived from the cross of Megalobrama amblycephala (BSB) × Culter alburnus (TC) generates the variations in phenotypes and genotypes, but it is still a challenge to deeply mine biological information in the transcriptomic profile of this lineage owing to its genomic complexity and lack of efficient data mining methods. In this paper, we establish an optimization model by non-negative matrix factorization approach for deeply mining the transcriptomic profile of the sub-genomes in hybrid fish lineage. A new so-called spectral conjugate gradient algorithm is developed to solve a sequence of large-scale subproblems such that the original complicated model can be efficiently solved. It is shown that the proposed method can provide a satisfactory result of taxonomy for the hybrid fish lineage such that their genetic characteristics are revealed, even for the samples with larger detection errors. Particularly, highly expressed shared genes are found for each class of the fish. The hybrid progeny of TC and BSB displays significant hybrid characteristics. The third generation of TC-BSB hybrid progeny (BTF3 and TBF3) shows larger trait separation.


Taxonomy aims to define and name groups of biological organisms on the basis of their shared similarity in morphological structure and physiological functions (Tautz et al., 2002). It plays an important role in understanding the relationship and evolution between different groups (Tautz et al., 2003). From classical morphology to new achievements in modern molecular biology, taxonomy also involves the comprehensive application of biological multidisciplinary, which can be used as a basis for classification, such as chromosome-based cell taxonomy (or chromosomal taxonomy), serum taxonomy based on serum reaction, chemical composition-based chemical taxonomy, and DNA taxonomy, with the sequence analysis of a uniform target gene (Stoeckle, 2003).

In the past two decades, with an increasing number of genome-wide sequencing and fine mapping, extensive data on transcriptomics, proteomics and metabolomics are available in the literature (Liu et al., 2016; Ren et al., 2016; Ren et al., 2017a; Ren et al., 2017b; Floriou-Servou et al., 2018; Li et al., 2018; Wang L. et al., 2018; Wang M. et al., 2018; Wang N. et al., 2018; Ye et al., 2018; Chen et al., 2019; Liu et al., 2019; Ning et al., 2019). To mine more and more biological information from these data, many computational models have been established to classify different species or examine their genetic relationships (Yang et al., 2015;Tan et al., 2019). For example, in (Wang L. et al., 2018; Wang M. et al., 2018; Wang N. et al., 2018; Yu et al., 2015; Wang et al., 2017; Hu et al., 2012), some statistical methods and statistical softwares have been used for biological classification by analyzing the data of protein sequences. However, to our best knowledge, there exists no research result on classification of distant multi-generation hybrid fishes in virtue of transcriptome data and optimization techniques.

Distant hybridization is a hybrid between two different species (Lou and Li, 2006). For this interspecific hybridization, it may be a hybrid of different species of the same genus, or between different genus, between different subfamilies, between different families, and even between different individuals Zhang et al., 2014). Since distant hybridization can transfer a set of genomes from one species to another, it can effectively change the genotype and phenotype of hybrid progeny (Liu et al., 2001). In terms of genotype, distant hybridization can lead to changes in the genomic level and sub-genome levels of the offspring, and the formation of these different hybrid progeny often depends on the genetic relationship of the parent. In terms of phenotype, the distant hybridization can integrate the genetic characteristics of the parents, which may make hybrid progeny show heterosis in aspects of shape, growth rate, survival rate and disease resistance (Hu et al., 2012). It has been shown that the distant hybridization occurs widely in fishes and has become an effective tool to integrate existing natural species and quickly cultivate more excellent traits in fisheries development. For more details, readers are referred to recently published article (Qin et al., 2014; Hu et al., 2019) and the references therein.

Different from protein (DNA) sequences, the transcriptome of a cell or a tissue is the collection of RNAs transcribed in it, and is often dynamic and a good representative of the cellular state (Carnes et al., 2018). Ease of genome-wide profiling using sequencing technologies further makes the transcriptome analysis an important research tool of bioinformatics, where the information content of an organism is recorded in the DNA of its genome and expressed through transcription (Kaletsky et al., 2018). Therefore, full-length transcriptome analysis of distant multi-generation hybrid fishes seems to be a more useful tools to provide a more profound explanation for the biological performance of distant multi-generation hybrid fishes. However, on the one hand, cultivating new generation of hybrid fishes often needs more than one and a half years, hence collection of the relevant experimental data is difficult, such that only the small-size sample inference can be made (Rogoza, 2019). On the other hand, owing to a lack of effective classic statistical methods to analyze the small-size and full-length transcriptome sample data, genomic research on similarity of this species and its descendants based on optimization models is unavailable in the literature. Actually, since the full-length transcriptome data is associated with expressed levels of ten thousands genes, classification of small-size sample data becomes impossible by using existing statistical methods. In this paper, combining the RNA sequencing group data of distant hybrid progeny and parental types, we intend to develop a new method for the genetic regulation of the whole transcriptome to statistically analyze the distant hybrid progeny and its excellent germplasm selection.

Basically, our new research method originates from optimization techniques, called a nonnegative matrix factorization method (NMF). By this method, we attempt to approximately factorize the small-size and full-length transcriptome sample data of the distant multi-generation hybrid fishes such that their classification and the gene-expression characteristic of each class can be revealed. As a result, it is associated with solution of large-scale optimization problems with nonnegativity constraints. Therefore, we also aim to develop an efficient algorithm for solving this large-scale optimization problem.

Clearly, one of the challenges in this research lies in making statistical inference from the small-size samples. We have collected 24 samples (liver tissues) of the distant multi-generation hybrid fishes, which constitutes three different groups corresponding to the three sampling periods. Each group consist of 20093 genes expression levels of eight different fish. Actually, the classical statistical methods, such as k-mean clustering method and the principal component analysis (PCA), are inappropriate to analyze this type of data (8 samples with 20093 features). As stated in (El-Shagi, 2017; Ristic-Djurovic et al., 2018; Rogoza, 2019), if the size of samples is small, it is difficult to believe that the classical statistical methods cangive good prediction accuracy owing to bias of small-size samples. For the small-size samples, the existing main inference methods include: the probabilistic index models (Amorim et al., 2018), the bootstrapping U-statistics method (Jiang and Kalbfleisch, 2012), the Jackknife empirical likelihood inference (Zhao et al., 2015), the SVM-based methods (Cong et al., 2016), the grey-theory-based methods (Meng et al., 2017), and the neural network (Zhu et al., 2019). However, for the small-size samples with more than ten thousand features, such as the full-length transcriptome sample data of the distant multi-generation hybrid fishes, it is desirable to study new statistical inference methods to mine their statistical information.

The NMF has been regarded as a useful tool of unsupervised machine learning to classify the small-size samples with large-scale features (Pauca et al., 2006;Wan et al., 2018). It can integrate the functions of k-mean clustering method and PCA. However, the performance of NMF depends significantly on the development of efficient algorithms to solve the generated large-scale optimization problem such that the deviation of nonnegative matrix (sample data) factorization is minimized. Especially, if we need to classify 8 full-length transcriptome data of distant multi-generation hybrid fishes, it is necessary to factorize a matrix in R20093×8. Suppose that there are r classes of fishes, then the number of design variables is 20093 × r + 8. For solving such a large-scale optimization model, it is still a challenge to develop an efficient algorithm. In this research, we intend to modify the spectral conjugate algorithm in (Deng et al., 2013) to solve the generated large-scale optimization problems. Our goal is to reveal the relationship between multi-generation hybrid fishes on the basis of their gene expression profile described by their transcriptome data.

Materials and Methods

Samples and Transcriptome Sequencing

The Megalobrama amblycephala or Bluntnose black bream (BSB, 2n = 48 ) and Culture Alburnus or Topmouth culter (TC, 2n = 48) at sexual maturity in natural waters of the Yangtse River in China were collected for creating the allodiploids BT (BSB (♀)× TC (♂) ) and TB ( TC (♀) × BSB (♂) ) F1 individuals through intergeneric reciprocal crosses of BSB and TC, respectively. Then, the allodiploid F2 − F3 (2n = 48) hybrid offspring were obtained by self-mating of F1 – F2 populations, respectively. The chimeric offspring was identified based on 45S rDNA sequencing characteristics (Xiao et al., 2016), had been used in our study.

Transcriptome Sequencing and Gene Expression Profiles

To sequence the transcriptomes of reciprocal cross hybrids and their inbred parents, total RNA was isolated and purified from the liver by a TRIzol extraction method (Rio et al., 2010). RNA concentration was measured using Nanodrop technology. Total RNA samples were treated with DNase I (Invitrogen) to remove any contaminating genomic DNA. The purified RNA was quantified using a 2100 Bioanalyzer system (Agilent, Santa Clara, CA, USA). After the isolation of 1 μg mRNA using the beads with oligo (dT) Poly (A), fragmentation buffer was added for interrupting mRNA to short fragments. The resulting short fragments were reverse transcribed and amplified to produce cDNA. An Illumina RNA-seq library was prepared according to a standard high-throughput method ephigh-throughput method (Dillies et al., 2013). The cDNA library concentration and quality were assessed by the Agilent Bioanalyzer 2100 system, after which the library was sequenced with paired-end setting using the Illumina HiSeq 2000/4000 platform. Then, the raw reads containing adapters, ploy-N and low quality were removed using in-house perl scripts. The high quality reads were used in our analysis. The transcriptome data was obtained from the NCBI database.

All Illumina reads of M. Amblycephala and C. alburnus were aligned to the M. Amblycephala and C. alburnus genome using Star (v 2.4.0) with the default parameters (Bennett et al., 2001), respectively. The other RNA-seq reads of reciprocal cross hybrids were aligned to the two reference genomes of M. Amblycephala and C. alburnus, respectively. The numbers of mapping counts in each gene were calculated with in-house perl scripts. Consequently, the two mapping results of aligning to two reference genomes were obtained in hybrid offspring, and the total expression value was normalized based on ratio of the number of mapped reads at each gene to the total number of mapped reads for the entire genome.

Data Download

The collected data of 24 samples (liver tissues) of the distant multi-generation hybrid fishes in this research have been uploaded to and can be downloaded freely such that the numerical experiments in this paper can be repeated by anyone. The last upload time is 07-20-2019(File name as 2019_7_8 Copy.xlsx).

An Optimization Model for Classifying the Hybrids Fishes

We first propose an optimization model for classifying the hybrids fishes on the basis of NMF. Mathematically, NMF is stated as follows. For a given matrix ARn × m, we need to decompose A into two nonnegative matrices W and H, i.e.


where WRn × r and HRr × m. In particular, if the matrix A in (2.1) is the full-length transcriptome data of the distant multi-generation hybrid fishes, and A = WH, then r can represent the number of classes for this classification of fishes in the case that each column of H has only a unique element 1, while the other elements are zeros. Clearly, in this ideal case, the k-th column of W stands for the gene expression level of the k-th class of fishes, and its elements show the expression levels of different genes for each class. Therefore, W in Model (2.3) is called a base matrix in view of its practical meanings, while H is called a coordinate matrix.

For real sample data, it is often difficult to obtain the above ideal result of factorization. Therefore, we relax A = WH by AWH. In this case, each column of the matrix A is approximately equal to the linear combination of all column vectors of the matrix W, and the combination coefficients are given by the corresponding column vector of the matrix H, i.e. A:,jW:,k×hk,j, where A:,j denotes the j-th column of the matrix A, W:,k stands for the k-th column of the matrix W, and hk,j represents the element of the k-th row and the j-th column in the matrix H. In other words, A=[A:,1,A:,m]Rn×m, W=[W:,1,W:,r]Rn×r, and H=[hk,j]Rr×m.

Thus, if we define a membership matrix RRr × m:


Clearly, the j-th column of R represents the membership degrees of the j-th sample being affiliated all the different classes. Therefore, for all the samples, distinct differences of all the elements in each column of R imply an approximate classification result. By definition, the matrix R shows the result of classification in term of membership degrees, while each column of the matrix H exactly stands for the coordinate of each sample in the r-dimensional space linearly expanded by the r columns of W. In the case that all the r elements in each row of W have the same orders of magnitude, the classification results by H or R are same.

Unfortunately, it is very difficult to solve Problem (2.1) when n is very large, let alone the requirement of finding the unknown optimal number of classes r. To solve Problem (2.1), we first transform (2.1) into the following optimization model:


where ‖·‖F is the Frobenius norm. It has been shown that (2.3) is non-convex and NP-hard (Vavasis, 2009). Then, similar to the technique of alternating non-negative least squares (ANLS) in (Chu et al., 2004), we solve (2.3) by finding the optimal solutions of the following two convex sub-problems:


It is noted that the above model of NMF was first proposed in (Paatero and Tapper, 1994). Summarily, there are two types of algorithms to solve Model (2.3) (Lin, 2007): the multiplicative update (MU) method (Cai et al., 2010; Shang et al., 2012; Huang et al., 2018; Deng et al., 2019) and the technique of alternating non-negative least squares (ANLS) (Chu et al., 2004). For the ANLS, a main focus is on development of efficient algorithms to solve the subproblems (2.4) and (2.5). For example, the projected gradient (PG) method (Lin, 2007), the projected Newton method (Gong and Zhang, 2012), and the projected quasi-Newton method (Zdunek and Cichocki, 2006) have been reported to be efficient for solving the large-scale optimization model (2.3), although no one method has overwhelming advantage compared with the others.

Recently, Deng et al. (2013)proposed an efficient algorithm to solve general large-scale unconstrained optimizations, and they demonstrated that the numerical performance of this algorithm outperforms the similar ones available in the literature. In this paper, we intend to extend it into solution of the subproblems (2.4) and (2.5), which are two large-scale optimization problems with nonnegativity constraints.

Development of Algorithm

We are now in a position to present an efficient algorithm to solve the subproblems (2.4) and (2.5). Since both of them are large scale (the size of the problem is over 80000), we will extend the spectral conjugate gradient algorithm in (Deng et al., 2013) to solve the subproblems (2.4) and (2.5). Actually, in our previous research, this algorithm has been implemented to solve more than 700 large-scale benchmark test problems, and has been shown that its numerical performance outperforms the similar ones available in the literature.

In need of modifying the developed algorithm in (Deng et al., 2013) such that it can be used to solve Model (2.3), we first define the gradients of F in (2.4) and (2.5) with respect to the matrices W and H, respectively.

By direct calculation, it is easy to see that for any i and j,


Then, we denote the following two matrices the gradients of F(W, H) with respect to the matrices W and H, respectively:


For two given matrices S and T with the same size, we define their inner product by


Then, for k = 0, a search direction of F at a given initial point W(0) is


And for k ≥ 1, we define four matrices:


where H(k), W(k) and W(k − 1) are two given matrices. Similar to (Deng et al., 2013), we compute the spectral parameter and conjugate parameter by

θk={Dk1,yk1Dk1,GW(k)GW(k),sk1GW(k)2Dk1,yk1¯,if Dk1,yk1¯>ηGW(k1)2,Dk1,yk1Dk1,GW(k)GW(k),GW(k1)GW(k)2Dk1,GW(k1),otherwise,(2.11)


βk={GW(k),yk1GW(k),sk1Dk1,yk1¯,if Dk1,yk1¯>ηGW(k1)2,GW(k),yk1GW(k1)2,otherwise,(2.12)

where Dk−1 is the search direction at W(k−1), determined by

Dk={D0,if  k=0,θkGW(k)+βkDk1,if  k>0.(2.13)

The following algorithm is developed to solve the subproblem (2.4) with the given H(k).


Algorithm 1 (Modified Spectral Conjugate Gradient Algorithm)

Similarly, to solve the subproblem (2.5), we only need replace W and H by H and W in Algorithm 1, respectively. Particularly, we need to compute



Dk={D0,if  k=0,θkGH(k)+βkDk1,if  k>0,(2.15)


θk={Dk1,yk1Dk1,GH(k)GH(k),sk1GH(k)2Dk1,yk1¯,if  Dk1,yk1¯>ηGH(k1)2,Dk1,yk1Dk1,GH(k)GH(k),GH(k1)GH(k)2Dk1,GH(k1),otherwise,


βk={GH(k),yk1GH(k),sk1Dk1,yk1¯if  Dk1,yk1¯> ηGH(k1)2,GH(k),yk1GH(k1)2,otherwise,

With the above preparation, we now develop an overall algorithm to solve Model (2.3) in the end of this section.


Algorithm 2


In this section, in virtue of Model (2.3) and Algorithm 2, we present the results on classification of the distant multi-generation hybrid fishes based on their transcriptome data.

Result Of Classification

With the given transcriptome data of the distant multi-generation hybrid fishes, we easily get Model (2.3). Then, we implement Algorithm 2 to solve this model by choosing the same values of algorithmic parameters as in (Deng et al., 2013):

ε=107, δ1W=δ1H=0.4, ηW=ηH=0.001, δ2W=δ2H=0.001, ρW=ρH=0.65.

In addition, for any choice of, ρW, ρH ∊ [0.05, 0.75] we can obtain the almost same results in our numerical experiments, which indicates our algorithms are robust for classifying the fishes.

All codes of the computer procedures are written in MATLAB and run in a MATLAB R2016b, and are carried out on a PC(CPU 2.40 GHz,8G memory) with the Windows 10 operation system environment. All the codes have been uploaded to

For the sake of better understanding the inherent characteristics of the data, we take the 2nd-group samples with superscripts L2 as a training set, which were from the liver tissue of eight different fish. Since it is unclear how many classes can be identified for the fish samples before our research, we make a trial setting on the number of classes r = 2, …, 7 such that the best number of classes is found.

In Table 1, we report all the numerical results corresponding to the different class numbers.


Table 1 Coordinate matrices for the 2nd-group samples.

Table 1 shows that when r = 6, all the samples are clearly classified owing to existence of greater deviation of elements in the same column of H. In contrast, when r is equal to the other values, there are at least one sample that can not be clearly classified.

As r = 6, Table 1 indicates that the eight fishes can be categorized into 6 classes: BSBL2, TCL2, TBF2L2 and BTF3L2 belong to different four classes, respectively. BTF1L2 and BTF2L2 consist in another class. TBF1L2 and TBF3L2 are the same class.

For the sake of better understanding the above classification result, we use r = 6 as the number of classes to calculate the membership matrix R defined by (2.2). The numerical results are listed in Table 1, while Figure 1 more intuitively describe the biological similarity for the fish of each class.


Figure 1 Membership degrees of the 2nd-group samples.

Table 2 and Figure 1 further indicate that by membership matrices, the same classification result is obtained as that by coordinate matrices: {BSBL2}; {TCL2}; {TBF2L2}; {BTF3L2}; {BTF1L2,BTF2L2}; {TBF1L2,TBF3L2}. Particularly, either by H or by R, BSBL2 and TCL2 always belong to two different classes, while their hybrids are divided into different classes from the parents’ ones. In Figure 1, Classes from 1 to 6 are described by the colors of yellow, blue, green, purple, gray and red, respectively. It follows from Figure 1 that larger proportion of the green color in BTF1L2 and BTF2L2 (that of the yellow color in TBF1L2 and TBF3L2) demonstrate that there exists greater degree of biological similarity between BTF1L2 and BTF2L2 (between TBF1L2 and TBF3L2).


Table 2 Membership matrix R of the 2nd-group samples.

To further test robustness of the above trained results, given r = 6, we choose the 1st-group and the 3rd-group samples (with superscripts L1 and L3, respectively) as two test sets to see whether the results are the same or not.

In Table 3 and Figure 2, we report the numerical results. The used colors in Figure 2 only be used to show the similarity of fishes within the same figure. In other words, the same color has no any relation in different figures.


Table 3 Results for the 1st-/3rd-group samples.


Figure 2 The membership degrees of 1st-/3rd-group samples.

From Table 3 and Figure 2, it is clear that 6 out of 8 samples in the 1st-group or the 3rd-group are correctly classified, compared with the trained result from the samples of the 2nd-group. The accuracy rate reaches 75%. In Table A3, we show that the elements in each row of the matrix W have different orders of magnitude for the 1st-group samples, which can explain inconsistence of the classification results by H and R for the 4 samples: BSBL1,BTF1L1,BTF2L1 and BTF3L1.

To further validate the proposed model and algorithms in this paper, we use them to classify more test samples generated by mixing the training set and the test sets.

We first mix the training set and the 3rd-group test set. The obtained results are listed in Table 4. Table 4 demonstrates that compared with the trained result, 13 out of 16 samples are correctly classified by both of the membership and coordinate matrices, which includes all the samples in the 2nd-group and the 5 samples in the 3rd-group: BSBL3,BTF2L3,BTF3L3,TBF2L3 and TBF3L3. The accuracy rate is as high as 81.25%. Additionally, for the 5 species of fish (BSB, BTF2, BTF3, TBF2 and TBF3), the replicated samples of each fish are correctly classified into the same class in our test experiments, which also validates the proposed model and algorithms in this paper.


Table 4 Results for the mixed samples of the 2nd/3rd group.

Next, we compute the classification result of all 24 samples (8 samples in the training set, 16 samples in the two test sets). The results are given in Table 5. From Table 5, we know that 17 out of 24 samples are correctly classified by the membership matrix or the coordinate matrix, which excludes BSBL1, BTF1L1, BTF2L1, BTF1L3, TBF1L1, TBF1L3 and TCL3. The accuracy rate achieves 70.83%, compared with the trained results. In this test, for the 4 species of fish (BSB, BTF3, TBF2 and TBF3), the replicated samples of each fish are correctly classified into the same class.


Table 5 Results for the mixed samples of all three groups.

In summary, by all of the above test experiments, the average accuracy rate is 75.52% even if there exists larger detection error of the input initial sample data (see our subsequent correlation analysis). These tests further verifies that the proposed model and algorithm in this paper can be used to efficiently classify the distant multi-generation hybrid fishes based on their transcriptomic profile.

Correlation Analysis

To find out the reasons why the replicated samples are incorrectly classified such that the accuracy rate may be reduced, we calculate the correlation matrix of the sample data to reveal possible detection errors of the input initial data. In Figure 3, the correlation coefficient matrix of the 24 samples is concisely plotted.


Figure 3 Correlation of the input 24 sample data.

From Figure 3, it is easy to see that the sample of BSBL1 is only weakly correlated with the two replicated samples BSBL2 and BSBL3. Their correlation degree is even less than that between the samples of different fish BSBL1 and TBF2. It can explain why BSBL1 can not be clearly classified into the same class of BSBL2 and BSBL3 (revisiting the results in Table 5). Conversely, Figure 3 shows that in the 1st-group, the sample BSBL1 has greater correlation with the other 3 samples: TBF2L1, TBF3L1 and TCL1, which answers why the class of BSBL1 can not be clearly identified in Table 3.

From Figure 3, we can also find out similar reasons for the unsatisfactory classification of BTF1L1, BTF1L2, BTF1L3 in Tables 4 and 5. Actually, (1) owing to lower correlation among BTF1L1, BTF1L2 and BTF1L3, they can not be classified into the same class even if they are the three replicated samples. (2) In the 3rd group, the class of BTF1L3 can not be clearly identified in Table 3 since its sample is more correlated with the other 5 samples: BSBL3, BTF2L3,TBF1L3,TBF2L3 and TCL3.

Similarly, because the sample of BTF2L1 is only little correlated with the two replicated samples BTF2L2 and BTF2L3, BTF2L1 can not be classified into the same class of BTF2L2 and BTF2L3 in Table 5.

For the same reason of weaker correlation, in Tables 4 and 5, the three replicated samples of TBF1 (TC) are also classified into the different classes. It is believed that if the detection errors of samples can be controlled to be small enough, the proposed model and algorithms in this paper can provide a more satisfactory result of classification. Actually, for the three species of fish: TBF2, TBF3 and BTF3, their three replicated samples can always classified into the respective same class (see Tables 4 and 5), which may be related with higher correlation between them as shown in Figure 3.

Genes Of High Expression

In the end of this section, based on our classification result from the 2nd-group samples, we answer what are the differently expressed genes in all the six classes. By definition, we know that each column of the base matrix W gives the feature of gene expression for each class of fish. Since the sample of each class consists of 20093 genes, we only list a part of the highly expressed genes for each fish. When r = 6, the highly expressed genes are reported in Table A1 and Figures 1(a), 1(b) and 1(c).

From the numerical results in Table A1 and Figures 1(a), 1(b) and 1(c), it follows that there exists stronger genetic similarity between the BSB (parents) and the hybrids. Actually, the BSB (the 6th class) has 3 shared highly expressed genes with TB (the 1st class), 45 shared highly expressed genes with BTF3 (the 2nd class) and 12 shared highly expressed genes with TB (the 4th class). In contrast, the TC (the 5th class) does not have any shared highly expressed genes with their hybrids, which implies that their hybrids seem to look more like BSB, rather than TC, regardless of reciprocal hybrids.

Apart from one-by-one comparison in Table A1, we also statistically analyze the numbers of shared highly expressed genes for more than three classes of fish. The reported results in Table A2 demonstrate that BSB (6-th class) has higher hereditary conservatism than TC (5th class). Actually, by comparing the numbers of shared highly expressed genes among BSB, TC and the hybrids, it is clear that the gene expression profile of their grandchildren looks more like BSB (6st class), rather than TC (5th class).

It is also noted that in Table A2, there are no shared expressed genes between BTF1 (3rd class) and TBF1 (1st class), or between BTF2 (3rd class) and TBF2 (4th class), and there only exist 3 shared highly expressed genes between the BTF3 (2nd class) and TBF3 (1st class). It suggests that the trait separation occurs between these hybrids.

In addition, from Table A2 and Figure A1, it follows that the hybrids have larger transcript intersection than that between the hybrids and the parents, since the number of shared highly expressed genes between the hybrids (offspring) is far more than that between them and their parents. Actually, there are 277 shared highly expressed genes among TBF3, TBF1 (1st class), BTF3 (2nd class) and TBF2 (4th class). In contrast, there are only less than 45 shared highly expressed genes between the parent (BSB) and the hybrids (BTF3).


In our numerical experiments, it is found that the nonnegative factorization of the matrix A is not unique. In particular, if we choose different initial matrices W0 and H0, the base and coordinate matrices W and H may be different. However, our numerical experiments show that for Algorithms 1 and 2, different choices of W0 and H0 do not affect the final result of classification. For example, as r = 6, the result of classification always is the same for any W0 and H0, which can show robustness of our classification method.

Hybridization is considered as the rapidly driving forces that shape epigenetic modifications in plants and parts of lower vertebrate (Liu et al., 2016; Mallet, 2005). The merge of divergent genome always results in a ‘genomic and transcriptome shock’ in newborn hybrid (Ren et al., 2017b; Wu et al., 2016; Ren et al., 2016). Analysis on the expression changes after hybridization, including expression dominance and expression bias related to specific function-regulated genes, always provides us insights into the molecule mechanism of various phenotypes including heterosis (Ren et al., 2016; Zhou et al., 2015). However, the multiple regulatory mechanism and complex protein interaction network restricted our ability to investigate the underlying regulation in hybrid.

It is noted that in this research, we choose the 2nd-group samples as the training set, instead of the 1st-group or 3rd-group, and the latter is regarded as test samples to verify the trained result. One of the reasons for our doing so lies in that correlation analysis of the three-group samples indicates that each sample in the second-group is better correlated with the other replicated samples than those in the other two groups.

The proposed model and algorithms in this paper can be extended to solve more practical engineering problems from other fields. For example, if we can collect sufficient transcriptome data of patients possibly suffering from breast cancer, we can apply the proposed model and algorithms to identify the classes of patients, even development of the relevant smart aided-system of diagnosis for the sufferers.


In this paper, we have constructed a classification model for the distant multi-generation hybrid fishes based on transcriptome data, and developed an efficient algorithm, called the modified spectral conjugate gradient algorithm, for solving such a model.

In virtue of our model and algorithm, we have obtained a satisfactory classification for a given full-length transcriptome data of fish samples, and the differently expressed genes of each class have been identified. Our results are first obtained by a training set of samples, then are tested by many test samples generated by different ways.

Main results are stated as follows.

1. Even for input data with larger detection error, the average accuracy rate of classification still achieves 75.52% in all the test experiments. It suggests that our model and algorithms are promising in classifying the distant multi-generation hybrid fishes.

2. Owing to the weakest intersection of highly expressed genes between BSB and TC, they are deterministically divided into two classes. However, there exists a higher transcript intersection between them and their hybrids. These findings have further deeply mined the biological genetic characteristics of distant hybridization generated by BSB and TC, based on optimization techniques and transcriptome data.

3. Although the hybrids of TC and BSB have been divided into different classes, the hybrids display higher transcript intersection. Since the transcript intersection of the hybrids and the parents is smaller than that among the hybrids, it can be concluded that the hybrid progeny of TC and BSB has significant hybrid characteristics, which may be useful to carry out trait improvement in practice.

4. Since BTF3 and TBF3 are classified to two different classes, where there only exist 3 shared genes of high expression, it is concluded that there exists larger trait separation in the third generation of TC and BSB hybrid progeny (BTF3 and TBF3). In other words, both BTF3 and TBF3 are a good variety for the reproduction of fish.

5. Since there are no shared genes of high expression between BTF1 and TBF1, they belong to two different classes (1st and 3rd classes). It implies that the reciprocal hybrids in the first generation of TC and BSB (BTF1 and TBF1) have larger biological distinction.

Data Availability Statement

The genome assembly used in this study was downloaded from NCBI BioProject database (BioProject: under accession numbers PRJNA269572. All raw mRNA-seq data were downloaded from the NCBI Sequence Read Archive ( under accession number SRP050891.

Author Contributions

ZW conceived and designed the study, wrote and revised the paper. JT designed and implemented the algorithm to analyze the data, and wrote the paper. LR, YX and SL did all the relevant experiments, collected the data and revised the paper.


This research was supported by the National Science Foundation of China (Grant 71671190), National Key Research and Development Program of China (2018YFD0901202), National Science Foundation of China (31772902), and State Key Laboratory of Developmental Biology of Freshwater Fish (2018KF003).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


The authors would like to express their thanks to the anonymous referees for their constructive comments on the paper, which have greatly improved its presentation.


Amorim, G., Thas, O., Vermeulen, K., Vansteelandt, S., De Neve, J. (2018). Small sample inference for probabilistic index models. Comput. Stat. Data Anal. 121, 137–148. doi: 10.1016/j.csda.2017.11.005

CrossRef Full Text | Google Scholar

Bennett, B. L., Sasaki, D. T., Murray, B. W., O’Leary, E. C., Sakata, S. T., Xu, W., et al. (2001). SP600125, an anthrapyrazolone inhibitor of Jun N-terminal kinase. PNAS 98 (24), 13681–13686. doi: 10.1073/pnas.251194298

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, D., He, X., Han, J., Huang, T. S. (2010). Graph regularized nonnegative matrix factorization for data representation. IEEE Trans. Pattern Anal. Mach. Intell. 33 (8), 1548–1560. doi: 10.1109/TPAMI.2010.231

PubMed Abstract | CrossRef Full Text | Google Scholar

Carnes, M. U., Allingham, R. R., Ashley-Koch, A., Hauser, M. A. (2018). Transcriptome analysis of adult and fetal trabecular meshwork, cornea, and ciliary body tissues by RNA sequencing. Exp. Eye Res. 167, 91–99. doi: 10.1016/j.exer.2016.11.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Cheng, L., Zhang, X., Cao, J., Wu, Z., Zheng, X. (2019). Transcriptomic and proteomic effects of (-)-epigallocatechin 3-O-(3-O-methyl) gallate (EGCG3” Me) treatment on ethanol-stressed Saccharomyces cerevisiae cells. Food Res. Int. 119, 67–75. doi: 10.1016/j.foodres.2019.01.061

PubMed Abstract | CrossRef Full Text | Google Scholar

Chu, M., Diele, F., Plemmons, R., Ragni, S. (2004). Optimality, computation, and interpretation of nonnegative matrix factorizations. SIAM J. Matrix Anal. Appl. doi:

Google Scholar

Cong, L., Xu, T. X., Wang, Q. (2016). Missile competing fault prediction based on degradation data and fault data. J. Beijing Univ. Aeronaut. Astronaut. 42 (3), 522–531. doi: 10.13700/

CrossRef Full Text | Google Scholar

Deng, S., Wan, Z., Chen, X. (2013). An improved spectral conjugate gradient algorithm for nonconvex unconstrained optimization problems. J. Optim. Theor. Appl. 157 (3), 820–842. doi: 10.1007/s10957-012-0239-7

CrossRef Full Text | Google Scholar

Deng, J. L., Xu, Y. H., Wang, G., Zhu, Y. S. (2019). Identification of potential crucial genes and key pathways in breast cancer using bioinformatic analysis. Front. Genet. 10, 695. doi: 10.3389/fgene.2019.00695

PubMed Abstract | CrossRef Full Text | Google Scholar

Dillies, M. A., Rau, A., Aubert, J., Hennequet-Antier, C., Jeanmougin, M., Servant, N., et al. (2013). A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief Bioinform. 14 (6), 671–683. doi: 10.1093/bib/bbs046

PubMed Abstract | CrossRef Full Text | Google Scholar

El-Shagi, M. (2017). Dealing with small sample bias in post-crisis samples. Econ. Model. 65, 1–8. doi: 10.1016/j.econmod.2017.04.004

CrossRef Full Text | Google Scholar

Floriou-Servou, A., von Ziegler, L., Stalder, L., Sturman, O., Privitera, M., Rassi, A., et al. (2018). Distinct proteomic, transcriptomic, and epigenetic stress responses in dorsal and ventral hippocampus. Biol. Psychiatry 84 (7), 531–541. doi: 10.1016/j.biopsych.2018.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, P., Zhang, C. (2012). Efficient nonnegative matrix factorization via projected Newton method. Pattern Recogn. 45 (9), 3557–3565. doi: 10.1016/j.patcog.2012.02.037

CrossRef Full Text | Google Scholar

Hu, F., Fan, J., Qin, Q., Huo, Y., Wang, Y., Wu, C., et al. (2019). The sterility of allotriploid fish and fertility of female autotriploid fish. Front. Genet. 10. doi: 10.3389/fgene.2019.00377

CrossRef Full Text | Google Scholar

Hu, J., Liu, S., Xiao, J., Zhou, Y., You, C., He, W., et al. (2012). Characteristics of diploid and triploid hybrids derived from female Megalobrama amblycephala Yih×male Xenocypris davidi Bleeker. Aquaculture 364, 157–164. doi: 10.1016/j.aquaculture.2012.08.025

CrossRef Full Text | Google Scholar

Huang, S., Wan, Z., Zhang, J. (2018). An extended nonmonotone line search technique for large-scale unconstrained optimization. J. Comput. Appl. Math. 330, 586–604. doi: 10.1016/

CrossRef Full Text | Google Scholar

Jiang, W., Kalbfleisch, J. D. (2012). Bootstrapping U-statistics: applications in least squares and robust regression. Sankhya B. 74 (1), 56–76. doi: 10.1007/s13571-012-0043-2

CrossRef Full Text | Google Scholar

Kaletsky, R., Yao, V., Williams, A., Runnels, A. M., Tadych, A., Zhou, S., et al. (2018). Transcriptome analysis of adult Caenorhabditis elegans cells reveals tissue-specific gene and isoform expression. PLoS Genet. 14 (8), e1007559. doi: 10.1371/journal.pgen.1007559

CrossRef Full Text | Google Scholar

Li, W., Liu, J., Tan, H., Luo, L., Cui, J., Hu, J., et al. (2018). A symmetric expression patterns reveal a strong maternal effect and dosage compensation in polyploid hybrid fish. BMC Genomics 19 (1), 517. doi: 10.1186/s12864-018-4883-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Wan, Z. (2019). New adaptive Barzilar-Borwein step size and its application in solving large scale optimization problems. ANZIAM J. 61 (1), 76–98. doi: 10.1017/S1446181118000263

CrossRef Full Text | Google Scholar

Lin, C. J. (2007). Projected gradient methods for nonnegative matrix factorization. Neural Comput. 19 (10), 2756–2779. doi: 10.1162/neco.2007.19.10.2756

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S., Liu, Y., Zhou, G., Zhang, X., Luo, C., Feng, H., et al. (2001). The formation of tetraploid stocks of red crucian carp×common carp hybrids as an effect of interspecific hybridization. Aquaculture 192 (2-4), 171–186. doi: 10.1016/S0044-8486(00)00451-8

CrossRef Full Text | Google Scholar

Liu, S., Luo, J., Chai, J., Ren, L., Zhou, Y., Huang, F., et al. (2016). Genomic incompatibilities in the diploid and tetraploid offspring of the goldfish×common carp cross. PNAS 113 (5), 1327–1332. doi: 10.1073/pnas.1512955113

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, F., Meng, Y., He, K., Song, F., Cheng, J., Wang, H., et al. (2019). Comparative analysis of proteomic and metabolomic profiles of different species of Paris. J. Proteomics. 200, 11–27. doi: 10.1016/j.jprot.2019.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Lou, Y. D., Li, X. Q. (2006). Distant hybridization of fish and its application in aquaculture in China. J. Fish Sci. China 13 (1), 151–158. doi: 10.1360/aps050066

CrossRef Full Text | Google Scholar

Mallet, J. (2005). Hybridization as an invasion of the genome. Trends Ecol. Evol. 20 (5), 229–237. doi: 10.1016/j.tree.2005.02.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, F., Liu, Y., Liu, L., Li, Y., Wang, F. (2017). Studies on mathematical models of wet adhesion and lifetime prediction of organic coating/steel by grey system theory. Materials 10 (7), 715. doi: 10.3390/ma10070715

CrossRef Full Text | Google Scholar

Ning, M., Wei, P., Shen, H., Wan, X., Jin, M., Li, X., et al. (2019). Proteomic and metabolomic responses in hepatopancreas of whiteleg shrimp Litopenaeus vannamei infected by microsporidian Enterocytozoon hepatopenaei. Fish Shellfish Immunol. 87, 534–545. doi: 10.1016/j.fsi.2019.01.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Paatero, P., Tapper, U. (1994). Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values. Environmetrics 5 (2), 111–126. doi: 10.1002/env.3170050203

CrossRef Full Text | Google Scholar

Pauca, V. P., Piper, J., Plemmons, R. J. (2006). Nonnegative matrix factorization for spectral data analysis. Linear Algebra Appl. 416 (1), 29–47. doi: 10.1016/j.laa.2005.06.025

CrossRef Full Text | Google Scholar

Qin, Q., Wang, Y., Wang, J., Dai, J., Xiao, J., Hu, F., et al. (2014). The autotetraploid fish derived from hybridization of Carassius auratus red var.(female) × Megalobrama amblycephala (male). Biol. Reprod. 91 (4), 93–91. doi: 10.1095/biolreprod.114.122283

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, L., Cui, J., Wang, J., Tan, H., Li, W., Tang, C., et al. (2017a). Analyzing homoeolog expression provides insights into the rediploidization event in gynogenetic hybrids of Carassius auratus red var.× Cyprinus carpio. Sci. Rep. 7 (1), 13679. doi: 10.1038/s41598-017-14084-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, L., Li, W., Tao, M., Qin, Q., Luo, J., Chai, J., et al. (2016). Homoeologue expression insights into the basis of growth heterosis at the intersection of ploidy and hybridity in Cyprinidae. Sci. Rep. 6, 27040. doi: 10.1038/srep27040

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, L., Tang, C., Li, W., Cui, J., Tan, X., Xiong, Y., et al. (2017b). Determination of dosage compensation and comparison of gene expression in a triploid hybrid fish. BMC Genomics 18 (1), 38. doi: 10.1186/s12864-016-3424-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Rio, D. C., Ares, M., Hannon, G. J., Nilsen, T. W. (2010). Purification of RNA using TRIzol (TRI reagent). Cold Spring Harb. Protoc. 2010 (6), pdb.prot5439. doi: 10.1101/pdb.prot5439

CrossRef Full Text | Google Scholar

Ristic-Djurovic, J. L., Cirkovic, S., Mladenovic, P., Romcevic, N., Trbovich, A. M. (2018). Analysis of methods commonly used in biomedicine for treatment versus control comparison of very small samples. Comput. Methods Programs Biomed. 157, 153–162. doi: 10.1016/j.cmpb.2018.01.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Rogoza, W. (2019). Method for the prediction of time series using small sets of experimental samples. Appl. Math Comput. 355, 108–122. doi: 10.1016/j.amc.2019.02.062

CrossRef Full Text | Google Scholar

Shang, F., Jiao, L. C., Wang, F. (2012). Graph dual regularization non-negative matrix factorization for co-clustering. Pattern Recogn. 45 (6), 2237–2250. doi: 10.1016/j.patcog.2011.12.015

CrossRef Full Text | Google Scholar

Stoeckle, M. (2003). Taxonomy, DNA, and the bar code of life. BioScience 53 (9), 796–797. doi: 10.1641/0006-3568(2003)053[0796:TDATBC]2.0.CO;2

CrossRef Full Text | Google Scholar

Tan, M., Long, H., Liao, B., Cao, Z., Yuan, D., Tian, G., et al. (2019). QS-net: reconstructing phylogenetic networks based on quartet and sextet. Front. Genet. 10, 607. doi: 10.3389/fgene.2019.00607

PubMed Abstract | CrossRef Full Text | Google Scholar

Tautz, D., Arctander, P., Minelli, A., Thomas, R. H., Vogler, A. P. (2002). DNA points the way ahead in taxonomy. Nature 418 (6897), 479. doi: 10.1038/418479a

CrossRef Full Text | Google Scholar

Tautz, D., Arctander, P., Minelli, A., Thomas, R. H., Vogler, A. P. (2003). A plea for DNA taxonomy. Trends Ecol. Evol. 18 (2), 70–74. doi: 10.1016/S0169-5347(02)00041-1

CrossRef Full Text | Google Scholar

Vavasis, S. A. (2009). On the complexity of nonnegative matrix factorization. Trends Ecol. Evol. 20 (3), 1364–1377. doi: 10.1137/070709967

CrossRef Full Text | Google Scholar

Wan, Z., Guo, J., Liu, J., Liu, W. (2018). A modified spectral conjugate gradient projection method for signal recovery. Signal Image Video P. 12 (8), 1455–1462. doi: 10.1007/s11760-018-1300-2

CrossRef Full Text | Google Scholar

Wang, L., Liu, P. F., Zhao, H., Zhu, G. P., Wuyun, T. N. (2018). Comparative transcriptome analysis between interspecific hybridization (Huaren apricot ♀ × almond ♂) and intraspecific hybridization (Huaren apricot) during young fruit developmentalstage. Sci. Hortic. 240, 397–404. doi: 10.1016/j.scienta.2018.06.038

CrossRef Full Text | Google Scholar

Wang, Y., Zhang, H., Lu, Y., Wang, F., Liu, L., Liu, J., et al. (2017). Comparative transcriptome analysis of zebrafish (Danio rerio) brain and spleen infected with spring viremia of carp virus (SVCV). Fish Shellfish Immunol. 69, 35–45. doi: 10.1016/j.fsi.2017.07.055

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M., Zhou, Z., Wu, J., Ji, Z., Zhang, J. (2018). Comparative transcriptome analysis reveals significant differences in gene expression between appressoria and hyphae in, colletotrichum gloeosporioides. Genetics 670, 63–69. doi: 10.1016/j.gene.2018.05.080

CrossRef Full Text | Google Scholar

Wang, N., Zhu, F., Chen, L., Chen, K. (2018). Proteomics, metabolomics and metagenomics for type 2 diabetes and its complications. Life Sci. 212, 194–202. doi: 10.1016/j.lfs.2018.09.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Y., Sun, Y., Wang, X., Lin, X., Sun, S., Shen, K., et al. (2016). Transcriptome shock in an interspecific F1 triploid hybrid of Oryza revealed by RNA sequencing. J. Integr Plant Biol. 58 (2), 150–164. doi: 10.1111/jipb.12357

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, J., Hu, F., Luo, K., Li, W., Liu, S. (2016). Unique nucleolar dominance patterns in distant hybrid lineage derived from Megalobrama Amblycephala×Culter Alburnus. BMC Genet. 17 (1), 150. doi: 10.1186/s12863-016-0457-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Sado, T., Hirt, M. V., Pasco-Viel, E., Arunachalam, M., Li, J., et al. (2015). Phylogeny and polyploidy: resolving the classification of cyprinine fishes (Teleostei: Cypriniformes). Mol. Phylogenet. Evol. 85, 97–116. doi: 10.1016/j.ympev.2015.01.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, H., Lin, Q., Luo, H. (2018). Applications of transcriptomics and proteomics in understanding fish immunity. Fish Shellfish Immunol. 77, 319–327. doi: 10.1016/j.fsi.2018.03.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, F., Zhong, H., Liu, G., Liu, S., Zhang, Z., Zhou, Y., et al. (2015). Characterization of vasa in the gonads of different ploidy fish. Genetics 574 (2), 337–344. doi: 10.1016/j.gene.2015.08.016

CrossRef Full Text | Google Scholar

Zdunek, R., Cichocki, A. (2006). Non-negative matrix factorization with quasi-newton optimization. Int. Conf. Artif Intell. Soft Comput. 870–879. doi: 10.1007/11785231

CrossRef Full Text | Google Scholar

Zhang, Z., Chen, J., Li, L., Tao, M., Zhang, C., Qin, Q., et al. (2014). Research advances in animal distant hybridization. Sci. China Life Sci. 57 (9), 889–902. doi: 10.1007/s11427-014-4707-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Meng, X., Yang, H. (2015). Jackknife empirical likelihood inference for the mean absolute deviation. Comput. Stat. Data Anal. 91, 92–101. doi: 10.1016/j.csda.2015.06.001

CrossRef Full Text | Google Scholar

Zhou, Y., Ren, L., Xiao, J., Zhong, H., Wang, J., Hu, J., et al. (2015). Global transcriptional and miRNA insights into bases of heterosis in hybridization of Cyprinidae. Sci. Rep. 5, 13847. doi: 10.1038/srep13847

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, Y., Zhao, T., Jiao, J., Chen, Z. (2019). The lifetime prediction of epoxy resin adhesive based on small-sample data. Eng. Fail Anal. 102, 111–122. doi: 10.1016/j.engfailanal.2019.04.007

CrossRef Full Text | Google Scholar


Results On Highly Expressed Genes


Table A1 A part of higly expressed genes of the six classes of fishes.


Table A2 The number of highly shared genes.


Table A3 A part of the base matrix W of the 1st-group samples L1.


Figure A1 Shared highly expressed genes for the second-group samples.

Keywords: transcriptomic profile, distant hybridization, optimization model, algorithm, classification, hybrids of fish, nonnegative matrix factorization

Citation: Wan Z, Tang J, Ren L, Xiao Y and Liu S (2019) Optimization Techniques to Deeply Mine the Transcriptomic Profile of the Sub-Genomes in Hybrid Fish Lineage. Front. Genet. 10:911. doi: 10.3389/fgene.2019.00911

Received: 02 May 2019; Accepted: 29 August 2019;
Published: 30 October 2019.

Edited by:

Wayne Aubrey, Aberystwyth University, United Kingdom

Reviewed by:

Jun Xiao, Hunan Normal University, China
Shihao Shen, University of California, Los Angeles, United States

Copyright © 2019 Wan, Tang, Ren, Xiao and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Zhong Wan,; Shaojun Liu,