Although the diffuse states generated by the above RWR process represent the underlying topological environment and intrinsic connectivity spectrum of each gene or disease node in the network, they may not be completely accurate due to the lowquality and high-dimensional nature of biological data. For example, a small number of missing or false interactions in the network can significantly affect the outcome of the diffusion process (Kim and Leskovec, 2011). It is often inconvenient to directly use high-dimensional diffusion states as topological features in prediction tasks.

To solve this problem, we used a dimensionality reduction method called DCA to reduce the dimensions of the feature space and obtain important topological features from the diffusion state. In addition, for multi-omics networks, DCA also performs very well. The key idea of DCA is to obtain an informative but low-dimensional vector representation. Similar to principal component analysis (PCA), which seeks the inherent low-dimensional linear structure of data to best interpret variances, DCA learns the low-dimensional vector representation of all nodes to best interpret their patterns of connection in heterogeneous networks. We will briefly describe the DCA framework below.

To achieve the purposes of noise reduction and dimensionality reduction, DCA uses the polynomial logic model represented by a low-dimensional vector to approximate the obtained diffusion state distribution, and it has far fewer dimensions than the original n-dimensional vector representing the diffusion state. Specifically, the probability of assigning node i to node j in the diffusion state is modeled as:

$$\hat{s}\_{ij} = \frac{\exp\left\{\mathbf{x}\_i^T \boldsymbol{\omega}\_j\right\}}{\sum\_{\mathcal{I}} \exp\left\{\mathbf{x}\_i^T \boldsymbol{\omega}\_{\mathcal{I}}\right\}}\tag{3}$$

In equation (3), xi, wi∈Rd , d ≪ n. We take wi as the context feature and xi as the node feature of node i, both of which describe the topological properties of the network. If xi and wi point in similar directions, we obtain a larger inner product. This means that node j may be frequently visited in a random walk starting from node i. DCA uses the obtained diffusion state S={ s1,⋯, sn } as input to optimize w and x of all nodes. The optimization method uses KL divergence, as shown in equation (4).

$$\min\_{\mathbf{w}, \mathbf{x}} \mathcal{C}(\mathbf{s}, \ \hat{\mathbf{s}}) = \min\_{\mathbf{w}, \mathbf{x}} \frac{1}{n} \sum\_{i=1}^{n} D\_{KL}(s\_i || \hat{\mathbf{s}}\_i) \tag{4}$$

DKL(⋅||⋅) is the KL divergence between the two distributions. We use w and x to represent this formula according to the definition of KL invergence and ^s.

$$C(s, \hat{s}) = \frac{1}{n} \sum\_{i=1}^{n}$$

$$\left[ H(s\_i) - \sum\_{j=1}^{n} s\_{ij} \left( \boldsymbol{w}\_i^T \boldsymbol{\mathcal{x}}\_j - \log \left( \sum\_{j'=1}^{n} \exp \{ \boldsymbol{w}\_i^T, \boldsymbol{\mathcal{x}}\_{j'} \} \right) \right) \right] \tag{5}$$

In equation (5), H(⋅) represents entropy. The objective function can find the low-dimensional vector representation of w and x using the standard quasi-Newton L-BFGS method. Although the obtained low-dimensional vector can effectively capture the network structure, we found that this optimization method is time-consuming.

To make the DCA framework more suitable for large biological networks, we use a more efficient method, clusDCA (Wang et al., 2015), which is based on matrix factorization, to decompose the diffusion states and obtain their low-dimensional vector representations. According to the definition, the following formula can be obtained:

$$\log \hat{s}\_{ij} = \boldsymbol{\kappa}\_i^T \boldsymbol{w}\_j - \log \sum\_{\boldsymbol{j'}} \exp \left\{ \boldsymbol{\kappa}\_i^T \boldsymbol{w}\_{\boldsymbol{j'}} \right\} \tag{6}$$

The first term corresponds to the low-dimensional approximation of ^sij. The second term is a normalization factor, ensuring that ^si is a well-defined distribution. By removing the second term, we relax the constraint that the elements in ^sij must add up to 1. Although the obtained low-dimensional approximation of the diffusion state is no longer a strictly valid probability distribution, it is found that these approximations are very close to the true distribution, and the effects of relaxation are negligible. Therefore, it can be simplified as:

$$\log \hat{s}\_{ij} = \mathbb{x}\_i^T \boldsymbol{\omega}\_j \,. \tag{7}$$

In addition, we use the sum of squared errors as the objective function, instead of minimizing the relative entropy between the original diffusion state and the approximate diffusion state.

$$\min\_{\mathbf{w},\mathbf{x}} \mathcal{C}(\mathbf{s},\ \hat{\mathbf{s}}) = \min\_{\mathbf{w},\mathbf{x}} \sum\_{i=1}^{n} \sum\_{j=1}^{n} \left(\boldsymbol{w}\_{i}^{T}\boldsymbol{x}\_{j} - \log \boldsymbol{s}\_{ij}\right)^{2} \tag{8}$$

The obtained objective function can be optimized by singular value decomposition (SVD). To avoid taking the logarithm of 0, we add a small positive number 1 n to sij. The calculation formula of the logarithm diffusion state matrix L is as follows:

$$L = \log(\mathcal{S} + \mathcal{Q}) - \log(\mathcal{Q})\,. \tag{9}$$

In equation (9), S∈Rn×n , Q∈Rn×n and Qij = 1 n, for ∀i, j. Using the singular value decomposition method, we decompose L into three matrices:

$$L = U\Sigma V^T\tag{10}$$

In equation (10), U∈Rn×n , V∈Rn×n , S∈Rn×n and S is a diagonal singular value matrix. To obtain the low-dimensional vectors wj and xi in d dimensions, we simply select the first d singular vectors Ud, Vd, and Sd. Each row of matrix X = ½x1, …, xn- T represents the low-dimensional eigenvector corresponding to each node in the network. In matrix W = ½w1, …, wn- T , each row represents the corresponding vector of the context feature. The formulas for calculating X and W are as follows:

$$X = U\_d \Sigma\_d^{1/2}, \quad W = V\_d \Sigma\_d^{1/2}.\tag{11}$$

#### Denoising Autoencoder

Autoencoder is an unsupervised neural network model. It learns the implicit features of input data, which is called "coding." At the same time, the original input data can be reconstructed with the learned new features, which is called "decoding." Intuitively, autoencoder can be used for reducing feature dimensionality, like principal components analysis (PCA), but with stronger performance than PCA because the neural network model can extract more effective new features.

The denoising autoencoder adds noise to the input x to obtain ~x, and after training, it obtains a noiseless output z, as shown in Figure 2.

This prevents the autoencoder from simply copying the input to the output, so as to extract useful patterns in the data and improve the weight robustness. Noise can be either pure gaussian noise added to the input or randomly discarding a feature at input layer, similar to dropout. The specific equation for calculating z is as follows:

$$\begin{aligned} \mathcal{Y} &= f\left(\tilde{\mathcal{X}}\mathcal{W}\_1 + b\_1\right) \\ \mathcal{Z} &= \mathcal{g}\left(\mathcal{Y}\mathcal{W}\_2 + b\_2\right) \end{aligned} \tag{12}$$

In addition, network parameters are trained to minimize reconstruction errors, namely:

$$\min L\_H(\varkappa, z) = \min \left| \left| \varkappa - z\_{\mathbb{P}} \right. \right. \tag{13}$$

#### Pathogenic Gene Prediction Method Standard Inductive Matrix Completion

In the gene–disease association matrix P ∈ RNgNd , each row represents a gene ID and the number of genes is Ng. Each column represents a disease phenotype and the number of diseases is Nd. If Pij = 1, this means that gene i is related to disease j, and Pij = 0 means that the relationship between gene i and disease j is uncertain. Based on the most successful and deeply studied matrix completion method in the recommender systems, the IMC algorithm was used to complete the task of learning gene– disease associations. The advantage of this is that this method is inductive, and it can achieve the prediction of new genes or diseases that have rarely been studied.

IMC assumes that the association matrix has a low rank, with the goal of recovering Z using the observed values of P and the eigenvectors of genetic diseases, as shown in Figure 3.

The eigenvector matrix of Ng genes is represented by X ∈ RNgfg , and the eigenvector of gene i is represented by xi ∈ Rfg .

Similarly, Y ∈ RNdfd is used to represent the eigenvector matrix of Nd diseases, and yi ∈ Rfd is used to represent the eigenvector of disease j. The inductive matrix completion problem is to recover a low-rank matrix Z by using the known association W+ from the gene–disease association matrix P. We established a bilinear function to learn the projection matrix Z between the gene space and the disease space to predict the interaction between unknown genes and diseases. We modeled the matrix P as XZYT ≈ P. Then, we used the following formula to measure the probability of pairwise interaction score between gene i and disease j, and the higher the score(i, j) value, the more likely gene i and disease j interact.

$$score(i, j) = \varkappa\_i Z y\_j^T \tag{14}$$

There is usually a significant correlation between spatially close eigenvectors of genes or diseases, which can greatly reduce the number of effective parameters needed to model gene– disease interactions in Z. To consider this problem, we applied a low-rank constraint on Z and learned only a few potential factors. Let Z = GHT , where G ∈ Rfgk , H ∈ Rfdk , and k is small. This low-rank constraint not only alleviates the overfitting problem, but also facilitates the process of optimizing the calculation (Wang et al., 2015). The optimization problem of low-rank constraint is NP-hard on the original matrix Z. One standard method of relaxing the low-rank constraint is to minimize the trace norm, that is, the sum of the singular values. Minimizing the trace norm of Z = GHT is equivalent to minimizing 1 2 (‖G‖2 F + ‖H‖2 F). The decomposition of Z into G and H solves the following optimization problems by alternating minimization. A common choice for the loss function ℓ is the square loss function. l is the regularization parameter.

$$\min \sum\_{(i,j)\in\Omega^{+}} \ell\left(P\_{ij}, \boldsymbol{\chi}\_{i}^{T} \boldsymbol{G} \boldsymbol{H}^{T} \boldsymbol{\mathcal{Y}}\_{i}\right) + \frac{\lambda}{2} \left(\|\boldsymbol{G}\|\_{F}^{2} + \|\boldsymbol{H}\|\_{F}^{2}\right) \tag{15}$$

#### Improved Inductive Matrix Completion

To optimize the objective function, we introduce the idea of PU-Learning. Although we predicted positive examples from unknown relationships, that is, candidate disease-causing genes, it was undeniable that these unknown genes-disease pairs may be unrelated. Therefore, unknown association relationship information was added to the learning process as a negative example, and the objective function was as follows:

$$\begin{split} \min \sum\_{(i,j)\in\Omega^{\*}} \ell\left(P\_{\vec{\eta}}, \mathbf{x}\_{i}^{T} GH^{T} \boldsymbol{\wp}\_{i}\right) + \alpha \sum\_{(i,j)\in\Omega^{\*}} \ell\left(P\_{\vec{\eta}}, \mathbf{x}\_{i}^{T} GH^{T} \boldsymbol{\wp}\_{i}\right) \\ + \frac{\lambda}{2} \left(\|G\|\_{F}^{2} + \|H\|\_{F}^{2}\right) \end{split} \tag{16}$$

We represent the unknown association in the gene–disease association matrix P as W− . The key parameter a < 1 because the penalty weight of the known relationship must be greater than the unknown relationship. Finally, equation (14) was still used to calculate the interaction score between gene i and disease j. The scores are sorted in descending order, and the first k genes were selected as candidate pathogenic genes for the corresponding disease.

#### DATA SETS AND FEATURES

The data sets used in this paper can be divided into three categories: gene–disease association data, gene characteristic data, and disease characteristic data.

#### Gene–Disease Associations

The known gene–disease association data that we used were from the OMIM database, which contained 12,331 genes, 3,209 diseases, and 3,954 known gene–disease associations (the total number of nonzero elements in the gene–disease association matrix). It can be seen that the data in the incidence matrix are very sparse, with more than 90% of the columns having only one nonzero item and 70% of the rows having no nonzero elements.

#### Gene Characteristics

Gene characteristics were obtained by processing four different data sources through compact feature learning (Compact Feature Learning). The first source of gene characteristics was gene microarray data, which contained 8,755 genes and 4,536 characteristics. First, we linearly transformed the expression range of each gene to between 0 and 1. Because these characteristics are highly correlated, we used four layers of denoising autoencoder to reduce the dimensionality of the data, and the number of cells in each hidden layer was 3,000-800-300-100, respectively. Moreover, gaussian noise with a noise factor of 0.2 was added to the input data, and sigmoid was used to activate each layer. The model was optimized with Adam, and epoch was 100.

The second source of gene characteristics was from homologous gene phenotypic associations in eight other species, which were more abundant than in studies of human genetic diseases. The data used in the experiment are shown in Table 1. The features were extracted by two-layer denoising autoencoder with the following specific parameters: the number of nodes in each layer is "200–100," the corruption level of data is 0.2, the activation function is sigmoid function, the batch size is set as 150, and the model is optimized by Adam.

In addition, the data on interactions between genes can also be used as a part of the characteristics of genes. We integrated two networks, HumanNet (Lee et al., 2011) and STRING (Szklarczyk et al., 2017), for unified analysis. These two sets of data represent gene–gene interaction networks, but there are differences between them (Kuang et al., 2018). The integrated analysis of different sets of data can verify each set, and they can help to validate each other and expand understanding the potential rules. We used the RWR and DCA methods to fuse two networks to extract gene features. We set the restart probability to 0.05 and extracted the 600 dimensional gene characteristics. Finally, the gene characteristics used in the model were 800 dimensions.

#### TABLE 1 | Species Details.

#### Disease Characteristics

The disease characteristics are mainly derived from two data sources: the disease similarity network MimMiner and clinical manifestation data of the disease, as well as a large amount of data from analysis of the medical literature.

MimMiner data are processed by literature (van Driel et al., 2006) and are freely available online. This data set has been applied in gene prioritization methods (Vanunu et al., 2010; Singh-Blom et al., 2013; Natarajan and Dhillon, 2014). RWR and DCA were used to extract 100-dimension disease features in the disease similarity network, and the restart probability was set as 0.05.

Another disease feature that we incorporated was from the OMIM disease webpage. We paid special attention to the clinical features and clinical management of webpages. We obtained disease features through text mining. We used PCA to reduce the dimensions of feature space and retained the first 100 principal components. Finally, we obtained 200-dimension disease characteristics.

#### EXPERIMENT

Evaluation Indexes and Methods introduces the evaluation indexes and methods of the experiment. Parameter Settings describes the influence of important parameters in the experiment. In Global Performance, the global performance of the experiment is compared. Prediction of New Genes and New Diseases compares the ability to predict new genes and new diseases. Newly Discovered Genes compares the ability to predict newly discovered associations.

#### Evaluation Indexes and Methods

In the experiment, to quantitatively evaluate our method and compare it with the most advanced disease-causing gene prioritization methods, we used a cross-validation strategy to measure gene recovery. We divided the known gene–disease pairs into three groups of the same size. The associations in one group were hidden, and the associations in the remaining two groups were used as training data, repeated three times to ensure that each group was hidden only once. For each disease in our data set, we ranked all of the genes according to the degree to which they were associated with the disease. The first r genes were taken as candidate pathogenic genes for corresponding diseases; namely, the top-r ranking method was used. The performance of the algorithm was analyzed by comparing the recall and precision of each method under different thresholds r, usually r ≤ 100. The formula for calculating this was as follows:

$$Recall = \frac{\text{TP}}{\text{TP} + \text{FN}} \tag{17}$$

$$Precision = \frac{\text{TP}}{\text{TP} + \text{FP}} \tag{18}$$

Recall rate refers to the proportion of positive cases correctly judged by the model relative to all positive cases (TP+FN) in the data set. FN represents the data that are mistaken as negative cases by the model but are actually positive cases. The precision rate is the proportion of true positive cases (TP) relative to all positive cases (TP+FP) judged by the model (Xiong et al., 2012; Xu et al., 2017; Cheng et al., 2019; Cheng et al., 2019).

To further confirm the value of our approach, we also used the mean percentile ranking (MPR), an evaluation index based on recall, to evaluate the performance of the algorithm. This evaluation index has been applied in recommendation algorithm and analyses of the performance for predicting drug-targets (Hu et al., 2008; Johnson, 2014; Li et al., 2015; Ding et al., 2017; Hao et al., 2019; Liu et al., 2019b; Liu et al., 2019c; Zeng et al., 2019) and disease biomarkers (Chen et al., 2016; Zeng et al., 2016; Hong et al., 2019; Xu et al., 2019). For each disease, the genes were ranked in descending order according to the calculated gene–disease predictive value. The average ranking of the true and established associations among them is the final result. Here, rankji can be used to represent the percentile ranking (PR) of gene j and disease i. rank ji = 0% indicates that disease i is most likely to interact with gene j. Similarly, rank ji = 100% indicates that disease i has the lowest probability of interacting with gene j. Therefore, the definition of MPR is as follows:

$$MPR = \frac{\sum\_{i=1}^{N\_D^t} R\_i}{N\_D^t} \tag{19}$$

Nt T represents the number of diseases in the test set, and the formula for calculating Ri is as follows:

$$R\_i = \frac{\sum\_{j=1}^{N\_T^t} \text{rank}\_{ji}}{N\_T^t} \tag{20}$$

Nt T represents the number of genes in the test set for current disease i. It is important to emphasize that lower MPR values are preferable because they indicate that our approach has a higher probability, which means that the model works better. Conversely, a higher MPR indicates a lower likelihood of gene interactions with disease. Clearly, the randomly generated list is expected to have an MPR of 50%. Using this measure, we can obtain a list of recommended candidate pathogenic genes, where the recommended optimal prediction is used for higher priority experimental validation.

#### Parameter Settings

The key parameters of PUIMCHIF are the rank k of matrix Z∈R800×200, the regularized parameter l, and the penalty weight a for the unknown relation. As can be seen from Figure 4, the performance of the PUIMCHIF method increases with the increase of k. When k = 100, 150, and 200, the three curves are very close. In the following experiment, the PUIMCHIF method uniformly set parameters as follows: k = 200, l = 0.02, and a=0.0035.

As mentioned earlier, our approach features four improvements over the original IMC approach. For Figure 5, recall, precision, and MPR were used to analyze the effect of our improved method. The four experimental results in the figure represent (a) the initial experimental results of the original IMC method, (b) the results

of extracting features by using RWR and DCA, instead of PCA, for the network data of diseases and genes, (c) the prediction results of adding STRING data to the gene interaction network, and (d) the experimental results of each index of the PUIMCHIF method.

We found that using RWR and DCA can better extract the gene–gene and disease–disease relationships, and helps to improve the prediction of candidate pathogenic genes. Meanwhile, it was also found that the protein interaction network STRING improved the prediction recall rate to 47.45%, and the MPR value also decreased significantly. Using denoising autoencoder to represent the characteristics of genes and diseases, and introducing the idea of PU-Learning into the inductive matrix completion can further improve the predictive performance.

#### Global Performance

In this experiment, the threefold cross-validation method was used to compare the overall performance of the proposed method with CATAPULT, Katz, and IMC. As shown in Figure 6A, the vertical axis gives the probability of recovering the true gene association in the top-r prediction of different r values on the horizontal axis. The experimental results show that the PUIMCHIF algorithm proposed in this paper has a much higher probability of recovering true gene associations under different thresholds than the other methods. Figure 6B presents the precision–recall curve.

In addition, Table 2 shows the results of three evaluation indexes for each method when the threshold r=100. It is worth mentioning that a smaller value of MPR is associated with a higher probability and a better effect. It can be seen that the MPR value of PUIMCHIF is the lowest and the recall rate reaches 50%, while the best method among other methods, IMC, is only 25%, that is, the recall is doubled. The precision rate was also twice that of Katz which is the best method of other methods, reaching 4.87%. The overall performance of PUIMCHIF has been further improved, confirming the superiority of our method.

# Prediction of New Genes and New Diseases

# Prediction of New Genes

One problem affecting prioritization assessments is that well-related genes and diseases tend to be more predictable and therefore tend to generate inflated recall rates. Here, we focued only on genes that are known to have a single association in the gene-disease association data set. In other words, we selected the gene corresponding to the row with only 1 non-zero element in the gene-disease association matrix as the validation set, and hided these known associations in the training process. After repeated three-fold cross validations, Figure 7A shows the predictive power of different methods within the threshold r < 100. The Y-axis represents the probability of a true known single gene association hidden during recovery training.

TABLE 2 | Experimental Results with Threshold r =100.

Table 3 shows the specific experimental results of each method when r = 100. For the prediction of new genes, although the precision rate was slightly lower than Katz, the recall rate of our PUIMCHIF was significantly higher than other methods, reaching 40.7% when the recall rate of IMC method was only 13.7%. At the same time, we found that using the MPR index to evaluate the results, the PUIMCHIF method was only 13.5%, much lower than Katz and CATAPULT. This also shows that our method is more reliable.

#### Prediction of New Diseases

Similar to the prediction of new genes, we only considered diseases with a single known association in the gene-disease association data set as the validation set, that is, diseases corresponding to the columns with only 1 non-zero element in the gene-disease association matrix, and hided these known associations during training. Similarly, a three-fold crossvalidation analysis was used, and the results are shown in Figure 7B. The probability that the proposed method could recover the true association of new diseases reached 48%, which was a significant improvement compared with other methods. Moreover, the MPR value of our method was lower than that of other methods, and the precision rate was nearly 2.7 percentage points higher than that of IMC method. As can be seen from Table 4, PUIMCHIF method is superior to other methods in three evaluation indexes.

#### Newly Discovered Genes

Cross-validation of retrospective data can lead to overly optimistic performance estimates. For example, certain gene interactions may be found because of associations with specific diseases being evaluated. Although the association itself is hidden, other features are contaminated by this information. Therefore, the use of recently reported associations to assess gene prioritization tools is unbiased in this assessment.

We trained all methods using all the gene associations of the 3,209 OMIM diseases collected. We found 162 newly discovered associations, of which 83 genes had no known associations previously. Thus, the assessment of new associations also helps determine the ability of methods to recommend new genes. The ranking performance of each method in 162 new associations is shown in Figure 8. We can see that the IMC method is superior to other methods in the range of threshold 6 ≤ r ≤ 100.

#### CONCLUSION

In this paper, a PU induction matrix completion algorithm based on heterogeneous information fusion, PUIMCHIF, was proposed to predict gene–disease associations. Based on the specific advantages of IMC method, PUIMCHIF can predict new genes and diseases, and has good predictive performance. In addition, because closely connected or functionally similar genes are more likely to cause the same or similar diseases, we

TABLE 3 | Prediction of New Genes with Threshold r =100.

TABLE 4 | Prediction of New Disease with Threshold r =100.

constructed low-dimensional feature representations of genes and diseases from various data sources such as STRING using the compact feature learning method, which effectively alleviated the impact of data sparsity. Although there is no evidence that genes are unrelated to diseases in the data set, it is clear that most of the unknown associations are negative. PUIMCHIF conducts biased learning by treating unlabeled data as negative cases and constraining the penalty weight of known relationships to be greater than that of unknown relationships. Compared with the existing prediction methods, the PUIMCHIF method can significantly improve the prediction results regarding recall rate, precision rate, and MPR. According to the evaluation

#### REFERENCES

index of MPR, the experimental results of the PUIMCHIF method that we proposed are the lowest; that is to say, the candidate genes given by our algorithm have a higher priority for validation by biological experiments.

#### DATA AVAILABILITY STATEMENT

The datasets analyzed in this article are not publicly available. Requests to access the datasets should be directed to<EMAIL_ADDRESS>
#### AUTHOR CONTRIBUTIONS

CW initiated the idea, conceived the whole process and drafted the manuscript. JZ and XW implemented the experiments and designed the figures. KH helped with data analysis and revised the manuscript. MG and finalized the paper. All authors have read and approved the final manuscript.

#### ACKNOWLEDGMENTS

The work was supported Natural Science Foundation of China (61872114, 61571163, 61532014, and 61871020), and the National Key Research and Development Plan Task of China (No. 2016YFC0901902). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. We thank Liwen Bianji, Edanz Group China (www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.

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.

Copyright © 2020 Wang, Zhang, Wang, Han and Guo. 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.

# Gene-Focused Networks Underlying Phenotypic Convergence in a Systematically Phenotyped Cohort With Heterogeneous Intellectual Disability

Yan Wang1,2,3† , Li-Na Zhu1,2,3† , Xiu-Wei Ma1,2,3† , Fang Yang4 , Xi-Lin Xu4 , Yao Yang1,2,3 , Xiao Yang1,2,3, Wei Peng1,2,3, Wan-Qiao Zhang1,2,3, Jin-Yu Liang5 , Wei-Dong Zhu5 , Tai-Jiao Jiang4,6, Xin-Lei Zhang7 and Zhi-Chun Feng1,2,3 \*

#### Edited by:

Tao Huang, Shanghai Institutes for Biological Sciences (CAS), China

fbioe-08-00045 February 6, 2020 Time: 16:29 # 1

#### Reviewed by:

Kuo Yang, Beijing Jiaotong University, China Shur-Jen Wang, Medical College of Wisconsin, United States

#### \*Correspondence:

Zhi-Chun Feng<EMAIL_ADDRESS><EMAIL_ADDRESS>
†These authors have contributed equally to this work

#### Specialty section:

This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Bioengineering and Biotechnology

Received: 18 November 2019 Accepted: 21 January 2020 Published: 07 February 2020

#### Citation:

Wang Y, Zhu L-N, Ma X-W, Yang F, Xu X-L, Yang Y, Yang X, Peng W, Zhang W-Q, Liang J-Y, Zhu W-D, Jiang T-J, Zhang X-L and Feng Z-C (2020) Gene-Focused Networks Underlying Phenotypic Convergence in a Systematically Phenotyped Cohort With Heterogeneous Intellectual Disability. Front. Bioeng. Biotechnol. 8:45. doi: 10.3389/fbioe.2020.00045 1 BaYi Children's Hospital, The Seventh Medical Center of PLA General Hospital, Beijing, China, 2 National Engineering Laboratory for Birth Defects Prevention and Control of Key Technology, Beijing, China, 3 Beijing Key Laboratory of Pediatric Organ Failure, Beijing, China, 4 Suzhou Institute of Systems Medicine, Chinese Academy of Medical Sciences, Suzhou, China, 5 The Second People's Hospital of Aohanqi, Inner Mongolia, China, 6 Center of Systems Medicine, Institute of Basic Medical Sciences, Chinese Academy of Medical Sciences & Peking Union Medical College, Beijing, China, 7 Suzhou Geneworks Technology Co., Ltd., Suzhou, China

The broad spectrum of intellectual disability (ID) patients' clinical manifestations, the heterogeneity of ID genetic variation, and the diversity of the phenotypic variation represent major challenges for ID diagnosis. By exploiting a manually curated systematic phenotyping cohort of 3803 patients harboring ID, we identified 704 pathogenic genes, 3848 pathogenic sites, and 2075 standard phenotypes for underlying molecular perturbations and their phenotypic impact. We found the positive correlation between the number of phenotypes and that of patients that revealed their extreme heterogeneities, and the relative contribution of multiple determinants to the heterogeneity of ID phenotypes. Nevertheless, despite the extreme heterogeneity in phenotypes, the ID genes had a specific bias of mutation types, and the top 44 genes that ranked by the number of patients accounted for 39.9% of total patients. More interesting, enriched co-occurrent phenotypes and co-occurrent phenotype networks for each gene had the potential for prioritizing ID genes, further exhibited the convergences of ID phenotypes. Then we established a predictor called IDpred using machine learning methods for ID pathogenic genes prediction. Using10-fold crossvalidation, our evaluation shows remarkable AUC values for IDpred (auc = 0.978), demonstrating the robustness and reliability of our tool. Besides, we built the most comprehensive database of ID phenotyped cohort to date: IDminer http://<IP_ADDRESS>: 3100/IDminer/, which included the curated ID data and integrated IDpred tool for both clinical and experimental researchers. The IDminer serves as an important resource and user-friendly interface to help researchers investigate ID data, and provide important implications for the diagnosis and pathogenesis of developmental disorders of cognition.

Keywords: intellectual disability, phenotypic convergence, gene-focused networks, co-occurrent phenotype, machine learning, pathogenic genes prediction

# INTRODUCTION

fbioe-08-00045 February 6, 2020 Time: 16:29 # 2

Intellectual disability (ID), also known as mental retardation, is characterized by significant impairment in cognition. The patients with ID usually have the obvious shortcomings of adaptive behavior before the age of 18, and a high incidence rate, 1–3%, making it a worldwide social problem (Maulik et al., 2011; Mefford et al., 2012). It can occur in isolation or in combination with congenital malformations or other neurological features such as epilepsy, congenital malformations, sensory impairment, and autism spectrum disorders (ASD), and its severity (mild, moderate, severe, and profound) is highly variable (Vissers et al., 2016). The heterogeneity of phenotypes poses additional challenges for understanding the complex etiology, with contributions by environmental factors, perinatal hypoxia, and genetic factors. In recent years, genetic factors including chromosomal abnormalities, single and multiple gene mutations have found to become increasingly prominent for the disease (Gilissen et al., 2014; Lelieveld et al., 2016; Reichenberg et al., 2016). With the increasing number of ID cases identified in clinics, its phenotypes have found to be extremely heterogeneous. Previous studies found that patients with identical mutations in a single gene could give rise to different phenotypes (Hoischen et al., 2014). As the limitations of detection technologies and the heterogeneity of ID genes and phenotypes, many patients still lack appropriate diagnosis.

In the past 10 years, a large number of studies have been carried out in order to explore the genetic mechanism of ID (Gécz et al., 2009; Ellison et al., 2013). In particular, the development of second-generation sequencing technology facilitates the rapid investigation of more DNA samples from ID cases (Rauch et al., 2012; Gilissen et al., 2014). This led to an expansion in the number of genes associated with ID. Having mass data about ID genes, clinical phenotypes, and pedigrees available in the public domain could shed insights into ID mechanisms. A previous report suggests that ID genes are substantially enriched with coexpression, protein-protein interactions, and specific biological functions. Furthermore, they also revealed combinations of typical phenotypes within process-defined groups of ID disorders by clusters of ID genes with significantly elevated biological coherence (Kochinke et al., 2016). This suggests that ID genes and phenotypes have their own characteristics, and these data can be used to define mechanisms of ID and may improve the diagnosis of patients.

In this study, the ID genes, phenotypes, and pedigrees were extracted manually and analyzed and then integrated to build a standard ID database IDminer, which analyzed the phenotypes, genes, families and their relationships based on the individual patient. Furthermore, the candidate pathogenic genes for ID patients could be prioritized based on the molecular feature of ID genes and the genes specific phenotypes and phenotypic pairs. Furthermore, the similarity between patients was also evaluated via clinical features and could help patients with effective intervention. Importantly, the curated data including ID phenotypes, genes and pedigrees, their integrated analysis and their applications are accessible online via http://<IP_ADDRESS>: 3100/IDminer/.

# MATERIALS AND METHODS

# Analysis of Specific Phenotypes and Phenotypic Pairs

Each pathogenic gene could be associated with multiple patient samples, and each patient may have different phenotypes. For each gene, the specific phenotypes were obtained with the enrichment analysis using the hypergeometric distribution. A gene could correspond to multiple patients. For each patient, any two of their phenotypes formed a phenotype pair, referred to as co-occurrence. A phenotype pair could appear in N patients (N represents the frequency of phenotypic pairs). In situations with a single gene affecting multiple cases, multi-phenotypic pairs and their frequencies were obtained. For each phenotypic pair, we analyzed whether the co-occurrence was enriched in the affected patients or not.

## Construction of Co-occurrence Network

For constructing a co-occurrence network, all phenotypic pairs with a P-value = 0.05 for at least one gene were built as a non-directional network. In this network, each node represents a phenotype and the node size indicates the frequency of the phenotype in the database, while the edges denote significantly enrichment between phenotype pairs. Then the modules were extracted with the R igraph package.

# Phenotype-Based Samples Similarity Analysis

The same phenotype may appear in different patient samples. Based on the number of the same phenotypes between these samples, similarity scores between pairs of patient samples was calculated.

# The Phenotype Converting Tool

The tool was used to calculate the similarity between the users' input phenotypes and the 2,075 standard phenotypes in this website. The python module named FuzzyWuzzy was used to calculate the similarity score [0,100]. The higher the score, the more similar the two phenotypes.

#### Supervised Machine Learning Prediction

In this study, the supervised machine learning method, Support Vector Machine (SVM), was employed for ID pathogenic genes prediction. The R language interface of LIBSVM was used to construct the SVM-based pathogenic predictors. The radial basis function was chosen as the kernel function, and the other parameters were set at the default. A prediction model was trained using repeated 10-fold cross-validation of the training dataset, and their predictive performance was evaluated in the independent test dataset.

# Web Interface Configuration

The interface has two main parts: one part displayed the ID knowledge base data and the search results, while the other displayed the input and results of the analysis tool. Through the search box on the main page, users could search for a gene

or a phenotype. Through the tools button in the main menu, users could enter the analysis interface, and according to the given phenotypes and genes, the ID genes were identified, and the association between the genes and their phenotypes were visualized. The web service was mainly based on java server pages, JavaScript, R, Python, Ajax, Apache, and MySQL.

#### RESULTS

#### Data Curation

We first employed the keywords, such as ID, mental retardation, developmental delay, cognitive impairment, developmental disability, and learning disability to accomplish the literature searches by using PubMed. Then the literature was filtered through the artificial proofing method, and the ID-related papers and genes were retained. The text mining method was used to mark phenotypes in the literature using the HPO1 database phenotypic information as a reference. Then the gene name, mutation site, and phenotypes were curated manually (**Figure 1**). Based on the sample description in the literature, the family information of the samples were also collated from the HGNC (HUGO Gene Nomenclature Committee) database according to the acquired ID-related gene name information, such as gene alias, chromosome localization, corresponding OMIM ID, and Ensembl ID, and the biological function and pathway information for these genes were marked simultaneously through GO2 and KEGG3 databases.

#### The Landscape and Convergence of ID Genes

Through 1174 ID papers, we obtained a total of 3803 samples with 2075 phenotypic descriptions, that were caused by 704 ID genes. Among these genes, there are 3848 mutations, containing 1793 missense/non-sense mutations, 182 splicings, and 610 indels. We found that the majority of the genes were identified in less than 10 patients, and 305 genes (43.3%) found in only one patient and 103 genes (14.8%) in two patients (**Figure 2A**). Also, a small set of genes caused more patients than other genes, as shown in **Figure 2B**, the top nine genes ranked by the number of patients accounted for 14.9% of the total patient group, and the top 44 genes included 39.9% of patients. Moreover, our analysis also showed some ID genes had the dominant mutation types (**Figure 2C**). For the top 57 genes ranked by the number of ID patients, the majority of mutations of patients harboring mutated MECP2, HUWE1, and CREBBP are gross insertions. In addition, the predominant mutation type of patients with mutated THOC2, KIF1A, KDM5C, IQSEC2, SLC6A8, TBC1D24, MAN1B1, YAP1, GRIN2B, PAK3, NALCN, CLPB, and GRIN1 genes are missense/non-sense mutations, while deletions are mainly found in patients harboring SOX4, NRXN1, FMR1, MEF2C, OPHN1, PQBP1, AUTS1, MYT1L, CNTNAP2, MAPT, and TUSC3 genes. Importantly, the mutation types of 47 of the top 57 genes contained gross insertions (most duplications) and missense/non-sense, suggesting that both deletion and overexpression of these genes were likely to cause ID disease. These findings suggested that despite the diversity of ID genetic variation, most ID patients are caused by a small number of genes based on its genetic bias and convergence.

#### The Heterogeneities of ID Phenotypes

Among the patient cohort, 637 (16.6%) patients have a unique phenotype, while 901 (23.7%) patients have more than ten phenotypes (**Figure 3A**). Also, our data showed that the number of phenotypes for each patient had a positive correlation with the number of the patients, which showed a significant linear relationship (Spearman P-value < 0.001, **Figure 3B**) and indicated the heterogeneity of the ID phenotypes. Additionally, HPO structure analysis found the accompanying phenotypes of ID were also widely distributed, including symptoms in many parts of the body (**Figure 3C**). For these phenotypes, as shown in **Figure 3D**, the top 50 phenotypes ranked by the number of patients exhibited that the ID was usually accompanied by other mental diseases, such as seizure, epilepsy, microcephaly, ataxia, microcephaly and autism, abnormal behaviors containing hypotonia, strabismus, sleep disturbance, constipation, delayed or absent speech, motor delay, hyperactivity, feeding difficulties and inability to walk, and dysmorphism about spine, face, stature, and cryptorchidism. These results showed that the phenotypes of ID patients had extreme heterogeneity.

#### The Convergences of ID Phenotypes

The phenotypes that were converged for each gene based on the fact that intra-similarity between patients caused by one gene were more than inter-similarity between different genes' patients (**Figure 4A**) and the phenotypes in patients caused by the identified mutations in the same family had more similarity than other families (**Figure 4B**). To better understand the convergence of the ID clinic features, we first obtained the specific phenotypes for each ID gene with enrichment analysis. A total of 143 phenotypes, appearing in at least five patients caused by the same gene, were enriched in some genes' patients (**Figure 4C**). Importantly, among the phenotypes, 47 appeared in only single gene's patients and accounted for 30 genes, which could help to diagnosis the patients caused by the genes (**Figure 4C**). To illustrate the relationships between phenotypes, we also investigated the situation of two phenotypes could be co-occurred in one patient, and the co-occurrence phenotypes were recorded as "phenotypic pairs." We analyzed these phenotypic pairs presented in patients with an enrichment analysis. Interestingly, we found that most enriched phenotypic pairs were specific for a single gene. Like single phenotype analysis, phenotypic pairs made it easy to diagnosis patients with 82 ID genes (**Figure 4D**).

#### Gene-Focused Network for Phenotype Enrichment

Then we analyzed the network diagram of the phenotypic pairs for each gene, which revealed the gene-focused network (**Figure 5A**) and three typical sub-networks (**Figure 5B**). The

1http://human-phenotype-ontology.github.io

2http://geneontology.org/

3http://www.genome.jp/kegg/

first type of sub-network was radial, indicating that most of the phenotypes co-occurred with another one phenotype (like gene ZNF711). The pathogenic genes with the first type of sub-network may have a core phenotype, or an important phenotype that appeared more frequently, and it illustrated that there are strong association between the core phenotype and the biological function. The second type of sub-network was dense, and the phenotypes co-occurred with each other (like gene PIGO). The pathogenic genes with the second type of subnetwork often result in a set of concurrent phenotypes. In this case, the prediction of pathogenic genes by phenotype may be more accurate. The third type of sub-network was the mixed state of the above two types (like gene MECP2). With the third type of sub-network of pathogenic genes, the mutations are usually more extensive, the phenotypes are complex, and one independent group phenotypes is often insufficient to reveal the pathogenic genes information. Our analysis showed that the cooccurrence network of each gene had its own characteristics, and the phenotypes in the co-occurrence network of each gene are different. And the co-occurrence networks of different genes had commonality in their structural similarity. Analysis of co-occurring networks further illustrated the phenotypic conservation relative to genes, despite the heterogeneity of phenotypes. Based on the above discoveries, we inferred that the pathogenic genes for patients could be achieved by analyzing specific phenotypes and phenotypic pairs. Our analyses indeed showed that the more the patients' phenotypes, the more accurate the prediction of pathogenic genes (**Figure 5C**). Furthermore, given more phenotypes, the predicted pathogenic genes incline to have a more significant P-values (**Figure 5D**). These results showed that phenotypic analysis could reveal the convergences of ID phenotypes and be used for clinical pathogenic gene analysis.

# Pathogenic Gene Prediction

Support Vector Machine is one of the most widely used machine learning algorithms in computational biology. It was previously used for predicting virulent proteins in bacterial pathogens (Garg and Gupta, 2008), the clinical outcome from cancer patients (Yeoh et al., 2002) and gene interactions in genetic diseases (Upstill-Goddard et al., 2013). As shown in **Supplementary Figure S1**, developed SVM-based predictor, a 10-fold crossvalidation was employed on the training datasets for model selection purpose (**Figure 6A**), and the final performance of the predictor was measured on the independent testing dataset (Ortiz-Gonzalez et al., 2018) compared with other ID pathogenic gene prediction models (Yang et al., 2015; Stelzer et al., 2016; **Figures 6B,C**). The receiver operating characteristic curve (sensitivity against 1-specificity) was used to measure the prediction performance under different decision thresholds, and the area under the curve (AUC) was calculated as the main performance evaluation metric. For calculating variable importance for prediction, 100 sets of independent training were performed using different random seed. The median of variable importance obtained in each training was used as a representative value (**Supplementary Figure S2**).

# Database and Tool for ID Research and Diagnosis

In order to represent the ID data and the analysis tools for ID research and diagnosis, the IDminer system was designed. The database included a number of components, including a knowledge base for intellectual disabilities, specific phenotypes and phenotypic pairs for genes, co-occurrence networks, and analysis tools for converting phenotypes to standard phenotypes

and exploring the expressions of interesting genes in the brain (**Figure 7**). IDminer was built on open sources software systems, such as MongoDB database, Express web development framework, Nginx web server, and Ubuntu operating system. Python and R were used for data collection processing and analysis. A user-friendly web interface was provided to help users search and analyze the data online at http://<IP_ADDRESS>:3100/ IDminer/. The interface consists of seven parts: Home, Browser, Tools, Statistics, Download, Help, and Q&A. On the Home page, an introduction to the IDminer outlines a description statistic about all the data integrated into the database and the search module for gene and phenotype. There are two analysis tools for converting phenotypes and prioritizing candidate genes, respectively. Converting phenotypes is to help user mapping their clinical descriptions to our standard ID phenotypes, while coexpression analysis can be based on the brain gene expression data to study the expression profile of the interesting genes and its related genes. In the Document and Q&A pages, the guidelines for the database, and frequently asked questions and answers were showed. Furthermore, our database could be easily updated with the latest published information. For gene query, we provided basic gene information and linked it to multiple external databases, such as containing Ensemble, UniProtKB, GO, KEGG, and OMIM. Reported mutations, ID phenotypes, and patient information were also represented. Additionally, the gene's phenotypic pairs were also interactively visualized. When users entered a phenotypic item in the input box, we listed its basic information such as HPO ID, synonyms and phenotype definitions, reported patients with this phenotype, reported causative genes causing the patients, and its co-occurrence network. For reported genes, in addition to displaying detailed mutation information of these genes, we also annotated the genes' functions and performed PPI network analysis. Importantly, the query clinic feature could be enriched for some genes, and the

genes were also listed. Finally, the top co-occurred phenotypic pairs ranked by their frequencies were shown as a network and the enriched genes for each pair were shown by clicking the edge.

# A Use Case for the IDpred

The case of a real patient with the pathogenic gene AAR2 and the standardized phenotypes [Microcephaly (HP:0000252), Cochlear malformation (HP:0008554), Hypoplasia of the corpus callosum (HP:0002079), Ventricular septal defect (HP:0001629), Global developmental delay (HP:0001263), Anteriorly placed anus (HP:0001545), Macule (HP:0012733), Patent foramen ovale (HP:0001655)] was selected based on the previous studies (Charng et al., 2016). The other input candidate genes were randomly selected from the gene list in our database. Then, the query genes list consisted of MXRA8, DMBX1, AAR2, CLIC2, PLA2G6, and phenotypes list consisted of all the standardized phenotypes of this patient (genes and phenotypes are separated by semicolons) were entered into the corresponding box on the page of the website. Then the selection of the models (for example, SVM) with the appropriate parameters should be submitted (**Supplementary Figure S3A**). The result page contains seven columns (GeneSymbol, PathogenicGeneRank, PathogenicScore, Pathogenicity, SimilarRank, SimilarScore, and Phenotypes) would be displayed. On the result table, PathogenicGeneRank is the rank of the input pathogenic genes compared to all deposited genes in our database, PathogenicScore is the score of the pathogenic genes, and Pathogenicity is defined as "Probably" (PathogenicScore > 0.5) or "Less likely" (PathogenicScore = 0.5). SimilarRank refers to the rank of similarity between gene and phenotypes, and SimilarScore refers to the calculated score of similarity between gene and phenotypes. Phenotypes listed the phenotypes related to the GeneSymbol. As shown in the result of this case, AAR2 was predicted as the pathogenic gene with the highest pathogenic score of 0.788 (**Supplementary Figure S3B**).

# DISCUSSION

Our work manually extracted a large number of genes, clinic phenotypes and basic information of the patients from published ID literature. By integrating these data for comprehensive analysis, we have provided a holistic view of the current genetic research of ID and made the correlation of various clinic factors of ID patients, prompting researchers to further explore the mechanisms causing ID. The mutation spectrum delineated in our datasets provided essential information for molecular diagnosis in ID patients. Though most genes had its major mutation types, the spectrum showed that all mutation types were identified in ID cases. This combination of mutation types

raises the need of using several clinical detection methods for ID diagnoses such as Array Comparative Genomic Hybridization (aCGH), target panel sequencing, whole exon sequencing, and even whole genome sequencing (De Ligt et al., 2013; Redin et al., 2014). Notably, because a small number set of genes accounted for most ID patients, targeted panel sequencing may be favorable than other methods in consideration of cost, time and the difficulty of the data analysis.

The phenotypes of ID patients were extremely diverse and heterogeneous. Unlike the previous study of phenotypebased clustering (Kochinke et al., 2016), we mapped the phenotypes of ID patients to HPO items and found the 2075 phenotypes in total 3803 patients. We confirmed not only mutations in different genes could lead to various phenotypes, but defects in a single gene had been implicated in different phenotypes. Interestingly, there was also considerable phenotypic heterogeneity even among individuals who have identical mutations in the disease gene. We speculated that, besides various genes, the heterogeneity of phenotypes could be affected by other factors, such as mutation types, genetic background, and environment. Though the phenotypes of ID patients were heterogeneous, the specific phenotypes for genes could be analyzed and used for prioritizing caused genes. A previous report suggests that, for tubulinopathies, each mutated gene has

other families. (C) The oncoprint-like representation of specific phenotypes for genes. (D) The oncoprint-like representation of specific phenotypic pairs for genes.

number of given phenotypes.

an associated predominant pattern of cortical dysgenesis (Bahi-Buisson et al., 2014). Additionally, the previous studies in ID found that convergent molecular pathways result in common phenotypes (Kochinke et al., 2016), allowing some phenotypegenotype correlation. However, the common phenotypes for each gene could be achieved until recently the applications of NGS, aCGH, target sequencing, WES, and WGS to ID patients, which lead to an increase of diagnosis. This larger sample size could raise the power of the statistical significance test. Then, for some genes, a large number of patients are sufficient to statistically to find the specific phenotypes, phenotypic pairs and co-occurrence networks for the genes. These features were extracted with enrichment in patients subgroup caused by each gene, confirming the phenotype-genotype correlations and the convergence of ID phenotypes among their extreme heterogeneities.

With the deepening of ID research and the increase of reported patients, it also requires the development of analytical tools for ID researchers to understand the data. Therefore, providing online friendly and easy-to-use analysis tools will also greatly assist in the research of the entire ID field. So, our website not only provides

FIGURE 6 | Performance comparison of pathological gene prediction between IDpred and other algorithms. (A) ROC curve derived from IDpred model based on 10 fold cross validation. (B) the percentage of predicted pathogenic gene derived from IDpred, phenolyzer, and varElect. (C) cumulative distribution of TopN rate base on the rank of the pathogenic gene derived from IDpred, phenolyzer, and varElect.

a knowledge base of ID but also aggregates tools commonly used in ID analysis. And more analysis tools for ID will be added in the future to promote ID research as much as possible.

Overall, our data and analysis showed the convergences of ID genes and phenotypes among their extreme heterogeneities. For genes, the convergence was characterized by the fact that a small percentage of genes could explain the majority of ID phenotypes. And for phenotypes, it was represented as genes' specific phenotype and phenotypic pairs. Importantly, we provided analysis tools based on ID genes and phenotypes in hopes of establishing the standard ID gene and phenotype libraries and, in turn, aiding in clinical diagnosis. Overall, the findings and tools could contribute to the understanding of the genetic basis of ID disease and ultimately improve the diagnosis and treatment of the disease.

# CONCLUSION

Our analysis provided evidence to support, though the ID genes and phenotypes were extremely heterogeneous, the genetic bias and phenotypic convergence deserved our more attention, which may help to help us to quickly diagnose ID patients and further promote the studies of disease mechanisms. Moreover, our curated data, analysis, and developed tools were integrated to build a standard ID database IDminer, which could be accessed through http://<IP_ADDRESS>:3100/IDminer/. The database and interface are user-friendly for geneticists and clinicians, and a very wide range of ID researchers.

#### DATA AVAILABILITY STATEMENT

Publicly available datasets were analyzed in this study. This data can be found here: http://<IP_ADDRESS>:3100/IDminer/.

#### REFERENCES

# AUTHOR CONTRIBUTIONS

Z-CF, YW, L-NZ, and X-WM conceived the project, analyzed and interpreted the data, and wrote the manuscript. FY, X-LX, YY, XY, WP, W-QZ, J-YL, W-DZ, and X-LZ acquired data and performed bioinformatics analyses. T-JJ edited the manuscript. All authors read and approved the final manuscript.

#### FUNDING

This work was supported by the National Key Research and Development Program of China (Project Nos. 2016YFC1000707 and 2017YFC1001701), and the Sanming Project of Medicine in Shenzhen (SZSM201606088).

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe. 2020.00045/full#supplementary-material

FIGURE S1 | The flowchart of pathogenic gene prediction model specification. The framework for training and test data generation, model construction, and validation.

FIGURE S2 | Feature importance scores derived from IDpred. Feature importance is defined as the average gain of the feature in trees from XGBoost in IDpred\_XGBoost.top, IDpred\_XGBoost.random, and IDpred\_XGBoost.low model.

FIGURE S3 | Dpred interface and direct mode example. (A) The user enters three types of input: gene symbol list, phenotype expression and modeling type. (B) Output results presented in a tab with seven columns which were defined as GeneSymbol, PathogenicGeneRank, PathogenicScore, pathogenicity, SimilarRank, SimilarScore, and Phenotypes.

causes neurodegeneration and aberrant autophagy. Ann. Neurol. 83, 153–165. doi: 10.1002/ana.25130

**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.

Copyright © 2020 Wang, Zhu, Ma, Yang, Xu, Yang, Yang, Peng, Zhang, Liang, Zhu, Jiang, Zhang and Feng. 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.

# Integrative Analysis of Methylation and Gene Expression in Lung Adenocarcinoma and Squamous Cell Lung Carcinoma

Hao Zhang1† , Zhou Jin1,2† , Ling Cheng3 and Bin Zhang1 \*

1 Department of Respiratory and Critical Care Medicine, Second Affiliated Hospital of Zhejiang University School of Medicine, Hangzhou, China, 2 Department of Respiration, Hospital of Traditional Chinese Medicine of Zhenhai, Ningbo, China, 3 Shanghai Engineering Research Center of Pharmaceutical Translation, Shanghai, China

#### Edited by:

Tao Huang, Shanghai Institutes for Biological Sciences (CAS), China

#### Reviewed by:

Yunguang Tong, University of California, Los Angeles, United States Hong Liu, Temple University, United States

#### \*Correspondence:

Bin Zhang<EMAIL_ADDRESS><EMAIL_ADDRESS>
†These authors have contributed equally to this work

#### Specialty section:

This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Bioengineering and Biotechnology

> Received: 31 October 2019 Accepted: 03 January 2020 Published: 07 February 2020

#### Citation:

Zhang H, Jin Z, Cheng L and Zhang B (2020) Integrative Analysis of Methylation and Gene Expression in Lung Adenocarcinoma and Squamous Cell Lung Carcinoma. Front. Bioeng. Biotechnol. 8:3. doi: 10.3389/fbioe.2020.00003 Lung cancer is a highly prevalent type of cancer with a poor 5-year survival rate of about 4–17%. Eighty percent lung cancer belongs to non-small-cell lung cancer (NSCLC). For a long time, the treatment of NSCLC has been mostly guided by tumor stage, and there has been no significant difference between the therapy strategy of lung adenocarcinoma (LUAD) and squamous cell lung carcinoma (SCLC), the two major subtypes of NSCLC. In recent years, important molecular differences between LUAD and SCLC are increasingly identified, indicating that targeted therapy will be more and more histologically specific in the future. To investigate the LUAD and SCLC difference on multi-omics scale, we analyzed the methylation and gene expression data together. With the Boruta method to remove irrelevant features and the MCFS (Monte Carlo Feature Selection) method to identify the significantly important features, we identified 113 key methylation features and 23 key gene expression features. HNF1B and TP63 were found to be dysfunctional on both methylation and gene expression levels. The experimentally determined interaction network suggested that TP63 may play an important role in connecting methylation genes and expression genes. Many of the discovered signature genes have been supported by literature. Our results may provide directions of precision diagnosis and therapy of LUAD and SCLC.

Keywords: lung adenocarcinoma, squamous cell lung carcinoma, methylation, gene expression, Boruta, Monte Carlo Feature Selection

# INTRODUCTION

Lung cancer, considered to be a highly prevalent type of cancer, is a leading cause of cancer-related mortality worldwide, resulting in 1.6 million deaths each year with poor 5-year survival rate of about 4–17% (Hirsch et al., 2017; Altorki et al., 2019).   It can also be typically divided into three main subtypes, lung adenocarcinoma (LUAD), squamous cell lung carcinoma (SCLC), and large cell cancer (LCC), according to standard pathology methods

(Socinski et al., 2016; Swanton and Govindan, 2016; Herbst et al., 2018). Compared with squamous lung cancer, adenocarcinoma was associated with better prognosis. Despite the advances in diagnostic and therapeutic technology, lung cancer remains a serious global public health concern.

For a long time, the treatment of NSCLC has been mostly guided by tumor stage, and there has been no significant difference between the therapy strategy of LUAD and SCLC. Most lung cancers are usually diagnosed at an advanced stage and are treated primarily with systemic chemotherapy, typically with platinum-based regimens (Bishop et al., 2010). Recent progress in characterization of NSCLC by molecular typing, especially in adenocarcinomas of the lung, have brought new investigation of therapeutic agents that target dominant oncogenic mutations, such as epidermal growth factor receptor (EGFR)-targeted therapies, which have showed improved response rates in patients with NSCLC (Shigematsu et al., 2005).

Currently, progress in molecular biology of lung cancer has resulted in the identification of multiple potential biomarkers that may be related to the clinical management of NSCLC patients. In recent years, with the emergence of next-generation sequencing technologies, important molecular differences between LUAD and SCLC are increasingly identified, indicating that targeted therapy will be more and more histologically specific in the future (Kim et al., 2005; Sun et al., 2007; Li et al., 2014). Several studies have identified multiple gene expression subtypes that differ in prognosis, genomic alterations, clinical characteristics, including tumor differentiation, stage-specific survival, underlying drivers, and potential responses to treatment within LUAD and SCLC (Wilkerson et al., 2010; Thomas et al., 2014; Lu et al., 2016). For example, LUAD patients that harbor EGFR, ALK, ROS1, or BRAF mutations were discovered to benefit the most (Villalobos and Wistuba, 2017; Herbst et al., 2018). Targeted therapies for gene abnormalities of HER2, MET, RET, and NTRK1 appear to be an effective approach to treat LUAD (Dearden et al., 2013; Mazieres et al., 2013). SCLC shows different mutation spectrum from that of adenocarcinoma, and the mutation targeted therapy for SCLC has not been thoroughly studied to obtain approved treatment (Bunn et al., 2016; Soldera and Leighl, 2017).

A series of imaging studies suggested that NSCLC may progress rapidly between occurrence and primary treatment (Koh et al., 2017). Therefore, it is necessary for clinicians to identify between these two subtypes of NSCLC in a convenient and rapid way. With the improvement of the above clinical and molecular levels, growing evidences have shown that immunohistochemistry (IHC) is an effective tool for differentiating adenocarcinoma from squamous cell carcinoma (Bass et al., 2009; Weiss et al., 2010).

It is reported that the formation and development of lung cancer are related to the accumulation of permanent genetic changes and dynamic epigenetic changes. Therefore, enhancing our understanding of tumor biology and gene expression profiles will be critical for cancer treatment and diagnosis. In this study, an integrative analysis of lung cancer methylation data and gene expression data was performed, and mixed features were also screened out for analysis.

# MATERIALS AND METHODS

## The Joint Methylation and Expression Profiles of Lung Cancer Patients

The methylation and gene expression profiles of lung cancer patients were obtained from GEO (Gene Expression Omnibus)1 . The data were originally generated by Karlsson et al. (2014). They used the data to cluster the patients into five groups, and these groups showed different overall survival (Karlsson et al., 2014). We were more interested in how the methylation and expression differ from well-known subtypes, especially LUAD and SCLC. Therefore, we analyzed the 77 LUAD and 22 SCLC patients who had both methylation and expression data.

The methylation profiles were measured with Illumina HumanMethylation450 BeadChip while the gene expression profiles were measured with Illumina HumanHT-12 V4.0 expression BeadChip. The probe expression levels were averaged onto 20,178 genes. The 354,251 methylation sites within genes were analyzed. Therefore, each patient was represented with 20,178 genes and 354,251 methylation sites.

# Screen for the Relevant Methylation and Expression Features

Since the number of methylation and expression features was very large, it was difficult to analyze directly. We applied the Boruta method (Kursa and Rudnicki, 2010) to screen the combined data and identify the relevant methylation and expression features. The Boruta method was based on random forest classification, and the relevance of features to sample classes was measured by the ensemble of the random forest classifier's stochasticity.

## Evaluate the Importance of Relevant Methylation and Expression Features

After the irrelevant features were removed, the relevant methylation and expression features were ranked based on their importance evaluated with MCFS (Monte Carlo Feature Selection) (Draminski et al., 2008). The MCFS was a widely used method to rank features based on classification trees (Chen et al., 2018, 2019; Pan et al., 2018, 2019a,b; Li et al., 2019). First, for the d features, we selected s subsets and each subset included m features (m was much smaller than d). Then, for each subset, t trees were constructed. Based on the s × t trees, we can estimate a feature's importance by considering how many times it appeared in these trees and how well it performed in these trees as a node. By comparing the permutation results, the significance of features was evaluated.

## Perdition Performance of the Mixed Methylation and Expression Signature

The MCFS can find the significant top-ranking features by comparing with permutations. To objectively evaluate the significant top-ranking features' prediction performance, we performed LOOCV (Leave One Out Cross Validation) using

1https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE60645

SVM (Support Vector Machine) classifier (Li et al., 2018; Sun et al., 2018; Pan et al., 2019a). Each time, one sample was chosen as test samples and all other samples were used to train the SVM predictor. After all samples were tested once, we compared the actual sample classes with predicted sample classes and calculated the sensitivity, specificity, accuracy, and Mathew's correlation coefficient (MCC) based on the confusion matrix (Huang et al., 2011, 2013; Cai et al., 2012).

#### RESULTS AND DISCUSSION

#### Rank the Methylation and Expression Features

The methylation and gene expression data were combined and, therefore, each lung cancer patient was represented with mixed methylation and gene expression features. The number of mixed features (20,178 gene expression features and 354,251 methylation features) was too large to conduct sophisticated statistical analysis. So, we removed irrelevant features using the Boruta method (Kursa and Rudnicki, 2010). At last, 711 relevant features were remained.

Then, these 711 Boruta selected features were further ranked with the MCFS method (Draminski et al., 2008). As a classification tree-based ensemble learning algorithm, MCFS can rank the features based on how many times and how much it contributed to the sample classification in s × t trees. By comparing with permutation results, it can evaluate the significance of features.

#### Identify the Methylation and Expression Signature

The 136 significant top-ranking features were identified using the latest dmLab version 2.3.0 software downloaded from2

2https://home.ipipan.waw.pl/m.draminski/mcfs.html

TABLE 1 | The 136 methylation and gene expression signature identified with the MCFS method.

with default parameters. These 136 methylation and expression signatures are given in **Table 1**.

It can be seen that within these 136 signature features, there were 113 methylation features and 23 gene expression features. The annotations of the 113 methylation features based on GPL135343 are provided in **Supplementary Table S1**. We plotted the heatmaps of LUAD and SCLC lung cancer patients with 113 methylation features and 23 gene expression features in **Figures 1**, **2**, respectively. Both the 113 methylation features and 23 gene expression features can successfully group almost all samples with only three misclassified SCLC samples. They did not show difference on cluster results.

To more objectively and carefully compare the performance of the 113 methylation features and 23 gene expression features, we conducted LOOCV with SVM classifier. The LOOCV prediction performances of the 136 mixed features, 113 methylation features and 23 gene expression features are listed in **Tables 2**–**4**. It can be seen that the prediction results of 113 methylation features

3https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL13534

TABLE 2 | The confusion matrix using 136 mixed methylation and gene expression features.

TABLE 3 | The confusion matrix using 113 methylation features.

were the same as the 136 mixed features and better than the 23 gene expression features. The 23 gene expression features had one more misclassified SCLC samples. It seemed that methylation had better performance.

#### Comparison With CNV Signature

Comparing with the 136 LUAD and SQCLC CNV signatures identified by Li et al. (2014), we found that the methylated genes HORMAD2, KLHL3, LPP, and PTPN3 are also CNAs genes.  Lipoma preferred partner (LPP) may be an important candidate molecular marker for the classification of NSCLC tissue subtypes. PTPN3 can inhibit lung cancer by regulating EGFR signal (Li et al., 2015). However, there are no reports of KLHL3 in lung cancer, which also suggests a new idea of candidate molecular markers for the identification of lung cancer subtypes.

#### The Relationship Between Methylation and Expression Signature Genes

The 113 methylation features can be mapped onto 93 genes. We overlapped the selected methylation feature genes and

TABLE 4 | The confusion matrix using 23 gene expression features.

expression feature genes and found that HNF1B and TP63 were dysfunctional on both methylation and gene expression levels. HNF1B was one of the DNA methylated markers of the same subtype (Matsuo et al., 2014; Shi et al., 2017). TP63, also known as P63, was considered to be the most common marker for SCLC (Bishop et al., 2012; Van de Laar et al., 2014).

We downloaded the 66 lung cancer genes from KEGG hsa05223 NSCLC4 and mapped them and the overlapped two genes: HNF1B and TP63, onto STRING network (Szklarczyk et al., 2018). TP63 interacted with 39 KEGG lung cancer genes: AKT1, AKT3, ALK, BAK1, BAX, CASP9, CCND1, CDK4, CDK6, CDKN1A, CDKN2A, DDB2, E2F1, E2F2, E2F3, EGF, EGFR, EML4, ERBB2, FHIT, FOXO3, GADD45A, GRB2, HRAS, KRAS, MAP2K1, MAPK1, MAPK3, NRAS, PIK3CA, PIK3CB, PIK3R1, RB1, STAT3, STAT5A, STAT5B, STK4, TGFA, and TP53. HNF1B interacted with 14 KEGG lung cancer genes: AKT1, AKT2, CCND1, CDKN1A, CDKN2A, EGF, HRAS, KRAS, MAPK1, MAPK3, PIK3CA, RXRA, STAT3, and TP53.

What's more, we searched the methylation genes and expression genes in STRING database (Szklarczyk et al., 2018) and extracted the experimentally determined interaction and plotted the network in **Figure 3**. The light-yellow nodes were methylation genes, the light-blue nodes were expression genes. The overlapped methylation and expression genes were marked in red, the overlapped methylation and CNV genes from Li et al. (2014) were marked in pink. It can be seen that TP63 played an important role in connecting methylation genes and expression genes. The methylation genes and expression genes were closely connected to form a dense functional module on the network.

4https://www.genome.jp/dbget-bin/www\_bget?pathway+hsa05223

# The Biological Significance of the Identified Signature

To develop more specific and individualized targeted therapy, there is an urgent need to improve our knowledge on the molecular basis, in addition to different phenotypes. It is noteworthy that adenocarcinoma and squamous cell carcinoma show marked differences in expression profiles, DNA methylation, and lesion location. In this study, the features containing methylation and expression data were screened by Boruta and then further sorted by MCFS. After comparing the selected features with related literatures, a certain correlation was found between these features and lung cancer subtypes.

In this study, 113 methylation features were screened and mapped to 93 genes. We inquired about the functions of these genes and their relationship with lung cancer to discuss whether they have the potential as molecular markers to recognize LUAD and SQCLC. Many genes have been proved to promote or inhibit the progression of lung cancer. For instance, FOXK1

was expressed in many malignant tissues (Huang and Lee, 2004) and Ma et al. (2018) also found that FOXK1 plays a carcinogenic role in lung cancer. MAD1L1 is a checkpoint gene, with its mutation been proved to play a pathogenic role in lung cancer (Tsukasaki et al., 2001). Some genes have been reported to be related with the prognosis of NSCLC, such as HORMAD2 and ANO1. The overexpression of ANO1 is related to the high expression of EGFR, which can be used as a predictor of recurrence after NSCLC (He et al., 2017). In addition, according to Zhang et al. (2014) HORMAD2 gene polymorphism has great potential prognostic value in Chinese patients with NSCLC. Other genes are associated with NSCLC subtypes, such as another member of the FOX family, FOXK2, which was reported to be closely related to the overall survival of LUAD (Chen et al., 2017). DOK1 and HOPX were found to serve as lung tumor suppressors for LUAD (Berger et al., 2010; Chen et al., 2015). In the study of Zhou et al. (2017) the methylation locus of PARD3 gene was positively correlated with the expression of PARD3 and suppression of PARD3 intensified chemoresistance in LUAD cells. SFTA3 was found obviously overexpressed in LUAD, and its expression in LUAD and SQCLC was quite different. Therefore, the sensitivity and specificity of using SFTA3 to distinguish the two subtypes will be relatively high (Zhan et al., 2015). ARHGEF1 aliased p114RhoGEF and its expression might help to predict progression and survival of SQCLC patients (Song et al., 2013). Notably, LPP has multiple functions of actin binding protein and transcriptional coactivator (Kuriyama et al., 2016). Ngan et al. (2017) proved that the expression of LPP reduces the number of circulating tumor cells and inhibits lung cancer metastasis. Kang et al. (2009) used high-resolution array-CGH to find that the difference in genomic imbalance patterns between SQCLC and LUAD was most significant in 3q26.2-q29, while LPP (3q28) was significantly targeted in SQCLC, suggesting that LPP may be an attractive candidate molecular marker for histological subtype classification of NSCLC and may be involved in the pathogenesis of SQCLC.

We also investigated 23 expressed genes in lung cancer, and found that many studies clearly indicated that some genes were associated with LUAD or SQCLC. DSC3 (Han et al., 2014; Lv et al., 2015) and KRT5 (Xu et al., 2014; Travis et al., 2015) have been proved to be an effective marker of SQCLC. ANXA8 (Chao et al., 2006) and DSG3 (Savci-Heijink et al., 2009) were significantly over-expressed in SQCLC, and DSG3 could be an effective ancillary marker to identify SQCLC (Sanchez-Palencia et al., 2011; Gómez-Morales et al., 2013). VSNL1, also known as VILIP-1, was a tumor suppressor gene specific to SQCLC (Fu et al., 2008). KRT6A, KRT6B, and KRT6C, members of the keratin protein family, are specific to squamous cells and associated with epidermis of squamous epithelium (Fujii et al., 2002; Hawthorn et al., 2006; Chang et al., 2011). In addition, we also identified several genes primarily associated with LUAD. According to Balabko et al. (2014) RORC is a specific transcription factor in the tumor area of lung tissue in patients with LUAD. DLX5 (Kato et al., 2008; Balabko et al., 2014), MUC1 (Mashima et al., 2005; Molina-Pinelo et al., 2014), and

TABLE 5 | The GO enrichment results of the identified signature.

KRT17 (Erdogan et al., 2009; Liu et al., 2018) were found to be overexpressed in LUAD.

## The GO Enrichment Analysis of the Identified Signature

In order to further analyze the relationship between mixed characteristics and lung cancer, we carried out GO enrichment analysis. The results suggest that characteristic genes are mainly related to keratinization, epidermal cell differentiation, tissue development, and cytoplasm. The GO enriched results with FDR (False Discovery Rate) smaller than 0.05 are listed in **Table 5**. P63 appears to be useful in differentiating SQCLC from LUAD in small biopsies with no keratosis or glandular differentiation, helping to establish different treatments (Camilo et al., 2006). The expression of keratinocyte transglutaminase and cytokeratin 10 was measured as markers of squamous differentiation (Lokshin et al., 1999). Epidermal cell differentiation is related to EGFR signal pathway, which can inhibit the proliferation and metastasis of cancer cells, while EGFR mutation is largely limited to LUAD (Ladanyi and Pao, 2008). The expression of Promyelocytic leukemia zinc finger (PLZF) in SQCLC was weak or absent, which was significantly lower than that in LUAD (Xiao et al., 2015).

To sum up, most of the 113 methylated genes and 23 expressed genes we found are closely related to lung cancer, and some of them have the possibility of distinguishing SQCLC from LUAD, which is helpful for the targeted selection of lung cancer treatment and provide more research support for lung cancer molecular markers.

#### DATA AVAILABILITY STATEMENT

All datasets generated for this study are included in the article/**Supplementary Material**.

# AUTHOR CONTRIBUTIONS

fbioe-08-00003 February 6, 2020 Time: 16:33 # 8

HZ, ZJ, LC, and BZ contributed to the study design. HZ, ZJ, and LC conducted the literature search. HZ, ZJ, and BZ acquired the data. ZJ and LC wrote the manuscript. HZ and BZ performed the data analysis. All authors gave the final approval of the version to be submitted, read, and approved the final manuscript.

#### REFERENCES

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe. 2020.00003/full#supplementary-material

TABLE S1 | The annotations of the 113 methylation features.

senescence in human lung cancer. J. Pathol. 235, 397–407. doi: 10.1002/path. 4469

a new therapeutic target. J. Natl. Cancer Inst. 97, 765–777. doi: 10.1093/jnci/ dji133

coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 47, D607–D613. doi: 10.1093/nar/gky1131

**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.

Copyright © 2020 Zhang, Jin, Cheng and Zhang. 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.

# Analysis of the Clinicopathologic Characteristics of Lung Adenocarcinoma With CTNNB1 Mutation

Chao Zhou1,2, Wentao Li 2 , Jinchen Shao3 , Jikai Zhao3 and Chang Chen1\*

1 Department of Thoracic Surgery, Shanghai Pulmonary Hospital, Tongji University School of Medicine, Shanghai, China, 2 Department of Thoracic Surgery, Shanghai Chest Hospital, Shanghai Jiao Tong University, Shanghai, China, 3 Department of Pathology, Shanghai Chest Hospital, Shanghai Jiao Tong University, Shanghai, China

Introduction: Lung adenocarcinoma with CTNNB1 mutation is relatively uncommon, and its clinicopathologic characteristics, disease course, and prognosis have not been well-studied.

#### Edited by:

Tao Huang, Shanghai Institutes for Biological Sciences (CAS), China

#### Reviewed by:

Liang Liu, Fudan University Shanghai Cancer Center, China Deli Liu, Cornell University, United States

> \*Correspondence: Chang Chen<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Genetics

Received: 17 October 2019 Accepted: 12 December 2019 Published: 07 February 2020

#### Citation:

Zhou C, Li W, Shao J, Zhao J and Chen C (2020) Analysis of the Clinicopathologic Characteristics of Lung Adenocarcinoma With CTNNB1 Mutation. Front. Genet. 10:1367. doi: 10.3389/fgene.2019.01367 Methods: A total of 564 lung adenocarcinoma patients were enrolled in this study. The relationship between CTTNB1 mutational status and clinicopathologic parameters, the rates of relapse-free survival (RFS) and overall survival (OS), and the mutational status of other genes commonly mutated in lung adenocarcinoma were analyzed.

Results: Of 564 lung adenocarcinoma patients, 30 (5.3%) harbored CTNNB1 mutations. Univariate analyses revealed that gender, smoking history, pleural invasion, and histological subtype were all significant predictors of RFS and OS. Pleural invasion and histological subtype remained significant predictors of RFS and OS in a multivariate analysis. There were no significant differences in RFS (p = 0.504) or OS (p = 0.054) between lung adenocarcinoma patients with CTNNB1 mutation and those without CTNNB1 mutation. However, patients with CTNNB1 mutation tended to have a worse OS.

Conclusions: Female patients and nonsmokers are likely to harbor CTNNB1 mutation and primary lung adenocarcinoma with mutated CTNNB1 has a poor prognosis.

Keywords: adenocarcinoma, lung cancer, CTNNB1, mutation, prognosis

# INTRODUCTION

Lung cancer, which has the highest incidence of all cancers and the highest rate of disease-related fatalities, is the main cause of cancer-related death worldwide (Torre et al., 2015; Gu et al., 2017a; Gu et al., 2018). Lung adenocarcinoma is the most common pathological subtype, accounting for nearly 70% of all lung tumors (Sun et al., 2010). With the introduction of low-dose computed tomography, which enables earlier detection, the incidence of lung cancer, especially early-stage lung cancer, has risen sharply in recent years (Field et al., 2012).

**109**

With the advent of genomics, molecular or genetic variants affecting disease risk can be identified (Field, 2008).  Woenckhaus et al. (2008) identified a number of differentially expressed genes in smoke-exposed bronchial epithelium and nonsmall cell lung cancers (NSCLCs), they found in adenocarcinomas, the cytoplasmic expression of beta-catenin was associated with shorter survival (p = 0.012). Shigemitsu et al. (2001) found CTNNB1 is infrequently mutated in lung cancer. Akyol et al. (2019) defined an immunohistochemical algorithm to dissect Wnt pathway alterations in formalin-fixed and paraffin-embedded neoplastic tissues and found all six colon adenomas of the 126 total adenomas studied for the altered/ mutant b-catenin staining pattern had presumptively pathogenic point mutations or deletions in CTNNB1. The N-terminus of b-catenin, with contains conserved phosphorylated threonine/ serine amino acid residues, is the most frequent location of cancer-related CTNNB1 mutations (Dar et al., 2017). The level of free b-catenin in the cytoplasmic pool is regulated by ubiquitination and proteasomal degradation (Akyol et al., 2019). b-catenin is a member of the Wnt signaling cascade and is associated with cadherin-mediated cell–cell adhesion systems (Woenckhaus et al., 2008). In lung tumors, the immunohistologic loss of b-catenin membrane staining along with a corresponding increase cytoplasmic or nuclear staining has been reported (Nozawa et al., 2006).

Although CTNNB1 mutation occurs in many tumors types, it has not been well-studied in the context of lung adenocarcinoma, and the clinicopathologic characteristics and prognosis of lung adenocarcinoma with mutated CTNNB1 has not been described. Therefore, we compared the clinicopathologic characteristics of 30 lung adenocarcinomas with CTNNB1 mutations with those of 534 lung adenocarcinomas with wild-type CTNNB1.

# MATERIALS AND METHODS

From July 2008 to April 2013, resected primary lung adenocarcinomas were collected at the Department of Thoracic Surgery of Shanghai Chest Hospital, Shanghai Jiaotong University. To confirm the diagnosis of primary lung cancer, all the patients received thorough preoperative testing at our hospital, including physical exams, serological tests, pulmonary function tests, chest/brain computed tomography (CT), technetium bone scanning, and abdominal ultrasound. Biopsies were done by bronchoscopy or endobronchial ultrasound-guided transbronchial needle aspiration, and in some cases, positron emission tomography CT was used to exclude mediastinal lymph node metastases (Gu et al., 2017b). The lung adenocarcinoma subtype was determined by light microscopy intraoperatively, using frozen sections, and confirmed postoperatively, using paraffin-embedded sections. All surgical samples had at least 5% tumor content. Each case was reviewed by at least two junior pathologists and a senior pathologist to confirm the histologic subtype of resected lung neoplasms. The combination of routine preoperative examination and intra-/postoperative pathological diagnosis is recommended to make an exact lung cancer diagnosis.

In total, 601 patients with primary lung adenocarcinoma were identified. Of these, 17 and 20 patients were excluded because they received neoadjuvant chemotherapy or were lost to follow-up, respectively. The remaining 564 patients were enrolled in this study.

Informed consent was given by all patients or their legal representatives. The study was initiated after obtaining Institutional Review Board approval. The medical records for all patients were reviewed to collect corresponding clinicopathologic data, including sex, age, smoking status, pathologic tumor, node, and metastasis (TNM) stage [according to the staging system of the 7th edition of the American Joint Committee on Cancer (Edge et al., 2010)], thyroid transcription factor-1 status, and treatment information. Data on disease recurrence and survival were obtained from follow-up clinic visits or by telephone.

#### Bioinformatics Analysis

Data of The Cancer Genome Atlas (TCGA) were analyzed by Gene Expression Profiling interactive Analysis (http://gepia. cancer-pku.cn/) and Kaplan–Meier Plotter (http://kmplot.com/ analysis/index.php?p=service&cancer=lung). Gene CTNNB1 were further analyzed by Gene Expression Profiling interactive Analysis and the survival curves were draw and compared by Kaplan–Meier Plotter.

## Mutational Analysis

The mutational status of EGFR, KRAS, and CTNNB1 was determined by targeted sequencing and verified by DNA sequencing analysis. Relevant primers were designed to amplify all known ALK fusion variants by quantitative real-time reverse transcriptase PCR of cDNA. ALK fluorescent in situ hybridization was used to confirm the presence of ALK gene fusions (Wang et al., 2012).

# Statistical Analysis

Clinicopathologic data was analyzed using the SPSS 22.0 software package (SPSS Inc, Chicago, IL). Relapse-free survival (RFS) and overall survival (OS) were estimated by the Kaplan– Meier method, and differences were compared by log-rank testing using Prism 6 (GraphPad Software, La Jolla, CA). A p value of <0.05 was considered statistically significant.

# RESULTS

#### Mutational Status of Lung Adenocarcinomas

Of the 564 lung adenocarcinoma patients examined, 30 (5.3%) harbored CTNNB1 mutations (Table 1). The distributions of specific mutation types are shown in Figure 1.

#### Relationship Between Clinicopathologic Factors and CNTTB1 Mutational Status

Of the 30 patients with CTNNB1 mutations, there were 21 (70%) female patients and 9 (30%) male patients, ranging in age from 22 to 81 years (median, 59 years). Histologically, 263 of the tumors were acinar-predominant (47%), 115 were micropapillary-predominant (20%), 94 were papillarypredominant (17%), 49 were lepidic-predominant (9%), 32 were mucinous adenocarcinoma-predominant (5%), and 11 were solid-predominant (2%). Most of the patients had earlystage lung cancer in both the CNTTB1 mutation group (stage I: 18/30, 60%) and the CNTTB1 wild-type group (stage I: 241/534, 45%). Age (p = 0.851), tumor size (p = 0.256), lymph node status (p = 0.184), pathologic stage (p = 0.322), and the presence of pleural invasion (p = 0.459) were similar between lung adenocarcinomas with CTNNB1 mutation and lung adenocarcinomas without CTNNB1 mutation, but the former group tended to have more female patients (p < 0.001) and more smokers (p = 0.019) (Table 2).

## Relationship Between CNTTB1 Mutational Status and Survival

Univariate analysis revealed that gender, smoking history, pleural invasion, and histological subtype were all significant predictors of RFS and OS (Table 3). Pleural invasion and histological subtype were still significant predictors of RFS and OS in a multivariate analysis (Table 4).

During follow-up, 19 (63.3%) patients with lung adenocarcinomas with mutated CTNNB1 and 259 (48.5%) patients with lung adenocarcinomas with wild-type CTNNB1 experienced a relapse, and 10 (33.3%) and 111 (20.8%) patients died, respectively. There were no statistically significant differences in RFS (p = 0.504) or OS (p = 0.054) between patients with CTNNB1 mutation and patients without CTNNB1 mutation (Figure 2). However, patients with CTNNB1 mutation tended to have a worse OS.

Zhou et al. Lung Adenocarcinoma With CTNNB1 Mutation

TABLE 2 | Features of patients with lung adenocarcinoma harboring CTNNB1 mutations.

As for lung adenocarcinomas from TCGA, there was no significant differences in the distributions of CTNNB1 mRNA expression among different lung adenocarcinoma stages (Supplementary Figure 1). Besides, between lung adenocarcinoma patients with and without CTNNB1 mutation, there was no significant differences in RFS (p = 0.49), while significant differences were found in OS (p = 8.9e−05). (Supplementary Figures 2 and 3).

#### DISCUSSION

Lung cancer remains the leading cause of cancer-related death worldwide (Gu et al., 2017). Low-dose computed tomography TABLE 3 | Independent predictors of overall survival.

#### TABLE 4 | Independent predictors of relapse-free survival.

FIGURE 2 | Survival curves for relapse-free survival and overall survival according to CTNNB1 status. (A) Relapse-free survival between the two groups. (B) Overall survival between the two groups.

screening reduces the mortality of lung cancer by as much as 20% in high-risk patients (National Lung Screening Trial Research Team, 2011). Early detection and diagnosis increases the number of patients who are newly diagnosed with lung cancer while it is still early-stage, improving the prognosis of lung cancer patients as a whole. For individuals at a high risk of developing lung cancer, periodic screening could have a survival benefit.  With the development of gene mutation testing, targeted therapy has changed the treatment strategy for lung cancer. In this study, we describe the clinicopathological characteristics of lung adenocarcinoma with CTNNB1 mutation.

b-catenin is important for the establishment and maintenance of the epithelial layer and is a key downstream component of the canonical Wnt signaling pathway. The WNT/b-catenin pathway is involved in cancer and pluripotent stem cell signaling, which may suggest the mechanism underlying cancer stem cells. In this study, out of 564 patients, 30 (5.3%) patients with CTNNB1 mutations were identified. Kase et al. (2000) conducted an immunohistochemical analysis of 331 lung cancer specimens and reported that b-catenin expression was reduced in 122 (37%) of the samples, which was associated with significantly worse patient survival. Similarly, Woenckhaus et al. (2008)reported that reduced membrane staining of b-catenin and its abnormal accumulation in the cytoplasm and/or nuclei of lung adenocarcinoma cells was associated with shorter survival (p = 0.012). Another study also suggested that reduced b-catenin expression in surgically resected non-small cell lung cancer specimens was associated with lymph node metastasis and a poor prognosis (Retera et al., 1998). These studies suggest that decreased expression of b-catenin is associated with an unfavorable prognosis in lung cancer.

In our study, during follow-up, 19 patients (63.3%) with lung adenocarcinomas with CTNNB1 mutations and 259 patients (48.5%) with lung adenocarcinomas with wild-type CTNNB1 relapsed, and 10 (33.3%) and 111 (20.8%) patients died, respectively. Patients with CTNNB1 mutations therefore tended to have a worse prognosis, although this did not reach statistical significance. When compare with data from TCGA, patients with CTNNB1 mutation in TCGA also had worse OS. Our findings therefore correspond well to the results of previous studies and common directory (Sunaga et al., 2001).

In Cox proportional hazards models, univariate analyses revealed that gender, smoking history, the presence of pleural invasion, and histological subtype were all significant predictors of RFS and OS. Pleural invasion and histological subtype remained significant predictors of RFS and OS in a multivariate analysis. With respect to histological subtype, adenocarcinoma patients with micropapillary or solid subtypes, which are defined as highrisk subtypes in the 2011 classification proposed by the International Association for the Study of Lung Cancer/ American Thoracic Society/European Respiratory Society (Travis et al., 2011), had significantly worse prognosis. As for pleural invasion, pleural invasion, as well as visceral invasion, is considered an aggressive and invasive factor in NSCLC and has been included in the TNM staging system as a factor that should upstage the T factor (Rami-Porta et al., 2007; Travis et al., 2008; Butnor and Travis, 2012). Shimizu et al. (2005) demonstrated that velopharyngeal insufficiency (VPI) is a significant and independent predictor of a poor prognosis regardless of tumor size or N status, and as a result, VPI is a good indicator of the degree of invasion and aggressiveness of NSCLC. As more earlystaged lung neoplasms are detected, whether VPI has impact on survival of patients with early-staged lung cancer is unknown. Therefore, Jiang et al. (2015) published a meta-analysis and found VPI together with tumor size has a synergistic effect on survival in patients with N0 disease. Patients with stage IB NSCLC and larger tumor size with VPI might be considered for adjuvant chemotherapy after surgical resection and need careful preoperative evaluation and postoperative follow-up.

There are several limitations to this study. First, the sample size was relatively small. Contributing to the small sample size, there were several patients with CTNNB1 gene mutations who could not be included in the data analysis because of incomplete clinicopathological records. Finally, the patients' outcomes could have been influenced by the use of different treatment strategies, which may confound the survival analysis.

In summary, our results suggest that female patients and nonsmokers are likely to harbor CTNNB1 mutation and primary lung adenocarcinoma with mutated CTNNB1 has a poor prognosis. Further research is needed to verify our results. However, these data suggest that b-catenin could be a potential therapeutic target for advanced-stage lung cancer.

#### DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

#### ETHICS STATEMENT

The studies involving human participants were reviewed and approved by Ethics committee of Shanghai chest hospital. The

#### REFERENCES

patients/participants provided their written informed consent to participate in this study.

#### AUTHOR CONTRIBUTIONS

Conception and design: CZ and CC. Development of methodology: CZ, WL, and CC. Acquisition of data (provided surgical samples, gene detection, pathological diagnosis, acquired and managed patients, provided facilities, etc.): CZ, WL, JS, JZ, and CC. Analysis and interpretation of data: CZ. Writing, review, and/or revision of the manuscript: CZ and CC.

#### FUNDING

This work was supported by Shanghai Pulmonary Hospital Innovation Team Program (fkcx1906).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019. 01367/full#supplementary-material

lung cancer. J. Thoracic Oncol. 3 (12), 1384–1390. doi: 10.1097/JTO.0b013 e31818e0d9f

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.

Copyright © 2020 Zhou, Li, Shao, Zhao and Chen. 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.

# Transcriptome Changes of Mycobacterium marinum in the Process of Resuscitation From Hypoxia-Induced Dormancy

Jun Jiang1† , Chen Lin1† , Junli Zhang1 , Yuchen Wang1 , Lifang Shen1 , Kunpeng Yang1 , Wenxuan Xiao1 , Yao Li 1\*, Lu Zhang1,2\* and Jun Liu1,3\*

#### Edited by:

Tao Huang, Shanghai Institutes for Biological Sciences (CAS), China

#### Reviewed by:

Li Zhang, East China Normal University, China Yongyong Yang, Northwestern University, United States

#### \*Correspondence:

Yao Li<EMAIL_ADDRESS>Lu Zhang<EMAIL_ADDRESS>Jun Liu<EMAIL_ADDRESS>
† These authors have contributed equally to this work

#### Specialty section:

This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Genetics

Received: 14 October 2019 Accepted: 11 December 2019 Published: 07 February 2020

#### Citation:

Jiang J, Lin C, Zhang J, Wang Y, Shen L, Yang K, Xiao W, Li Y, Zhang L and Liu J (2020) Transcriptome Changes of Mycobacterium marinum in the Process of Resuscitation From Hypoxia-Induced Dormancy. Front. Genet. 10:1359. doi: 10.3389/fgene.2019.01359 1 State Key Laboratory of Genetic Engineering, Institute of Genetics, School of Life Science, Fudan University, Shanghai, China, 2 State Key Laboratory of Genetic Engineering, MOE Engineering Research Center of Gene Technology, School of Life Science, Fudan University, Shanghai, China, 3 Department of Molecular Genetics, University of Toronto, Toronto, ON, Canada

Nearly one-third of the world's population is latently infected with Mycobacterium tuberculosis (M. tb), which represents a huge disease reservoir for reactivation and a major obstacle for effective control of tuberculosis. During latent infection, M. tb is thought to enter nonreplicative dormant states by virtue of its response to hypoxia and nutrientdeprived conditions. Knowledge of the genetic programs used to facilitate entry into and exit from the nonreplicative dormant states remains incomplete. In this study, we examined the transcriptional changes of Mycobacterium marinum (M. marinum), a pathogenic mycobacterial species closely related to M. tb, at different stages of resuscitation from hypoxia-induced dormancy. RNA-seq analyses were performed on M. marinum cultures recovered at multiple time points after resuscitation. Differentially expressed genes (DEGs) at each time period were identified and analyzed. Co-expression networks of transcription factors and DEGs in each period were constructed. In addition, we performed a weighted gene co-expression network analysis (WGCNA) on all genes and obtained 12 distinct gene modules. Collectively, these data provided valuable insight into the transcriptome changes of M. marinum upon resuscitation as well as gene module function of the bacteria during active metabolism and growth.

Keywords: transcriptional regulation, resuscitation, M. marinum, hypoxia, latency

# INTRODUCTION

Mycobacterium tuberculosis (M. tb), the causative agent of tuberculosis (TB), is the leading cause of death due to an infectious disease globally, with an estimated 10 million new cases and 1.3 million deaths in 2017. There were an additional 300,000 deaths from TB among HIV-positive people. The success of M. tb as a leading pathogen is associated with its ability to infect and persist in the host. About 1.7 billion people, 23% of the world's population, are estimated to have a latent TB infection (LTBI), which is asymptomatic but can persist for decades (Stewart et al., 2003; North and Jung., 2004). About 5–10% of LTBI will eventually develop active disease, and host immunosuppression (e.g., HIV coinfection) markedly increases the risk of reactivation (Corbett et al., 2007). LTBI poses a

**116**

major challenge to the effective control of TB because of the difficulty in treatment and the fact that LTBI represents a huge disease reservoir.

During LTBI, M. tb is thought to enter nonreplicative 'dormant' states by virtue of its lowered or altered metabolism in response to hypoxia, nitrosative stress, and/or nutrient deprivation (Boshoff and Barry, 2005). Accordingly, much research has been focused on environmental conditions and genetic programs that induce bacteriostasis, and the most extensively studied culture condition is hypoxia (Wayne and Sohaskey, 2001; Rustad et al., 2009). It was shown that an immediate bacterial response (2 hr) was the coordinated upregulation of 47 M. tb genes under the control of the response regulator (DosR) and two sensor kinases (DosS and DosT), known as the DosR regulon (Sherman et al., 2001; Boon and Dick, 2002; Park et al., 2003). A second set of 230 genes, induced by longer hypoxia exposure (7 days), was also identified (Rustad et al., 2008). These genes, collectively known as the enduring hypoxic response (EHR), were DosR-independent genes (Rustad et al., 2008).

During the reactivation of LTBI, the dormant bacteria are believed to resuscitate and resume active growth and metabolism. A few recent studies have used reaeration of hypoxic cultures for in vitro modeling of reactivation or resuscitation (Veatch and Kaushal, 2018). Several regulatory proteins, such as transcription factor ClgR and sigma factors SigH and SigE, were found to play a role in M. tb resuscitation from hypoxia (Mcgillivray et al., 2015; Iona et al., 2016; Veatch et al., 2016).

Despite the progress, knowledge of the genetic programs used to facilitate entry into and exit from the nonreplicative dormant states remains incomplete. In this study, we examined the transcriptional changes of Mycobacterium marinum (M. marinum) at different stages of resuscitation from hypoxia-induced dormancy. M. marinum is a pathogenic Mycobacterium and the closest genetic relative of the M. tb complex. M. marinum is an excellent model through which to understand various aspects of host–pathogen interactions in M. tb pathogenesis. For example, M. marinum and M. tb share many virulence determinants, such as the ESX-1 secretion system (Tobin and Ramakrishnan, 2008) and lipid virulence factors phthiocerol dimycocerosates and phenolic glycolipids (Yu et al., 2012). As such, findings from the current study of M. marinum may be applicable to M. tb.

# RESULT

## RNA-Seq Analysis of M. marinum Recovered From Hypoxia

Larry Wayne and co-workers were the first to develop an in vitro model to mimic the hypoxic environment of the human granuloma (Wayne, 1977; Wayne and Hayes., 1996; Wayne and Sohaskey., 2001). In the Wayne model, a sealed, standing culture is incubated over an extended period while the bacteria deplete the available oxygen. The gradual depletion of oxygen leads to nonreplicating persistence states with a concomitant shift in gene expression and metabolism.

To gain insight into the genetic mechanisms that facilitate the exit of mycobacteria from the nonreplicative state, we grew M. marinum under hypoxia for 7 days using the Wayne model and then reaerated the cultures. At different time points thereafter (0, 0.5, 4, 12, 24, and 48 hr), M. marinum cultures were collected and subjected to RNA-seq analysis. The growth curve of the bacteria is shown in Supplementary Figure 1. A total of 18 samples were collected (three biological replicates at each time point) and analyzed.

The RNA-seq reads showed a high mapping ratio for all samples (>96%) (Table 1), supporting the overall sequencing accuracy. Transcripts of more than 4,900 genes were detected in each sample. We compared the RNA-seq data of cultures recovered at different time points under aerobic conditions. As expected, results showed that the correlation coefficient decreased as the interval between two samples increased (Figure 1). This result also suggested that the recovery from hypoxia is a gradual but dynamic process.

#### Dynamic Changes of Gene Expression at Different Stage of Resuscitation

To analyze transcriptome changes of M. marinum, we focused on genes with RPKM ≥ 10 and compared samples from adjacent intervals: between 0.5 and 0 hr, 4 and 0.5 hr, 12 and 4 hr, 24 and 12 hr, as well as 48 and 24 hr. Differentially expressed genes (DEGs) were identified, and this was defined as a fold change greater than 2 and false discovery rate P value less than 0.05.

At the earliest time point after resuscitation (0.5 hr), 136 DEGs were detected, of which 71 were upregulated and 65 were downregulated. Between 4 and 0.5 hr, most of the DGEs were downregulated (81 out of 88 total DEGs). The numbers of DGEs found between 12 and 4 hr, 24 and 12 hr, as well as 48 and 24 hr

TABLE 1 | RNA-seq reads of 18 samples (three biological replicates at each time point).

were 72, 85, and 172, respectively. The heat map of the five groups of DEGs is shown in Figure 2.

We performed a Venn analysis of the five DEG groups (Figure 3A). The proportion of unique genes in each group was high: 68.4% (93/136, between 0.5 and 0 hr), 60.2% (53/88, between 4 and 0.5 hr), 46.7% (35/72, between 12 and 4 hr), 47.1% (40/85, between 24 and 12 hr), and 74.4% (128/172, between 48 and 24 hr). This suggests that a variety of genes were involved at different stages of resuscitation.

Combing the DEGs of different stages, a total of 440 genes were identified (Supplementary Table 1), and their expression underwent dynamic changes during resuscitation (Figure 3B). For example, the expression of MMAR\_0922, MMAR\_3562, and MMAR\_1654 were significantly changed at the early stage of resuscitation, suggesting that they may play an important role in this period. Some genes had changed in multiple time periods. For example, the expression of MMAR\_3403 was changed in last three periods, suggesting that this gene may be associated with the late stage of resuscitation.

#### Validation of RNA-seq Results by RT-qPCR

To validate the RNA-seq results, a real-time quantitative (RTqPCR) analysis was performed. Three biologically independent samples at each time point were used for this experiment. For each DEG group, we selected the top 10 upregulated genes and 10 downregulated genes for this analysis (Figure 4). Between 0 and 0.5 hr, MMAR\_4852 and MMAR\_5170 (whiB4) were significantly downregulated, and six genes, MMAR\_5122 (lipX), MMAR\_0548 (espG3), MMAR\_0547 (esxR), MMAR\_0551 (eccE3), MMAR\_0546 (esxG), and MMAR\_0550 (mycP3), were significantly upregulated (Figure 4A). Between 0.5 to 4 hr, six genes, including MMAR\_1656, MMAR\_1658 (hycQ), MMAR\_1653 (Rv0081), MMAR\_1655, MMAR\_5122 (lipX), and MMAR\_5170 (whiB4), were significantly downregulated and four genes MMAR\_0845 (hemB), MMAR\_5484, MMAR\_1908 (ATC1), and MMAR\_3776 (rpfE) were significantly upregulated (Figure 4B). Between 4 and 12 hr, MMAR\_2343 (papA1), MMAR\_3555, and MMAR\_2320 (wecE) were significantly upregulated and MMAR\_4903 were significantly downregulated (Figure 4C). Between 12 and 24 hr, six genes, MMAR\_0335, MMAR\_0602, MMAR\_4903, MMAR\_4899, MMAR\_2009, and MMAR\_0615 (iniA), were significantly downregulated (Figure 4D). Between 24 and 48 hr, six genes, MMAR\_3465 (PPE51), MMAR\_4824, MMAR\_4482 (cypM), MMAR\_4750, MMAR\_2944, and MMAR\_1790 (PPE2), were significantly downregulated, and seven genes, MMAR\_2651, MMAR\_2914 (katG), MMAR\_2649, MMAR\_5315 (lpqH), MMAR\_5319, MMAR\_2839 (mpt63), and MMAR\_0656, were significantly upregulated (Figure 4E).

There is a good agreement between the RNA-seq and qPCR data, evident from the scatter plot using the expression levels of all 97 genes that were analyzed by both RNA-seq and qPCR (R2 = 0.784) (Figure 4F). Based on this result, we consider that the RNA-seq data is reliable.

number of DEGs in each period is shown on the right.

#### Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Enrichment Analysis

(F) Scatter plot of RT-PCR data of all genes analyzed in (A–E), comparing them to RNA-seq data of the same genes.

To gain insight into the biological consequence of the observed transcriptome changes, we performed GO and KEGG pathway analyses of the five DEG groups (Figure 5). A GO analysis was applied to identify the functional categories of DEGs. Between 0.5 and 0 hr, more than half of DEGs were involved in membrane and cell wall processes and were significantly enriched (Figure 5). This is consistent with the notion that, upon reaeration, the bacteria resumed cell division, which involved cell wall and membrane biogenesis. A KEGG analysis consistently revealed that DEGs involved in metabolic pathways and biosynthesis were significantly enriched and accounted for the largest number. A similar trend was observed for later periods, in which genes involved in cell wall and membrane biogenesis were highly enriched and accounted for the majority of DEGs at these stages. These results provide snapshots into the recovery of the bacteria from hypoxia and active growth under aerobic conditions.

#### Co-Expression Analysis Between Transcriptional Regulators and mRNAs

Co-expression networks can show relationships between genes. To explore the regulatory mechanisms at different stages of resuscitation, we constructed a co-expression network between transcriptional regulators and mRNAs. For this, we selected known transcriptional regulators, including transcription factors and sigma factors from the five DEGs, and calculated the correlation coefficients between these transcription factors and the remaining DEGs in the same group. We considered that a relationship existed between a given transcriptional regulator and other genes if the absolute value of the correlation coefficient was greater than 0.9, which included both positive and negative correlation. Based on these results, we constructed five

co-expression modules and integrated them into a large network (Figure 6). The dark blue nodes in the figure represent transcription factors, and the green nodes denote DEGs in the same period. The size of the node is determined by the degree of connectivity. Greater degrees of connectivity are indicated by larger points. If there is a line between two nodes, then there is a relationship between them.

In the first co-expression module (between 0 and 0.5 hr), three transcription factors, MMAR\_4874 (CosR), MMAR\_1653 (Rv0081), and MMAR\_4852 (KmtR), formed the major regulatory hubs, and MMAR\_4874 (CosR) was the largest hub and interacted with other hubs in the network (Figure 6). The MMAR\_4874 (CosR) and MMAR\_1653 (Rv0081) hubs remained in the second co-expression module (between 0.5

FIGURE 6 | Co-expression networks in five periods. Each dot represents a gene. Transcription factors are labeled in blue, and other genes are labeled in green. Line between dots represents co-expression relationship between genes. The size of dot is proportional to the level of connectivity.

and 4 hr), in addition to three new hubs formed by MMAR\_4254, MMAR\_1725 and MMAR\_1132. In the third period (between 4 and 12 hr), MMAR\_0229 and MMAR\_4902 formed the hubs. In the fourth period (between 12 and 24 hr), MMAR\_2003 (SigB) and MMAR\_4219 formed the hubs. In the final period (24 to 48 hr), MMAR\_2651, MMAR\_1555, and MMAR\_0249 formed the hubs.

#### Weighted Gene Co-Expression Network Analysis (WGCNA)

The DEG analysis focused on partial dynamic changes during resuscitation. While the co-expression network of transcription factors provides an overview of the regulatory programs enabling the resuscitation of M. marinum, our knowledge on the overall genetic changes is still missing. Therefore, in this section, we analyzed the expression of all genes from cultures at different stages of resuscitation.

Weighted gene co-expression network analysis (WGCNA) (Langfelder and Horvath, 2008) is a method for analyzing the gene expression patterns of multiple samples. It clusters genes into modules by similar expression trends and reveals the relationship between gene modules and specific traits or phenotypes. We applied this method to analyze the RNA-seq data of M. marinum at different stages of resuscitation.

The gene modules were identified by the WGCNA package in R software. We first determined the appropriate "soft-thresholding" value, which emphasizes strong gene–gene correlations at the expense of weak correlations. An optimal parameter (power = 20) was determined by plotting the strength of correlation against a series (range 2 to 20) of soft threshold powers (Figure 7A).

An unsigned pairwise correlation matrix was calculated, and the WGCNA algorithm was used to transfer the correlation coefficient between genes into the adjacent coefficient. Then, the dissimilarity of the topological overlap matrix was calculated based on the adjacent coefficient. Using the calculated dissimilarity, we carried out a hierarchical analysis by using agglomerative hierarchical clustering, also known as the bottomup method. Other assumptions were made: (i) distances between

module fusion.

different classes were measured by the average connectivity, and (ii) there should be at least 30 genes in each gene module.

Based on these analyses, we initially obtained 48 gene modules. The hierarchical cluster treewas then treated using the dynamic tree cut algorithm in the WGCNA package. A total of 13 gene modules were obtained. The "gray" module was the default module, which included discarded genes that could not be clustered. Thus, we focused on the analysis in the remaining 12 gene modules. The process offusion is shown in Figure 7B. The number of genes varied among these 12 modules, and the detailed information is listed in Table 2 and Supplementary Table 2.

The first principal component analysis (PCA) was performed on the 12 gene modules (Figure 8). The PCA results reflected the main trend of gene expression in the modules. Module 20 played an important role in the early stage (0 to 0.5 hr) of recovery, module 35 played a role mostly in the middle stage (4 to 12 hr), and module 12 was only involved in the last stage (24 to 48 hr). Other modules appeared to play roles in more extended periods.

#### Identification of Key Gene Modules Associated With Different Stages of Resuscitation

Thus far, we have identified 5 DEGs and 12 gene modules. We then performed an enrichment analysis between them. When the P value of Fisher's exact test was less than 0.001, we considered that these gene modules were significantly enriched in the DEG sets. The results are shown in Figure 9. Interestingly, module 20 was significantly enriched in DEGs of the early stage (0 to 0.5 hr) and module 12 was significantly enriched in DEGs of the last stage (24 to 48 hr) of the recovery. This is consistent with the result that these two modules were only involved in the early and last stages of resuscitation, respectively (Figure 8).

#### DISCUSSION

In this study, we examined the transcriptome changes of M. marinum recovered from hypoxia-induced dormancy. To gain a comprehensive view, multiple time points, including shortly after resuscitation (0.5 hr) to more extended periods up to 48 hr, were included. For each time point, three biologically independent

samples were analyzed. Transcripts of the whole genome were analyzed by RNA-seq, and the quality of the RNA-seq data was reflected by the high genome mapping ratio andfurther validated by qPCR analysis of close to 100 genes. With these high-quality sequence data, we performed in-depth analyses, which included the identification of DEGs and the construction of co-expression network of transcription factors in each period. The availability of transcriptomes of independent samples at multiple time points also allowed us to employ a weighted gene co-expression network analysis to identify gene modules of M. marinum. Collectively, these data provide valuable insight into not only the genetic changes of the bacteria upon resuscitation but also the genemodule function of M. marinum during active metabolism and growth.

A total of 136 DEGs were identified in M. marinum upon resuscitation from dormancy (0 to 0.5 hr), including eight transcription factors (Figure 6). Interestingly, all of these transcription factors were significantly downregulated. Among them, MMAR\_1653 is a homolog of Rv0081 in M. tb. Previous studies have shown that Rv0081 is a member of the DosR regulon and is induced at the early stage of hypoxia (Sherman et al., 2001). Rv0081 is a major regulator of M. tb response to hypoxia and forms a large regulatory hub (Galagan et al., 2013; Chawla et al., 2018; Sun et al., 2018). Rv0081 plays an important role connecting the early and enduring hypoxic responses (Sun et al., 2018). WhiB4 (MMAR\_5170) is an oxygen-sensitive transcription factor and has been shown to regulate PE/PPE family proteins, and it plays a role in M. marinum virulence (Chawla et al., 2012; Wu et al., 2017). CosR (MMAR\_4874) is a copper-inducible transcriptional regulator, and the loss of cosR resulted in a hypoxia-type response with the induction of the DosR regulon (Talaat et al., 2004; Ward et al., 2008; Marcus et al., 2016). Given their roles in the hypoxic response, it is not surprising that Rv0081, WhiB4, and CosR underwent dynamic changes in expression upon resuscitation by reaeration. Two other transcription factors that were downregulated at this stage, MMAR\_5405 (EthR) and MMAR\_1394 (Rv3176c), belong to the TetR family transcription factors (Leiba et al., 2014; Sharma et al., 2017).

 These proteins are components of the Type VII secretion systems, and many of them are important T cell antigens and play a critical role modulating the host–pathogen interactions (Abdallah et al., 2007).

From 0.5 to 4 hr after reaeration, 11 transcription factors were downregulated (Figure 6), including Rv0081, CsoR, and WhiB4 as mentioned above. Notably, the expression of two other WhiB family proteins (Averina et al., 2012), WhiB3 and WhiB5, were also significantly altered. WhiB3 responds to dormancy signals, including hypoxia and NO, and controls redox homeostasis of the bacteria (Priscille et al., 2004). WhiB5 responds to oxygen and controls the expression type VII secretion systems (Priscille et al., 2004). Consistently, we found that whiB3 was downregulated while whiB5 was upregulated at this stage of resuscitation.

FIGURE 9 | Enrichment analysis between 12 gene modules and 5 DEGs. The above bar chart represents the number of DEGs in each period, and the bar chart on the right represents the number of genes in each gene module. Blue squares represent a significant enrichment between the row set and the column set.

During the period of 4 to 12 hr, we identified 72 DEGs, including 7 PPE family genes (MMAR\_5121, MMAR\_1095, MMAR\_1235, MMAR\_1847, MMAR\_1905, MMAR\_2669, and MMAR\_3989). The PE/PPE family proteins play a critical role in mycobacterial pathogenesis (Fishbein et al., 2015). In addition, several genes of the Mce family were upregulated, including MMAR\_3865 (mce2B), MMAR\_3868 (mce5E), MMAR\_3867 (mce5D), MMAR\_3866 (mce5C). The Mce family genes comprise four mammalian cell invasion factor (mce) operons (mce1-4), and some of these are involved in the invasion of host cells (Zhang and Xie, 2011).

A picture appears to have emerged from these analyses; in the early stage of resuscitation from the hypoxia-induced dormancy, transcription factors critical for a hypoxia-induced response are downregulated, and, as the recovery continues, genes important for virulence and host interactions are upregulated.

A WGCNA analysis revealed 12 distinct gene modules. Of particular interest is gene module 20, which was involved in the early stage of resuscitation only (Figure 8). This module comprises of ~300 genes, many of which have unknown functions or annotations. Future studies focusing on genes in this module may help to understand the molecular machinery enabling the exit of the bacteria from dormancy.

#### METHODS AND MATERIALS

#### Bacterial Strain, Media, and Growth Conditions

M. marinum 1218R (ATCC 927) was grown in Middlebrook 7H9 broth to OD600~0.5, at which point they were aliquoted and cultured in screw-capped conical flasks at 30°C without additional oxygen. The hypoxic culture conditions were described previously by Wayne and Hayes (Wayne and Hayes, 1996). After 7 days in hypoxic conditions, the screw cap was replaced with a permeable membrane, and the rest of the conditions were unchanged. After aeration, samples were taken at 0 h, 0.5 h, 4 h, 12 h, 24 h, and 48 h, and an aliquot was used to measure the OD value (Supplementary Figure 1). The remaining samples were collected and snap frozen in liquid nitrogen for RNA sequencing.

#### RNA Extraction, Illumina Sequencing, and RT-qPCR

M. marinum cultures were centrifuged at 4,500 × g for 5 min at room temperature and frozen on dry ice. The frozen cell pellets were resuspended in 1 mL TRIzol reagent (CW Bio). RNA extraction and illumina sequencing were performed as previously described (Wu et al., 2017). Raw data of RNA sequencing have been uploaded to the GEO database (BioProject ID : PRJNA588556).

For RT-qPCR validation of RNA-seq data, 1 μg RNA was reversed-transcribed to cDNA, which was then used as the template for RT-PCR analysis. The primers for analyzing the selected genes were listed in Supplementary Table 3.

#### Transcriptome and Bioinformatics Analysis

The RNA-seq analysis and identification of differentially expressed genes (DEGs) were performed as previously described (Lee et al., 2019).

#### GO and KEGG Analysis