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

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.


INTRODUCTION
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 Ren et al., 2016;Ren et al., 2017a;Ren et al., 2017b;Floriou-Servou et al., 2018;Li et al., 2018;Wang M. 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 Tan et al., 2019). For example, in Wang M. 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 subgenome 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 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, fulllength 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 multigeneration 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 . However, for the small-size samples with more than ten thousand features, such as the fulllength 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 largescale 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 R 20093×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 multigeneration hybrid fishes on the basis of their gene expression profile described by their transcriptome data.

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 highthroughput 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 multigeneration hybrid fishes in this research have been uploaded to https://github.com/TJY0622/TJY 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 A ∈ R n × m , we need to decompose A into two nonnegative matrices W and H, i.e. 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.

A WH
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: (2.5) 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, (2.7) Then, we denote the following two matrices the gradients of F(W, H) with respect to the matrices W and H, respectively: . (2.8) 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 And for k ≥ 1, we define four matrices: (2.10) 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 The following algorithm is developed to solve the subproblem (2.4) with the given H (k) .
Step 2 (Step length). Determine a step length α ρ k l l l a a l = = = max{ | , , , , ,} 0 1 2  such that α k satisfies the following inequality: Step 3  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 With the above preparation, we now develop an overall algorithm to solve Model (2.3) in the end of this section.
Step 5 (Projection of H). Replace H (k + 1) by Step 6 (Update). Set k := k + 1. Go to Step 1. Remark 1 Compared with the similar algorithms available in the literature , Algorithms 1 and 2 present a different computational procedure to solve Problem (2.1). Since the existing nonnegative matrix factorization methods depends on development of efficient solution algorithms, one of our contributions in this paper lies in developing Algorithms 1 and 2 to solve a sequence of subproblems like (2.4) and (2.5). Especially, in the section of result, we will implement them to solve the classification problem of distant multi-generation hybrid fishes based on their transcriptome profiles. Remark 2 In order to improve efficiency of Algorithm 2, before factorization of A, we conduct normalization of the sample data of fishes as follows.

RESULTS
In this section, in virtue of Model (2.3) and Algorithm 2, we present the results on classification of the distant multigeneration 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): 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 https://github. com/TJY0622/TJY.
For the sake of better understanding the inherent characteristics of the data, we take the 2nd-group samples with superscripts L 2 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 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: BSB L 2 , TC L 2 , TB F To further test robustness of the above trained results, given r = 6, we choose the 1st-group and the 3rd-group samples (with superscripts L 1 and L 3 , 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.
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 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: BSB BT BT . The accuracy rate is as high as 81.25%. Additionally, for the 5 species of fish (BSB, BT 2 F , BT F 3 , TB 2 F and TB F 3 ), 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.
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 BSB BT BT BT TB TB 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.
From Figure 3, it is easy to see that the sample of BSB L 1 is only weakly correlated with the two replicated samples BSB L 2 and BSB L 3 . Their correlation degree is even less than that between the samples of different fish BSB L 1 and TB F 2 . It can explain why BSB L 1 can not be clearly classified into the same class of BSB L 2 and BSB L 3 (revisiting the results in Table 5). Conversely, Figure 3 shows that in the 1st-group, the sample BSB L 1 has greater correlation with the other 3 samples: TB TB , and TC L 1 , which answers why the class of BSB L 1 can not be clearly identified in Table 3.
From Figure 3, we can also find out similar reasons for the unsatisfactory classification of BT BT BT    Table 5.
For the same reason of weaker correlation, in Tables 4 and 5, the three replicated samples of TB TC F 1 ( ) 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  , and BT F 3 , 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 F 1 (the 1st class), 45 shared highly expressed genes with BT F 3 (the 2nd class) and 12 shared highly expressed genes with TB F 2 (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 BT F 1 (3rd class) and TB F 1 (1st class), or between BT F 2 (3rd class) and TB F 2 (4th class), and there only exist 3 shared highly expressed genes between the BT F 3 (2nd class) and TB F 3 (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 TB TB F F 3 1 , (1st class), BT F 3 (2nd class) and TB F 2 (4th class). In contrast, there are only less than 45 shared highly expressed genes between the parent (BSB) and the hybrids ( BT F 3 ). Hybridization is considered as the rapidly driving forces that shape epigenetic modifications in plants and parts of lower vertebrate 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 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.

CONCLUSIONS
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 BT F 3 and TB F 3 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 (BT F 3 and TB F 3 ). In other words, both BT F 3 and TB F 3 are a good variety for the reproduction of fish. (5) Since there are no shared genes of high expression between BT F 1 and TB F 1 , they belong to two different classes (1st and 3rd classes). It implies that the reciprocal hybrids in the first generation of TC and BSB (BT F 1 and TB F 1 ) have larger biological distinction.

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.