Circular RNAs (circRNAs) are a novel class of endogenous noncoding RNAs that have wellconserved sequences. Emerging evidence has shown that circRNAs can be novel biomarkers or therapeutic targets for many diseases and play an important role in the development of various pathological conditions. Therefore, identifying potential disease-related circRNAs is helpful in improving the efficiency of finding therapeutic targets for diseases. Here, we propose a computational model (PreCDA) to predict potential circRNA–disease associations. First, we calculated the circRNA expression similarity based on circRNA expression profiles. The circRNA functional similarity is calculated based on cosine similarity, and the disease similarity is used as the dimension of each circRNA vector. The associations between circRNAs and diseases are defined based on the circRNA functional similarity and expression similarity. We constructed a disease-related circRNA association network and used a graph-based recommendation algorithm (PersonalRank) to sort candidate disease-related circRNAs. As a result, PreCDA has an average area under the receiver operating characteristic curve value of 78.15% in predicting candidate disease-related circRNAs. In addition, we discuss the factors that affect the performance of this method and find some unknown circRNAs related to diseases, with several common diseases used as case studies. These results show that PreCDA has good performance in predicting potential circRNA–disease associations and is helpful for the diagnosis and treatment of human diseases.

#### *Edited by:*

*Lei Deng, Central South University, China*

#### *Reviewed by:*

*Leyi Wei, The University of Tokyo, Japan Hui Ding, University of Electronic Science and Technology of China, China*

*\*Correspondence:*

*Tianyi Zang<EMAIL_ADDRESS>Yadong Wang<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Statistical Genetics and Methodology, a section of the journal Frontiers in Genetics*

> *Received: 24 June 2019 Accepted: 13 August 2019 Published: 12 September 2019*

#### *Citation:*

*Wang Y, Nie C, Zang T and Wang Y (2019) Predicting circRNA-Disease Associations Based on circRNA Expression Similarity and Functional Similarity. Front. Genet. 10:832. doi: 10.3389/fgene.2019.00832*

Keywords: circRNA, disease, circRNA expression similarity, circRNA functional similarity, PersonalRank

## INTRODUCTION

Circular RNAs (circRNAs) are a type of RNA molecule that forms a covalently closed continuous loop from exon circularization (Motieghader et al., 2017; Xu, 2017). In recent years, advances in high-throughput sequencing technology have greatly facilitated the study of circRNAs (Jeck and Sharpless, 2014). When compared to other ncRNAs (Danan et al., 2012), circRNAs are highly stable. Circular RNAs have evolutionarily conserved sequence features across species, tissues, and developmental stages (Jens, 2013; Conn et al., 2015; Rybak-Wolf et al., 2015). Therefore, circRNAs have become hotspots in transcriptomics research.

Recent studies have shown that alterations in the expression of circRNAs play important roles in human disease and other biological processes (Xu, 2017; Zhao and Shen, 2017; Xia et al., 2018). For example, the best-known circRNA, CDR1as, as the inhibitor of miR-7, is a critical ncRNA known to be involved in cancer, neurodegenerative diseases, diabetes, and atherosclerosis (Li et al., 2015; Xu et al., 2018).

**40**

Researchers found that the circRNA ciRS-7 may be a promising target for neurodegenerative disorder (Lukiw, 2013) and myocardial infarction (Lin et al., 2018). The circRNA CircCCDC66 has been demonstrated to regulate colon cancer growth and metastasis as a miRNA sponge (Hsiao et al., 2017). The circRNA hsa\_circ\_0001895 is involved in the expression of cancer-related proteins in gastric cancer (Shao et al., 2017). The circRNA CircHIPK3 plays an important role in cell growth by sponging multiple miRNAs (Zheng et al., 2016). Moreover, circRNAs can be found in exosomes, cell-free saliva, and plasma (Li Y et al., 2015). Circular RNAs are emerging as novel biomarkers or therapeutic targets for many diseases due to their conservation, cell type–specific expression, and tissue-specific expression, and they play roles in the development of various pathological conditions (Meng et al., 2017; Vo et al., 2018).

Although a large number of circRNAs have been discovered, the mechanisms of circRNAs in many diseases remain unclear (Xu et al., 2018). To enable research on circRNAs and diseases, several databases have been constructed, such as circRNADisease (Zhao et al., 2018), CircR2Disease (Fan et al., 2018), and Circ2Disease (Yao et al., 2018). They provide important data support for circRNA– disease association analyses. Some methods have been proposed to provide the most promising disease-related biomarkers, including those involving lncRNAs (Chen et al., 2015; Gu et al., 2017; Cheng et al., 2018a; Cheng et al., 2019), miRNAs (Peng et al., 2019b; Shao et al., 2018), genes (Cheng et al., 2016; Hu et al., 2019; Peng et al., 2019a), and drugs (Jiang et al., 2017; Zhang et al., 2017), for further experimental validation. These methods can decrease the time and cost of biological experiments. However, very few methods have been developed to predict potential circRNA– disease associations (Lei et al., 2018), and both disease functional similarity and semantic similarity were not considered in these methods. Improved knowledge has suggested that exploring both the semantic and functional associations of diseases, which are two types of significant associations, is beneficial in measuring disease similarity (Cheng et al., 2014; Peng et al., 2018).

In this study, we proposed a computational model (PreCDA) for potential disease-related circRNA identification. In view of the limited number of circRNA–disease associations, we introduced disease similarity to solve possible sparse problems and built a disease-related circRNA similarity network. However, relying entirely on circRNA-related diseases greatly limits the utility of the method because many circRNAs still have very few or no associated diseases. To overcome this limitation, we calculated the circRNA expression similarity based on the existing data resources. Subsequently, we built a new disease-associated circRNA network by fusing circRNA functional associations and expression similarities. To assess the practicability and accuracy of this method, we designed a validation process with different datasets of circRNA–disease associations, as good computational models must perform well on different data sources. Finally, PreCDA proved successful in predicting potential disease-related circRNAs.

#### MATERIALS AND METHODS

#### Workflow

A flowchart of the PreCDA workflow is shown in **Figure 1**. We preprocessed circRNA and disease data because of the lack of

uniform identification of circRNAs and diseases. We extracted the synonym vocabulary from the two circRNA databases, including circRNADisease (Zhao et al., 2018) and circBase (Glažar et al., 2014). Then, we unified different representations of the same circRNA in different databases. Additionally, the identification of the Human Disease Ontology (DO) (Kibbe et al., 2015) was used as the unified marker of diseases in the computational model. We measured the similarity between circRNAs in two ways, including the circRNA expression similarity and functional similarity. We extracted circRNA expression profiles from circBase (Glažar et al., 2014) and CIRCpedia (Dong et al., 2018). The circRNA expression similarity was calculated based on the Spearman correlation coefficient. The disease similarity was used as the dimension of each circRNA vector, and the circRNA functional similarity was calculated based on cosine similarity. A disease-related circRNA association network was built based on the circRNA expression similarity and functional similarity. Finally, we identified potential candidate disease-related circRNAs based on the PersonalRank algorithm (PR) (Haveliwala, 2002).

## Data Preprocessing

#### circRNA Data

In this study, we used three circRNA databases for experiments and validations. The circRNADisease database is a manually curated database of experimentally supported circRNA and disease associations, which collected 330 circRNAs and 48 diseases in 354 associations. Each entry in the circRNADisease database includes detailed information on a circRNA–disease association, including the circRNA and disease name, the circRNA expression pattern, literature references, and other annotation information. CircR2Disease is a database for experimentally supported circRNA–disease associations and provides a platform for investigating the mechanism of diseaserelated circRNAs. The present version of CircR2Disease collected 661 circRNAs and 100 diseases. Circ2Disease is a database that curates experimentally supported human circRNAs and provides comprehensive associations between circRNAs and human diseases. It contains 273 manually curated associations between 237 circRNAs and 54 human diseases from 120 studies. However, currently, the naming of circRNAs has not yet been unified (Xu et al., 2018), which leads to the underutilization of information from different public circRNA databases. Therefore, we designed and collected mappings among different circRNA names provided by different circRNA databases, including circRNADisease and circBase. circRNADisease contains circRNA synonyms, and circBase is a database that merged and unified datasets of circRNAs. We mapped circRNAs from the three circRNA databases to circBase referring to circRNA synonyms. Then, we used circRNA IDs from circBase as the unified IDs of circRNAs in this work.

#### Disease Data

Human Disease Ontology represents common and rare human disease concepts captured across biomedical resources. Each node in DO represents one disease term and is organized in a directed acyclic graph with the relationship of "is\_a". MEDIC (Davis et al., 2012) integrates OMIM (Online Mendelian Inheritance in Man) terms (Amberger et al., 2015), synonyms and identifiers with MeSH terms (Lipscomb, 2000), synonyms, definitions, identifiers, and hierarchical relationships.

We extracted disease terms and synonyms from MEDIC to annotate DO by the same external references in DO and MEDIC, as shown in **Figure 2**. If a disease term was recorded in both DO and MEDIC, the term and its synonyms in MEDIC were used to annotate DO. With this approach, a given disease name can be matched to DO to a great extent by string matching, considering that the naming rules for diseases in different disease-related circRNA databases are different. The diseases described by different names are considered to be the same disease that has a unique id in DO if these disease names can match the disease term or its extended synonyms in DO.

#### circRNA Expression Similarity

Considering that comprehensive circRNA expression data are still unavailable, we extracted circRNA expression profiles from circBase and CIRCpedia, including the expression profiles of 92488 circRNAs in 78 human cell types or tissues. We used the Spearman correlation coefficient between the expression profiles of each circRNA as the circRNA expression similarity, as shown in Formula 1.

$$\rho = 1 - \frac{\sigma \sum d\_i^2}{n\left(n^2 - 1\right)}\tag{1}$$

where *di* is the difference between the two ranks of the expression scores in the *i*th human cell type or tissue, and *n* is the number of the human cell types or tissues from circBase or CIRCpedia. Matrix *CB* and Matrix *CP* are, respectively, denoted as the circRNA expression similarity matrix of circBase and CIRCpedia, where *CB(i,j)* and *CP(i,j)* are the expression similarities between circRNA *c(i)* and *c(j)*. Then, to obtain reliable performance for circRNA expression data, we defined the expression similarity between circRNA *c(i)* and *c(j)* as shown in Formula 2 if circRNA *c(i)* and *c(j)* are included in both circBase and CIRCpedia.

$$\text{ExcSim}\{i, j\} = \begin{cases} \text{Max}\left(\text{CB}\{i, j\}, \text{CP}\{i, j\}\right) & \text{Max}\left(\text{CB}\{i, j\}, \text{CP}\{i, j\}\right) \succeq \pi \\ 0 & \text{otherwise} \end{cases} \tag{2}$$

To reduce the impact of data noise, we set a threshold τ to filter out those weak similarities between circRNAs. The threshold τ is set to 0.7 based on our experiments.

#### circRNA Functional Similarity

We extracted circRNA–disease associations from these above circRNA databases and defined a relational matrix of circRNAs and diseases. For each circRNA, all diseases in the matrix can be used to make a vector in a multidimensional space. Because of the limited number of available disease–circRNA pairs, there is a data sparsity problem in the matrix. Therefore, we calculated the circRNArelated disease similarity and filled this matrix with predicted association scores based on disease–circRNA associations and the disease similarity. Here, we use FNSemSim (Wang et al., 2017) to calculate disease similarity. This method, which combines disease functional similarity and semantic similarity, has good performance for calculating similarities between diseases. The workflow of calculating circRNA functional similarity is shown in **Figure 3**.

To calculate the association between one circRNA and any disease, the similarities between this disease and all diseases that are directly related to this circRNA are calculated by FNSemSim. *C* is defined as the set of disease-related circRNAs, and *D* represents the set of circRNA-related diseases. DisSet(c) is defined as the set of diseases directly related to circRNA *c*. The association score between disease *dis* and circRNA *c* is defined as follows:

$$\text{Discore}\{\text{dis},\mathcal{c}\} = \begin{cases} \text{Max}\left(\text{FNSemSim}\{\text{dis},\text{dis}\_{i}\}\right) & \text{dis}\_{i} \in \text{DisSet}\{\mathcal{c}\}, \quad \text{dis} \notin \text{DisSet}\{\mathcal{c}\} \\\\ 1 & \text{dis} \in \text{DisSet}\{\mathcal{c}\} \end{cases} \tag{3}$$

where DisSet(c) ⊆D, 1*≤*i *≤* |DisSet(c)|; |DisSet(c)| is denoted as the number of diseases in DisSet(c). If this disease belongs to DisSet(c), the score is 1; otherwise, the score is defined as the maximum of similarities between this disease and all the diseases related to

circRNA *c*. Therefore, circRNA *c* can be depicted by a vector that is composed of circRNA-related diseases in a multidimensional space. We can calculate the functional similarity between any two circRNAs based on cosine similarity. The functional similarity between circRNA *c(m)* and *c(n)* is defined as follows:

$$FnSim\left(m,n\right) = \frac{\sum\_{i=1}^{\|\mu\|} \text{Score}\left(\text{dis}\_i, \mathcal{L}\left(m\right)\right) \times \text{Score}\left(\text{dis}\_i, \mathcal{L}\left(n\right)\right)}{\sqrt{\sum\_{i=1}^{\|\mu\|} \text{Score}\left(\text{dis}\_i, \mathcal{L}\left(m\right)\right)^2} \sqrt{\sum\_{i=1}^{\|\mu\|} \text{Score}\left(\text{dis}\_i, \mathcal{L}\left(n\right)\right)^2}}\tag{4}$$

where |D| represents the size of the circRNA-related disease set *D*, and dis*i* is the *i*th disease in the circRNA-related disease set *D*.

#### Prediction of Candidate Disease-Related circRNAs

We take circRNA functional similarity and expression similarity as weights to construct a circRNA association network. In this network, the weight between circRNA *c(i)* and *c(j)* is defined as shown in Formula 5. If *ExSim(i,j)* is greater than 0, the weight between circRNA *c(i)* and *c(j)* is the average value of their

functional similarity and expression similarity; otherwise, the weight is defined as the functional similarity between them.

$$\text{Circ}\,\text{Weight}\left(i,j\right) = \begin{cases} \left(\text{FnSim}\left(i,j\right) + \text{ExSim}\left(i,j\right)\right) / 2 & \text{if } \text{ExSim}\left(i,j\right) > 0\\ \quad \text{FnSim}\left(i,j\right) & \text{otherwise} \end{cases} \tag{5}$$

To predict candidate disease-related circRNAs, the associations between diseases and circRNAs are also considered in this network. The weight between circRNA *c* and disease *dis* is defined as shown in Formula 6. If the disease is directly related to circRNA *c*, the weight between them is 1; otherwise, the weight is 0.

$$\text{CircDis Weight} \left(i, j\right) = \begin{cases} \quad 1 & \text{if } \text{dis} \in \text{DisSet} \left(\mathcal{c}\right) \\ \quad 0 & \text{otherwise} \end{cases} \tag{6}$$

In this network composed of circRNAs and diseases, we identify novel candidate disease-related circRNAs based on the PR. PersonalRank algorithm, as a recommendation algorithm based on random walking, can reveal more information between a target node and all the others in a specific network. PersonalRank algorithm is defined as follows:

$$PR\left(i\right) = \left(1 - d\right)r\_i + d\sum\_{j \in \text{in}\left(i\right)} \frac{PR\left(j\right)}{\left|\text{out}\left(j\right)\right|}\tag{7}$$

where PR(*i*) represents the possibility value that node *i* is accessed; *d* is the transfer probability; out(*j*) represents the outdegree of node *j*; in(*i*) is the in-degree of node *i*; and *ri* is defined as follows:

$$r\_i = \begin{cases} \quad 1 & \text{if } i = t \\ \quad 0 & \text{if } i \neq t \end{cases} \tag{8}$$

where *t* represents the target node. According to previous studies (Kang et al., 2014; Cheng et al., 2018b), *d* is set to 0.85. The target node *t* in the network randomly moves to adjacent nodes with the probabilities of the edges between these nodes. After enough iterations, the probabilities from the target node to all the other nodes will become stable. Eventually, the algorithm outputs the relevance degrees between all the nodes and this target node.

#### RESULTS

#### circRNAs and Diseases

We calculated similarities between 323 circRNAs from circBase and CIRCpedia based on circRNA expression profiles. Then, we obtained 11,281 circRNA pairs based on the preset threshold. Additionally, we found 507 relationships between 58 diseases and 445 circRNAs by mapping DO terms to the diseases in CircR2Disease. We matched 26 diseases based on DO terms and extracted 293 relationships between 277 circRNAs and these diseases from circRNADisease. In Circ2Disease, 218 relationships between 37 diseases and 199 circRNAs were found. Based on DO terms and the unification of circRNA naming, we analyzed the three circRNA databases and found the same circRNAs and diseases among these databases, as shown in **Figure 4**. This provided the test data for the performance evaluation of PreCDA.

We separately calculated the similarities between 445 circRNAs from CircR2Disease, 277 circRNAs from circRNADisease and 199 circRNAs from Circ2Disease. Three circRNA association networks were built that in turn contained 96,580 associations

TABLE 1 | Information on the three circRNA association networks.

between 440 circRNAs associated with 56 diseases; 38,226 associations between 277 circRNAs associated with 26 diseases; and 18,915 associations between 195 circRNAs associated with 36 diseases. The detailed statistics of the circRNAs and diseases are shown in **Table 1**.

#### Performance

We designed a test scheme to assess the performance of PreCDA. First, we selected two circRNA–disease databases, one to build the circRNA association network and the other to provide test data. Then, we extracted the same diseases from the circRNA association network and the reference database. For a given disease, if any circRNA related to this disease in the reference database exists in the network, but the association between the circRNA and the disease does not, the circRNA can be used as a test case for the disease to assess the performance of this circRNA association network. The test scheme is shown in **Figure 5**.

In this article, we used three circRNA–disease databases, including CircR2Disease, circRNADisease, and Circ2Disease. For example, both circRNA hsa\_circ\_0000284 and liver cancer (DOID: 3571) were recorded in Circ2Disease and CircR2Disease. The circRNA hsa\_circ\_0000284 was related to liver cancer (DOID: 3571) in Circ2Disease but not in CircR2Disease. Therefore, we built a circRNA association network based on CircR2Disease and calculated the relevance degrees between liver cancer and all circRNAs unrelated to the disease. We calculated the area under the receiver operating characteristic curve (AUC) according to the ranking of the circRNA hsa\_circ\_0000284 among these circRNAs to measure the prediction results. To validate the reliability of the computational model, we conducted nine validation experiments based on this scheme involving 18 diseases. We built three circRNA association networks based on the three different circRNA–disease databases. The three data sources were also used as the reference data. Additionally, we merged the known circRNA–disease associations in the three databases as an additional control data source.

PreCDA had an average AUC value of 78.15% in predicting candidate disease-related circRNAs. Furthermore, it had an outstanding performance on some diseases. For example, diabetes mellitus (DOID: 9351) in the network from Circ2Disease had an AUC of 98.48% based on the control data from circRNADisease and an AUC of 93.04% based on the control data from CircR2Disease. Based on the control data from Circ2Disease, the AUC of osteoarthritis (DOID: 8398) was 97.44% in the network from CircR2Disease and 98% in the network from circRNADisease. In the network from Circ2Disease, the AUC of stomach cancer (DOID: 10534) was 56.41% based on the control data from circRNADisease; it had an AUC of 73.88% in CircR2Disease. This shows that the networks from the different data sources have different results for a disease based on the same control database. However, the AUCs in the other validation experiments achieved more than 65%. Even so, the performance of PreCDA is excellent in predicting candidate disease-related circRNAs. The performance of PreCDA based on the different databases and the different control data sources is shown in **Figure 6**.

#### Case Study

To further evaluate the performance of PreCDA in predicting potential disease-related circRNAs, we conducted some case studies, including prostate cancer (DOID: 10283), liver cancer (DOID: 3571), breast carcinoma (DOID: 3459), Alzheimer disease, and pancreatic cancer (DOID: 1793). We integrated the known associations between circRNAs and diseases in the three databases and prioritized candidate disease-related circRNAs based on PreCDA.

In the ranking of candidate circRNAs related to liver cancer (DOID: 3571), hsa\_circ\_0001727 (Qiu et al., 2018) ranked 4th, hsa\_circ\_0001946 (Yu et al., 2016) ranked 7th, and hsa\_ circ\_0001141 (Guo et al., 2017) ranked 19th. They ranked in the

top 3% and were associated with liver cancer. For prostate cancer (DOID: 10283), hsa\_circ\_0001946 (Zhang et al., 2018) and hsa\_ circ\_0001649 (Yi et al., 2016) ranked 3rd and 5th in the ranking, respectively. They were documented to be related to prostate cancer. For pancreatic cancer (DOID: 1793), CircRNA\_100782 (Chen et al., 2017), which ranked 1st in the ranking, was

validated to regulate pancreatic carcinoma proliferation through the IL6-STAT3 pathway. We found that some candidate circRNAs related to these diseases were included by Circ2Traits (Ghosal et al., 2013), which is a comprehensive database for circRNAs potentially associated with disease and traits. For example, hsa\_circ\_0000118, which ranked 1st in the ranking of candidate

circRNAs associated with prostate cancer, was documented to be potentially related to this disease in Circ2Traits. The prediction results of the case studies are presented in **Table 2**.

#### DISCUSSION

Although functional associations between circRNAs are measured based on circRNA expression profiles, there are many weak connections among them. To reduce the impact of data noise, we set a threshold to filter out those weak connections between circRNAs. Based on the above validation strategy and different thresholds, we conducted nine groups of experiments in which these three databases were used as a reference to each other and to test the performance of PreCDA. As shown in

**Figure 7**, the average AUC of PreCDA varied with the change in the threshold, and the computational model worked best when the threshold was set to 0.7.

We calculated circRNA similarities by only cosine similarity and built a circRNA association network. Additionally, we merged the known circRNA–disease associations in these three databases

TABLE 3 | Performance differences of predicting circRNA–disease pairs based on different data sources.

as an additional control data source. Based on the validation strategy mentioned above, we used these three databases to test the performance of the network. As shown in **Figure 8**, the average AUC was 77.22%, the minimum AUC was 69.26%, and the maximum AUC was 88.85%. In comparison, PreCDA has a more stable performance, with an average AUC of 78.15%. Its minimum and maximum AUCs are 71.83% and 95.72%, respectively.

We found that the performance of predicting potential disease–circRNA pairs in the disease-related circRNA association network was impacted by different data sources. The result of predicting the associations between the same diseases and circRNAs was different based on the different data sources that were used to build networks. For example, referring to CircR2Disease, some of the data to be tested in the networks built based on circRNADisease and Circ2Disease were the same. However, the AUC values of predicting the associations between them were different. As shown in **Table 3**, we predicted the associations between colorectal cancer (DOID: 9256) and four circRNAs, including hsa\_circ\_0001649, hsa\_circ\_0000284, hsa\_circ\_0014717, and hsa\_circ\_0001141. The AUC value for the network of circRNADisease was 71.86%. The performance of identifying the associations between colorectal cancer and these four circRNAs based on Circ2Disease was improved, and its AUC achieved 82.17%.

### CONCLUSIONS

Circular RNA plays an important role in the development of various pathological conditions. Research on circRNA is invaluable in explaining the underlying pathogenesis. Therefore, we proposed a computational model to identify candidate disease-related circRNAs. First, we calculated the circRNA expression similarity with the circRNA expression profiles. Then, the disease similarity was used as dimensions of circRNA vectors, and the circRNA functional similarity was calculated based on cosine similarity. We defined the associations between circRNAs and diseases based on the circRNA expression similarity and functional similarity. A disease-related circRNA association network was built, and potential candidate diseaserelated circRNAs were ranked by the PR.

### REFERENCES

We evaluated the performance of PreCDA with the help of data differences among these three databases, including CircR2 Disease, circRNADisease, and Circ2Disease. The results showed that the average AUC of PreCDA was 78.15%, and it had good performance in predicting potential disease-related circRNA signatures. We discussed the selection of the threshold and the impact of different data sources on the performance of PreCDA. Then, we used several common diseases as case studies and found some unknown circRNAs that could be related to these diseases based on PreCDA. The findings of this study could be further applied in analyzing diseases in a system biology perspective (Cheng and Hu, 2018) and helpful for researchers to improve disease diagnostics and treatments.

### DATA AVAILABILITY

PreCDA is implemented using a combination of Java and scala, and it is freely available from the website at https://github.com/ wythit/PreCDA.

### AUTHOR CONTRIBUTIONS

YoW and CN did data collection and preprocessing. And with the guidance of TZ and YaW, YoW finished the algorithm design and validation. YoW was the major contributor in writing the manuscript. All authors have read and approved the final version of the manuscript.

### FUNDING

Publication costs were funded by the National Key Research and Development Program of China (grant no. 2016YFC0901605, 2016YFC1201702-01) and the National High-tech R&D Program of China (grant no. 2012AA02A604, 2015AA020108).

### ACKNOWLEDGMENTS

TZ and YaW are he corresponding authors. We thank them for their guidance. We also thank Ling Wang and Zhenxing Wang for their valuable suggestions on our work. We are grateful to Rongjie Wang and Yanshuo Chu for their helpful assistance in writing.

Cheng, L., Li, J., Ju, P., Peng, J. J., and Wang, Y. D. (2014). SemFunSim: A New Method for Measuring Disease Similarity by Integrating Semantic and Gene Functional Association. *PLoS One* 9(6), 11. doi: 10.1371/journal.pone.0099415

human and mouse. *Nucleic Acids Res.* 47(D1), D140-D144. doi: 10.1093/nar/ gky1051

**Conflict of Interest Statement:** 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 © 2019 Wang, Nie, Zang and Wang. 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.*

# iRO-PsekGCC: Identify DNA Replication Origins Based on Pseudo k-Tuple GC Composition

*Bin Liu1,2\*†, Shengyu Chen3†, Ke Yan4 and Fan Weng4*

*1 School of Computer Science and Technology, Beijing Institute of Technology, Beijing, China, 2 Advanced Research Institute of Multidisciplinary Science, Beijing Institute of Technology, Beijing, China, 3 School of Informatics, Computing, and Engineering, Indiana University Bloomington, Bloomington, IN, United States, 4 School of Computer Science and Technology, Harbin Institute of Technology, Shenzhen, China*

Summary: Identification of replication origins is playing a key role in understanding the mechanism of DNA replication. This task is of great significance in DNA sequence analysis. Because of its importance, some computational approaches have been introduced. Among these predictors, the iRO-3wPseKNC predictor is the first discriminative method that is able to correctly identify the entire replication origins. For further improving its predictive performance, we proposed the Pseudo k-tuple GC Composition (PsekGCC) approach to capture the "GC asymmetry bias" of yeast species by considering both the GC skew and the sequence order effects of *k*-tuple GC Composition (*k*-GCC) in this study. Based on PseKGCC, we proposed a new predictor called iRO-PsekGCC to identify the DNA replication origins. Rigorous jackknife test on two yeast species benchmark datasets (*Saccharomyces cerevisiae*, *Pichia pastoris*) indicated that iRO-PsekGCC outperformed iRO-3wPseKNC. It can be anticipated that iRO-PsekGCC will be a useful tool for DNA replication origin identification.

#### *Edited by:*

*Liang Cheng, Harbin Medical University, China*

#### *Reviewed by:*

*Xun Lan, Tsinghua University, China Qiwen Dong, East China Normal University, China*

#### *\*Correspondence:*

*Bin Liu<EMAIL_ADDRESS>
*†These authors have contributed equally to this work and share first authorship*

#### *Specialty section:*

*This article was submitted to Statistical Genetics and Methodology, a section of the journal Frontiers in Genetics*

> *Received: 05 July 2019 Accepted: 13 August 2019 Published: 18 September 2019*

#### *Citation:*

*Liu B, Chen S, Yan K and Weng F (2019) iRO-PsekGCC: Identify DNA Replication Origins Based on Pseudo k-Tuple GC Composition. Front. Genet. 10:842. doi: 10.3389/fgene.2019.00842*

Availability and implementation: The web-server for the iRO-PsekGCC predictor was established, and it can be accessed at http://bliulab.net/iRO-PsekGCC/.

Keywords: replication origin identification, pseudo *k*-tuple GC composition, random forest, web-server, DNA sequence analysis

### INTRODUCTION

In the process of the cell cycle, DNA replication is one of the most important steps (Shirahige et al., 1998). Since the DNA replication is initiated from a specific region, which is called replication origin, identifying the DNA replication origin is especially important for studying drug developments, cell life activities, genetic engineering, etc. (Méchali, 2010). Experimental methods detect the replication origins by using Chromatin immunoprecipitation (Chip) with high cost (Lubelsky et al., 2012). Therefore, researchers are seeking computational methods to efficiently predict the replication origins only based on the sequence information. Compared with non-replication origins, replication origins show uneven distribution of G (guanine) and C (cytosine) (Lobry, 1996), and the concept of "GC Skew" (Grigoriev, 1998) was proposed. Later, some computational methods incorporated these characteristics into the predictors based on the replication origins (Zhang and Zhang, 1991; Zhang and Zhang, 1994; Grigoriev, 1998; Roten et al., 2002; Thomas et al., 2007; Gao and Zhang, 2008; Luo et al., 2014; Bu et al., 2018). In order to further improve the predictive performance, the discriminative

**51**

methods were proposed by using both the information of the positive and negative samples (Chen et al., 2012; Li et al., 2015; Zhang et al., 2016), and all of these methods mentioned above achieved the-state-of-the-art performance. A recent method iRO-3wPseKNC incorporated the "GC asymmetry bias" (Lobry, 1996; Grigoriev, 1998; Lubelsky et al., 2012; Li et al., 2014) into the prediction by representing the entire replication origins based on three-window-based PseKNC (3wPseKNC) (Liu et al., 2018b). Feature extraction methods are the keys for the performance improvement. In this regard, many features have been proposed, which can be easily generated by some software tools.

These existing computational methods have significantly enhanced the development of this hot area, but they all suffer from certain disadvantages or limitations, for example, as discussed above the GC Skew is an important feature of replication origins, but all the existing discriminative methods failed to directly use GC Skew to construct the predictors. Furthermore, the existing feature extraction methods cannot reflect the uneven distribution of G and C. To solve these problems, we followed the framework of iRO-3wPseKNC (Liu et al., 2018b), and proposed an improved predictor called iRO-PsekGCC for replication origin identification. iRO-PsekGCC cannot only capture the CG asymmetry bias by using *k*-tuple GC composition (or *k*-GCC), but can also incorporate the GC Skew into the concept of PseKNC (Chen et al., 2014a; Chen et al., 2014b).

#### MANUSCRIPT FORMATTING

#### Benchmark Datasets

In order to evaluate the performance of the proposed method, two recently established benchmark datasets of the *Saccharomyces cerevisiae* and *Pichia pastoris* (Liu et al., 2018b) were employed in this study, because they showed clear CG asymmetry distributions, which can be represented as:

$$\mathbb{E}\mathbb{S}^{\tau} = \mathbb{S}^{+}\_{\tau} \bigcup \mathbb{S}^{-}\_{\tau}, \quad \mathbb{T} = \begin{cases} & \text{1 for } \text{Saccharowycz } c \text{2} \\ & \text{2 for } \text{Pichia } \text{pastoris} \end{cases} \tag{1}$$

where the symbol ∪ represents the union, and − + represents the positive dataset containing 340 replication origins, and 1 − represents the negative dataset containing 342 non-replication origins; 305 replication origins are in positive dataset 2 + , and 302 non-replication origins are in the negative dataset 2 − . For both of the two benchmark datasets, the redundant samples have been removed by using CD-HIT software tool (Li and Godzik, 2006) with the most stringent cut-off threshold (80%).

#### Pseudo *k*-Tuple GC Composition (PsekGCC)

One of the key steps for constructing machine-learning predictors for analyzing biological sequences is feature extraction. Following the framework of three-window-based PseKNC (3wPseKNC) (Liu et al., 2018b), we proposed a feature extraction method called "Pseudo k-tuple GC composition (PseKGCC)" to directly incorporate the CG asymmetry bias (Lobry, 1996; Grigoriev, 1998; Lubelsky et al., 2012; Li et al., 2014) and GC skew (Grigoriev, 1998) into the predictor. In the following sections, we will introduce how to represent DNA samples by using PseKGCC.

A DNA sequence D can be formulated as follow:

$$\mathbf{D} = \mathbf{N}\_1 \mathbf{N}\_2 \mathbf{N}\_3 \cdots \mathbf{N}\_i \cdots \mathbf{N}\_L \quad \text{( $i = 1, 2, 3 \cdots, L$ )}\tag{2}$$

where *L* denotes the length of **D**, and

$$N\_i \in \{\text{A(adending),}, \text{G(cyto sine),}, \text{G(guanine)},$$

$$\text{T(thymine)}, \quad (i = 1, 2, 3, \dots, L) \tag{3}$$

which represents the *i*-th nucleobase in the sequence, and fi ∈ denotes the "member of '" in the set theory. Following the study (Liu et al., 2018b), **D** is divided into three windows by two parameters ε and δ, including front window, middle window, and rear window respectively. ε and 1 − δ denote the percentage of total nucleobases of **D** in the front window and rear window, respectively. The front window, middle window and rear window can be represented as **D**[1,η], D[η + 1,ξ], and **D**[ξ + 1, *L*], respectively, where η and ξ are defined as (Liu et al., 2018b),

$$\begin{cases} \begin{array}{c} \mathfrak{n} = \operatorname{Int}^{\mathbb{C}}[L \times \mathfrak{c}] \\ \mathfrak{k} = \operatorname{Int}^{\mathbb{C}}[L \times \mathfrak{d}] \end{array} , \end{cases}, \qquad \text{(\$0 < \mathfrak{c} < \mathfrak{d} < 1.0)}$$

where the symbol IntC represents the ceiling operator, which means to return the smallest integer value greater than or equal to the float number .

According to (Liu et al., 2018b), if D is formulated by the *k*-tuple nucleotide (or *k*-mer) (Liu et al., 2019b; Liu, 2017) based on the three windows strategy, it can be represented as follow:

$$\mathbf{D} = \left[ f\_1^{(1)} \cdots f\_\nu^{(1)} \cdots f\_{\mathbf{4}^k}^{(1)} f\_{\mathbf{4}^k+1}^{(2)} \cdots f\_{\mathbf{4}^k+\nu}^{(2)} \cdots f\_{\mathbf{4}^k+\nu}^{(2)} \cdots \right]$$

$$f\_{2 \times 4^k}^{(2)} f\_{2 \times 4^k+1}^{(3)} \cdots f\_{2 \times 4^k+\nu}^{(3)} \cdots f\_{3 \times 4^k}^{(3)} \right]^\mathbf{T} \tag{5}$$

where in vector operations, symbol 'T' denotes the transformation symbol, and in the sample D, the normalized frequency values of the corresponding *k*-tuple nucleotides appearing in the front window, middle window and rear window are represented as *f* (1), *f* (2), *f* (3), respectively. The feature vector's dimension is 3 × 4*k* .

This strategy was proposed to capture the patterns of "GC asymmetry bias" in yeast species genomes, and it is able to improve the predictive performance for identifying replication origins among multiple yeast species genomes. However, this approach has the following disadvantages: 1) the three windows strategy can only capture the local GC asymmetry bias of replication origins, but it cannot incorporate the GC asymmetry bias in a global fashion; 2) for large *k* values of *k*-tuple nucleotide, the dimension of the resulting feature vectors is high, which will cause high dimension disaster.

In order to overcome these disadvantages, we proposed a new composition of DNA sequence called "*k*-tuple GC composition (or *k*-GCC)" to capture the GC preference in the replication origins and their global interactions. *k*-GCC treats A (adenine) and T (thymine) as one nucleotide type represented as \*. Therefore, the alphabet of *k*-GCC is

$$\mathbf{N}\_{i} \in \{ \text{G(guanine)}, \text{C(cytosine)}, \text{\*} \}, \quad \text{(\$i = 1, 2, 3, \dots\$, \$L)} \tag{6}$$

Therefore, by replacing the *k*-tuple by k-GCC, a DNA sequence D can be represented as:

$$\mathbf{D} = \left[ f\_1^{(1)} \cdots f\_{\nu}^{(1)} \cdots f\_{\mathbf{j}^k}^{(1)} f\_{\mathbf{j}^k+1}^{(2)} \cdots f\_{\mathbf{j}^k+\nu}^{(2)} \cdots f\_{\mathbf{2}\times\mathbf{3}^k}^{(2)} f\_{\mathbf{2}\times\mathbf{3}^k+1}^{(3)} \cdots \right] \tag{7}$$

$$f\_{\mathbf{2}\times\mathbf{3}^k+1}^{(3)} \cdots f\_{\mathbf{2}\times\mathbf{3}^k+\nu}^{(3)} \cdots f\_{\mathbf{3}\times\mathbf{3}^k}^{(3)} \right]^\mathbf{T} \tag{7}$$

Compared with Equation 5, the *k*-GCC can efficiently reduce the dimension of the feature vector from 3 × 4*k* to 3 × 3*k* by focusing on the GC composition.

The proposed Pse-KGCC incorporates both the *k*-GCC and GC skew into the framework of PseKNC (Chen et al., 2014a), which can be represented as:

$$\mathbf{D} = \begin{bmatrix} \phi\_1 \cdots \phi\_{\mathbf{3}^k} \cdots \phi\_{\mathbf{3}^k + \lambda} \ \phi\_{\mathbf{3}^k + \lambda + 1} \cdots \phi\_{(\mathbf{3}^k + \lambda) + \mathbf{3}^k} \cdots \phi\_{2 \times (\mathbf{3}^k + \lambda)} \ \phi\_{2 \mathbf{c} (\mathbf{3}^k + \lambda) + 1} \\\\ \cdots \phi\_{2 \times (\mathbf{3}^k + \lambda) + \mathbf{3}^k} \cdots \phi\_{\mathbf{3} \times (\mathbf{3}^k + \lambda)} \end{bmatrix}^\mathrm{T} \tag{8}$$

where

$$\Phi\_{u} = \begin{cases} \sum\_{j=1}^{j^{(1)}} f\_{1}^{(1)} + w \sum\_{j=1}^{\lambda} \theta\_{j}^{(1)} & 1 \le u \le 3^{\lambda} \\ \sum\_{j=1}^{\lambda^{\lambda}} f\_{1}^{(1)} + w \sum\_{j=1}^{\lambda^{\lambda}} \theta\_{j}^{(1)} & 3^{\lambda} + 1 \le u \le 3^{\lambda} + \lambda \\ \sum\_{j=1}^{\lambda^{\lambda}} f\_{1}^{(1)} + w \sum\_{j=1}^{\lambda^{\lambda}} \theta\_{j}^{(1)} & 3^{\lambda} + \lambda + 1 \le u \le 2 \times 3^{\lambda} + \lambda \\ \sum\_{j=3^{\lambda^{\lambda}}} \frac{w \theta\_{j}^{(2)}}{f\_{1}^{(2)} + w} + w \sum\_{j=1}^{\lambda} \theta\_{j}^{(2)} & 2 \times 3^{\lambda} + \lambda + 1 \le u \le 2 \times 3^{\lambda} + 2\lambda \\ \sum\_{j=3^{\lambda^{\lambda}}} \frac{w \theta\_{j}^{(3)}}{f\_{1}^{(2)} + w} + w \sum\_{j=1}^{\lambda} \theta\_{j}^{(3)} & 2 \times 3^{\lambda} + 2\lambda + 1 \le u \le 3 \times 3^{\lambda} + 2\lambda \\ \sum\_{j=3^{\lambda^{\lambda^{\lambda}}}} \frac{w \theta\_{j}^{(3)}}{f\_{1}^{(3)} + w} + w \sum\_{j=1}^{\lambda} \theta\_{j}^{(3)} & 3 \times 3^{\lambda} + 2\lambda + 1 \le u \le 3 \times 3^{\lambda} + 3\lambda \end{cases} \tag{9}$$

where λ denotes the highest tier correlation of the *k-*GCC nucleotides in each local window of **D**, whose the value is an integer. *w* is a float number that represents the weight factor, and the value of *w* is between 0 and 1. In the front window, the middle window and the rear window, the correlation factor of the *j*-th tier is represented as θ*j* ( ) 1 , θ*j* ( ) 2 , and θ*j* ( ) 3 , respectively. The GC skew value of the *k*-GCC nucleotides separated by *j* nucleotides is used to represent the correlation factor of the *j*-th tier in each local window. (**Figure 1**). θ*j* ( ) 1 , θ*j* ( ) 2 , and θ*j* ( ) 3 can be calculated by

$$\text{128}\begin{cases} \theta^{(i)}\_{j} = \frac{1}{\text{Int}^{c}[\frac{\mathbf{N}-k}{\beta}]+1} \sum\_{i=0}^{\text{Int}^{c}[\frac{\mathbf{N}-k}{\beta}]} \Theta(\mathbf{N}\_{\alpha j + 1}, \mathbf{N}\_{\alpha j + 2}, \dots, \mathbf{N}\_{\alpha j + k})\\ \qquad\qquad\qquad\qquad\qquad \qquad \qquad \qquad \text{128}\begin{cases} \Theta(\mathbf{N}\_{\alpha j + 1}, \mathbf{N}\_{\alpha j + 2}, \dots, \mathbf{N}\_{\alpha j + k})\\ \qquad\qquad\qquad \qquad \qquad \qquad \qquad \text{129}\begin{cases} \mathbf{N}\_{\alpha j} = \mathbf{0}, \mathbf{N}\_{\alpha j} \\ \end{cases} \end{cases} \\\ \Theta^{(i)}\_{j} = \frac{1}{\text{Int}^{c}[\frac{\mathbf{N}-k}{\beta}]+1} \sum\_{i=0}^{\text{Int}^{c}[\frac{\mathbf{N}-k}{\beta}]} \Theta(\mathbf{N}\_{\alpha j + k}, \mathbf{N}\_{\alpha j + 2}, \dots, \mathbf{N}\_{\alpha j + k})\\ \qquad\qquad \qquad \quad \quad \text{120}\begin{cases} \frac{\mathbf{N}-k}{\beta} \end{cases} \Theta(\mathbf{N}\_{\alpha j + k}, \mathbf{N}\_{\alpha j + 2}, \dots, \mathbf{N}\_{\alpha j + k}) \end{cases} \quad \text{121}\tag{10}$$

where Int C[ ] η − + *k j* 1 denotes the number of the *k*-GCC in the corresponding local window, and Θ(N*i* × *j* + 1N*i* × *j* + 2 ⋯ N*i* × *j* + *k*) is the GC Skew (Lobry, 1996; Grigoriev, 1998; Li et al., 2014) of the *i*-th *k*-GCC in the local window, which can be calculated by

$$\Theta\left(\mathbf{N}\_{\iota\times j+1}\mathbf{N}\_{\iota\times j+2}\cdots\mathbf{N}\_{\iota\times j+k}\right) = \frac{f\_G\left(\mathbf{N}\_{\iota\times j+1}\mathbf{N}\_{\iota\times j+2}\cdots\mathbf{N}\_{\iota\times j+k}\right) - f\_G\left(\mathbf{N}\_{\iota\times j+1}\mathbf{N}\_{\iota\times j+2}\cdots\mathbf{N}\_{\iota\times j+k}\right)}{f\_G\left(\mathbf{N}\_{\iota\times j+1}\mathbf{N}\_{\iota\times j+2}\cdots\mathbf{N}\_{\iota\times j+k}\right) + f\_G\left(\mathbf{N}\_{\iota\times j+1}\mathbf{N}\_{\iota\times j+2}\cdots\mathbf{N}\_{\iota\times j+k}\right)}\tag{11}$$

where *f*G(N*i* × *j* + 1 N*i* × *j* + 2 ⋯ N*i* × *j* + *k*) denotes the frequency of G in the subsequence N*i* × *j* + 1 N*i* × *j* + 2 ⋯ N*i* × *j* + *k f*C(N*i* × *j* + 1 N*i* × *j* + *k* ⋯ N*i* × *j* + *k*) denotes the frequency of C in the subsequence N*i* × *j* + 1 N*i* × *j* + 2 ⋯ N*i* × *j* + *k*, reflecting the CG asymmetry bias directly. Please note that for the terminal subsequence, if its length is less than *k*, then the GC skew will be calculated by all the available nucleotide residues.

#### Random Forest

Being widely used in bioinformatics (Zhao et al., 2014; Su et al., 2019), Random Forest (RF) (Ho, 1995; Barandiaran, 1998) is a machine learning classifier. Its training process can prevent overfitting (Hastie et al., 2008). The Random Forest model was implemented by calling the command line RandomForestClassifier ("max\_features='sqrt', min\_ samples\_leaf=1, min\_samples\_split=2, criterion = 'gini', = optimize-d value ") with the help of the Scikit-learn package (Pedregosa et al., 2011), where the values of represents the number of the trees in the forest, and it was set as 600 for both the two benchmark datasets (cf. Equation 1).

#### Ensemble Learning

Previous studies (Zou et al., 2015; Liu et al., 2016a; Chen et al., 2016b; Chen et al., 2017a; Chen et al., 2017b; Liu et al., 2018a) have demonstrated that fusing a series of individual predictors

by a voting strategy can improve the predictive performance. In this regard, in this study an ensemble predictors was constructed by fusing 10 top performing individual predictors constructed by different parameter combinations of PseKGCC (see **Supplementary Information S1**), which can be represented as (Liu et al., 2016a):

*k-*GCC (*k* = 4); (D) The coupling between the second most contiguous *k-*GCC (*k* = 4).

$$\mathbb{RF}^{\mathrm{E}} = \mathrm{RF}(1) \,\forall \, \mathrm{RF}(2) \,\forall \, \cdots \,\forall \, \mathrm{RF}(i) \,\forall \, \cdots \,\forall \, \mathrm{RF}(10) = \forall\_{i=1}^{10} \mathrm{RF}(i) \,\tag{12}$$

where E represents the ensemble classifier, ∀ represents the fusing operator, and RF(*i*) represents the basic Random Forest predictor.

The ensemble predictor is constructed based on the fusion score ß of the probabilities predicted by the 10 basic predictors, which can be calculated by

$$\mathfrak{A} = \sum\_{i=1}^{10} q\_i P\_i \tag{13}$$

where *qi* is the weight of the *i*-th basic RF predictor, which was optimized by the genetic algorithm (Mitchell, 1998), and their values were listed in **Supplementary Information S1**. If the value of ß is higher than 0.5, it is a replication origin, otherwise, it is a non-replication origin. The flowchart of the iRO-PseKGCC is illuminated in **Figure 2**.

FIGURE 2 | A flowchart illustration to show how the iRO-PseKGCC predictor works.

#### Cross Validation

Three widely used cross-validation strategies include: i) independent test, ii) K-fold cross validation, and iii) jackknife test. Among these methods, only the jackknife test can achieve the unique results for the same benchmark dataset. Therefore, in this study, the jackknife test was employed to give the final predictive results. However, considering its high computational cost, during the parameter optimization process, the 5-fold cross-validation was used to reduce the computational cost (see *Optimize Parameters* section).

#### Evaluation Method of Performance

To evaluate the quality of the classifier for prediction of the replication origins, the four metrics are used (Feng et al., 2013; Chen et al., 2016c; Chen et al., 2019): i) the sensitivity, Sn, ii) the specificity, Sp, iii) the overall accuracy of the predictive results, Acc, iv) the Mathew's correlation coefficient, MCC, and v) Arear under ROC Curve, AUC (Chen et al., 2016a), defined as:

where *N*+ denotes the number of all the positive samples (replication origins), *N*– denotes the number of all the negative samples (non-replication origins), *N*− + denotes the number of the positive samples (replication origins) incorrectly predicted as the negative samples (non-replication origins), *N*+ − denotes the number of the negative samples (non-replication origins) incorrectly predicted as the positive samples (replication origins). More information of these performance measures can refer to Liu et al. (2016b).

(14)

### RESULTS AND DISCUSSION

#### Optimize Parameters

There are five parameters in PseKGCC according to Equations 4–9. These parameters were optimized by the following equations:

The fivefold cross-validation was employed to search the optimal parameters by gridding method so as to reduce the time consumption, and the predictive results of the top 10 performing predictors, and their optimized parameters were listed in **Supplementary Information S1**.

#### Comparison With Other Methods

To the best knowledge of ours, iRO-3wPseKNC (Liu et al., 2018b) is the only existing predictor that is able to predict the entire replication origins. All the other predictors can only predict the fragments of replication origins. Therefore, the performance of the proposed iRO-PseKGCC was compared with iRO-3wPseKNC on the two benchmark datasets, and the results were listed in **Table 1**, from which we can see that iRO-PseKGCC obviously outperformed iRO-3wPseKNC in terms of the five performance measures (cf. Equation 14), indicating that the proposed PseKGCC feature is able to capture the GC asymmetry bias, and incorporate the GC skew into the predictor. Therefore, iRO-PseKGCC is an efficient approach for improving the predictive performance.

#### Feature Analysis

Random forest is a combination classifier model composed of decision tree classifiers. During the process of constructing each tree by the "Bootstrap" method (Efron, 1992), samples not extracted for training the corresponding tree can be used to make "Out Of Bag" (OOB) error estimate (Breiman, 1996) to evaluate the generalization performance of a predictor. Based on the OOB error, the Mean Decrease Accuracy (MDA) (Jiang et al., 2007) can

TABLE 1 | The results of the iRO-PseKGCC Predictor and comparison with iRO-PseKGCC on the two benchmark datasets (cf. Equation 1) obtained by using jackknife test.

*aThe parameters are listed in Supplementary Information S1.*

*bThe predictor reported in (Liu et al., 2018b) with parameter* ε *= 0.25,* δ *= 0.85, k = 5,* λ*= 6, w = 0.3, and = 700.*

*cThe predictor reported in (Liu et al., 2018b) with parameter* ε *= 0.15,* δ *= 0.55, k = 4,* λ *= 9, w = 0.3, and = 800.*

be used to estimate the importance of the features. The details of the process are (Jiang et al., 2007): 1) When training a Random Forest model, using the OOB samples to test the accuracy of each tree in the model; 2) Randomly disturb the value of the feature variable *v* in the OOB samples, and retest the accuracy of each tree; 3) Calculate the mean value of the decreasing accuracy between the two tests in all decision trees in the Random Forest model. The MDA value can reflect the importance of the corresponding feature.

As shown in previous studies (Liu and Zhu, 2019; Liu et al. 2019a), feature analysis is critical for exploring the characteristics of the predictors. To explore the reason why the proposed predictor iRO-PseKGCC works so well, we analyzed the features of the two top performing iRO-PseKGCC predictors (see **Supplementary Information S1**) on the two benchmark datasets (cf. Equation 1) by MDA approach, and the results are listed in the **Table 2**, from which we can see that: 1) for both the two RF-based predictors, their most important features are the "\*\*\*" and "\*\*\*\*\*," indicating the importance of the *k*-GCC; 2) The global sequence order effects measured by different λ values and GC skew values contribute to the performance improvement; 3) Features in certain local window show more discriminative powers than those in other windows, for examples, for *Pichia pastoris*, all the top 10 most important features are in the middle window, which is consistent with the previous observations that the nucleobase composition distribution is uneven along the replication origins (Lobry, 1996; Grigoriev, 1998; Frank and Lobry, 1999; Tillier and Collins, 2000; Liu et al., 2018b).

### Web Server and User Guide

Web-servers are important for the researchers to implement the corresponding computational predictors. In this regard, for the user's convenience, we established a web-server named

TABLE 2 | The top 10 most important features of the top two performing RF-based predictors on the two benchmark datasets (cf. Equation 1).

"iRO-PseKGCC." For users' convenience, a detailed user guide explaining how to use the web-server is given.

#### DATA AVAILABILITY

Publicly available datasets were analyzed in this study. These data can be found here: https://academic.oup.com/bioinformatics/ article-abstract/34/18/3086/4978052?redirectedFrom=fulltext

#### REFERENCES

#### AUTHOR CONTRIBUTIONS

BL provided the main idea of the manuscript and wrote the manuscript. SC did the experiments and revised the manuscript. KY revised the manuscript and did the typesetting. FW did the experiments.

#### FUNDING

This work was supported by the National Natural Science Foundation of China (No. 61672184, 61732012, 61822306), Fok Ying-Tung Education Foundation for Young Teachers in the Higher Education Institutions of China (161063), Shenzhen Overseas High Level Talents Innovation Foundation (Grant No. KQJSCX20170327161949608), Guangdong Natural Science Funds for Distinguished Young Scholars (2016A030306008), and Scientific Research Foundation in Shenzhen (JCYJ20180306172207178).

#### SUPPLEMENTARY MATERIAL

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

**Conflict of Interest Statement:** 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 © 2019 Liu, Chen, Yan and Weng. 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.*

# Variance-Preserving Estimation of Intensity Values Obtained From Omics Experiments

*Adèle H. Ribeiro1\*, Julia Maria Pavan Soler2 and Roberto Hirata Jr.1*

*1 Department of Computer Science, Institute of Mathematics and Statistics, University of São Paulo, São Paulo, Brazil, 2 Department of Statistics, Institute of Mathematics and Statistics, University of São Paulo, São Paulo, Brazil*

Faced with the lack of reliability and reproducibility in omics studies, more careful and robust methods are needed to overcome the existing challenges in the multi-omics analysis. In conventional omics data analysis, signal intensity values (denoted by *M* and values) are estimated neglecting pixel-level uncertainties, which may reflect noise and systematic artifacts. For example, intensity values from two-color microarray data are estimated by taking the mean or median of the pixel intensities within the spot and then subjected to a within-slide normalization by LOWESS. Thus, focusing on estimation and normalization of gene expression profiles, we propose a spot quantification method that takes into account pixel-level variability. Also, to preserve relevant variation that may be removed in LOWESS normalization with poorly chosen parameters, we propose a parameter selection method that is parsimonious and considers intrinsic characteristics of microarray data, such as heteroskedasticity. The usefulness of the proposed methods is illustrated by an application to real intestinal metaplasia data. Compared with the conventional approaches, the analysis is more robust and conservative, identifying fewer but more reliable differentially expressed genes. Also, the variability preservation allowed the identification of new differentially expressed genes. Using the proposed approach, we have identified differentially expressed genes involved in pathways in cancer and confirmed some molecular markers already reported in the literature.

Keywords: delta method, pixel-level uncertainty, spot quantification, optimal LOWESS normalization, two-color microarray, variability preservation, parameter selection

## INTRODUCTION

The growing number of omics datasets (e.g., genomics, transcriptomics, proteomics, metabolomics) and the recent advances in multi-omics integration approaches have contributed to the better understanding of biological mechanisms and also the emergence of the personalized medicine. However, the lack of reliability and reproducibility in omics studies stands as one of the biggest obstacles in bridging the gap between research and practice of personalized medicine (Alyass et al., 2015; Karczewski and Snyder, 2018). Considering that inflated variability and non-robust estimation may lead to inaccurate and misleading results, this paper proposes improvements to the conventional estimation and normalization of the intensity values obtained from omics experiments. Specifically, the proposal is to estimate the intensity values by a method that accounts for the variability due to pixel-level uncertainties and to normalize these values by using LOWESS with suitably selected

#### *Edited by:*

*Liang Cheng, Harbin Medical University, China*

#### *Reviewed by:*

*Tian Qing Zheng, Chinese Academy of Agricultural Sciences, China Tianyi Zhao, Harbin Institute of Technology, China*

> *\*Correspondence: Adèle H. Ribeiro<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Statistical Genetics and Methodology, a section of the journal Frontiers in Genetics*

> *Received: 05 April 2019 Accepted: 16 August 2019 Published: 20 September 2019*

#### *Citation:*

*Ribeiro AH, Soler JMP and Hirata R Jr (2019) Variance-Preserving Estimation of Intensity Values Obtained From Omics Experiments. Front. Genet. 10:855. doi: 10.3389/fgene.2019.00855*

parameter values, preserving variation that may be relevant to subsequent analyses.

Image processing and fluorescence analysis are the preferred approaches for data quantification in microarray technologies. Although microarrays have been predominantly used since the end of the nineties to measure gene expression levels, they remain widely used to detect other omics data types, including microRNA expression, DNA methylation, single-nucleotide polymorphisms (SNPs), and copy number variants (CNVs) (Goodwin et al., 2016). After hybridization and cleaning of the target molecules, the array is scanned by activation with lasers at different wavelengths (one for each of the fluorophores used), and each laser channel generates an image. The pixel intensities within each spot in these microarray images are summarized to represent the hybridization signal. Depending on the platform (e.g., gene expression array, DNA methylation array, SNP array, and comparative genomic hybridization [CGH] array), the interpretation of this signal is different (e.g., gene expression levels, methylation levels, allele frequencies, and copy number alterations).

The continuance of the microarray technology can be mainly explained by the availability of many datasets in public repositories, such as the Gene Expression Omnibus (GEO) (Edgar et al., 2002; Barrett et al., 2012) and ArrayExpress (Kolesnikov et al., 2015), by the existence of well-established strategies for data analysis and experimental design, and by the low cost compared with the next-generation sequencing technologies. However, given that microarray analysis is still facing reliability and reproducibility problems, more robust and rigorous methods are needed to account for the high variability and biases introduced in all steps of a microarray experiment.

Several preprocessing and normalization procedures have been proposed to remove biases due to the inhomogeneity of the background and the different fluorescence properties of the dyes. However, biases introduced in the image analysis step, which includes spot segmentation and signal extraction, have not received the same attention, and those may partially explain the existing reliability and reproducibility problems in omics studies. Particularly, several factors, including image resolution, scanner settings, effectiveness of the segmentation algorithm, and unexpected behaviors during hybridization, may lead to errors in spot localization and classification of the pixels (as foreground or background, depending on whether it is situated within or around the spot). Thus, spot intensities are usually noisy and that high pixel–level variability leads to uncertainty in microarray quantification and correlates with variability between replicate spots on duplicate slides (Brown et al., 2001).

Given that even state-of-art image processing tools are susceptible to errors that significantly influence the variability of the data derived from microarray images (Ahmed et al., 2004), new segmentation and intensity extraction algorithms are still being developed in order to improve precision in spot quantification (Li et al., 2017; Karthik and Manjunath, 2018; Shao et al., 2019). Usually, these tools combine sophisticated algorithms and pixellevel analyses in order to obtain an accurate estimate of the signal intensity in each spot. However, to allow subsequent analyses to take into account possible errors and uncertainties arising from the image processing, the method output usually includes not only statistical measures of location (e.g., mean and median) of the foreground and background intensities of each channel of each spot but also measures of dispersion, including standard deviation and covariance between both channels.

Despite the common use of pixel-level variability measures as data quality criteria for filtering purpose, the conventional microarray analysis is solely based on statistical measures of location of the spot intensities (Yang et al., 2002; Sun et al., 2011; Brady and Vermeesch, 2012). To improve robustness and reliability in microarray analysis, pixel-level uncertainties should be accounted for in the intensity log-ratio estimation and propagated to the next steps of the analysis.

Pixel-level uncertainties have been taken into account by many spot quantification algorithms in the literature, but requiring all pixel values to be available. Some of them are interested in improving the log-ratio estimator. Particularly, the method proposed by (Dodd et al., 2004) is a log-ratio estimator that corrects for signal saturation by regressing all pixel intensities at both test and control channels using a censored regression model. The META algorithm (Chan and Chang, 2009) estimates the intensity log-ratio by grouping the pixels according to their distance to the center of the spot and then weighting the log-ratio of each group in inverse proportion to its sample variance. A method that only uses pixel-level mean and variance summary statistics is the hierarchical maximumlikelihood estimator (Bakewell and Wit, 2005). However, it is not exactly based on the standard log-ratio representation of the spot intensity. It models the gene expression signal at control and treatment channels separately, incorporating the sample within-spot deviation and then performs the estimation using maximum likelihood. To the best of our knowledge, there is no intensity log-ratio estimator to be used after the image analysis phase (i.e., based solely on the pixel-level summary statistics) that takes into account pixel-level uncertainties.

The first contribution of this paper is a more robust estimator for the intensity log-ratio (*M*) and average log intensity (*A*) of a microarray spot that accounts for pixel-level variance and covariance between channels. For a spot *t*, these values are denoted by *Mt* and *At* , respectively (Dudoit et al., 2002). We derive these estimators by using the multivariate delta method (Casella and Berger, 1990). Specifically, we approximate the expected values of *Mt* and *At* by using their second-order Taylor's expansions, and the variance of *Mt* and *At* by using their firstorder Taylor's expansions. These expansions depend on the pixel-level variance and covariance between channels of the spot, whose sample estimates are readily accessible through standard output files of microarray image analysis tools.

After spot intensity estimation, it is necessary to perform a within-slide normalization to remove array-specific effects, intensity-dependent dye biases, and other systematic trends of the microarray data. The within-slide normalization based on the robust locally weighted regression (LOWESS) (Cleveland, 1979) is one of the most used techniques. The choice of the LOWESS parameters, particularly the smoothing parameter (also known as neighborhood size or bandwidth), dramatically affects the intensity and quality of the microarray data calibration. Although the smoothing parameter is still commonly set arbitrarily (around 0.2 and 0.4) (Dudoit et al., 2002; Smyth and Speed, 2003; Drăghici, 2012), some data-driven methods have been proposed to select its optimal value (Berger et al., 2004; Futschik and Crompton, 2004a; Lee et al., 2008). All these methods are similar in that they choose the smoothing parameter by minimizing a measure of error of the LOWESS fit. Berger et al. (2004) use the mean-squared difference between the LOWESS estimates and the corresponding normalization reference levels as cost function. These normalization levels are the true spot-specific calibration errors, which are usually unknown. Thus, Berger et al. suggest to estimate them from control transcripts and replicate slides. However, they are not always available for all genes in a typical microarray experiment, making it hard to reliably use the method. Futschik and Crompton's selection method, named OLIN (Futschik and Crompton, 2004a; Futschik and Crompton, 2004b), has the advantage of not relying on a reference level. Its optimization procedures use the generalized cross-validation (GCV) criterion, an estimator of the prediction mean square error (PMSE), as cost function. Lee et al. (2008) proposes to select the smoothing parameter by minimizing the bootstrap estimate of the mean integrated square error (MISE) and show that their results are comparable to OLIN.

Although all these methods have shown superiority over LOWESS normalization with a fixed arbitrarily chosen smoothing parameter, they lack in taking into account any heteroskedasticity in the data. In addition, they usually suffer from a poor bias–variance trade-off, tending to choose small smoothing values, which yield unnecessarily complicated (with high variance) LOWESS fits.

The second contribution of this paper is a data-driven method for selecting the smoothing parameter of the LOWESS normalization process. Inspired by the previous proposed methods, we choose the optimal smoothing value by minimizing a mean squared error criterion. However, our selection method also takes into account heteroskedasticity of the microarray data and offers a better bias–variance trade-off by selecting from among the low-MSE fits the one that is the most parsimonious. The parameter selection is obtained by solving a discrete optimization problem and is based on conventionally accepted ideas for analysis of M-plots—a graphical tool showing the curve of the MSE against the effective degrees of freedom of the estimate (Cleveland et al., 1988).

Given that the primary application of DNA microarrays has been to measure gene expression levels, we focus in this paper on variation-preserving estimation and normalization methods for gene expression levels from two-channel (or two-color) microarrays. However, it is straightforward to adapt the same ideas to improve analysis of other types of microarray data, even from single-channel technologies.

The proposed methods were evaluated by a differential gene expression analysis from real intestinal metaplasia and normal microarray samples. The proposed estimators for the *Mt* and *At* values were compared with the conventional estimators that neglect the pixel-level variability. In addition, we compared the proposed method for selecting the LOWESS smoothing parameter with OLIN, as it is conceptually similar to the other existing methods and can be applied even to microarray experiments with few or no replicates. Results show that a more robust and conservative analysis is performed when the LOWESS smoothing parameter is selected by our method, potentially reducing the number of false-positive differential expressions. Besides, both the pixel-level variabilities incorporated by the proposed estimators for the *Mt* and *At* values and the variability preserved by our more parsimonious normalization method contributed to the identification of new differentially expressed genes. Thus, the proposed methods may also reduce the falsenegative rate.

#### MATERIALS AND METHODS

Two procedures that critically affect the adequacy of microarray data analysis are the spot quantification, which extracts summarized quantitative measures of the pixel intensities within each spot of the microarray slide, and the within-slide normalization, which removes dye-specific biases and other systematic noises simultaneously from all logged spot intensities (*Mt* and *At* values).

In the section Intestinal Metaplasia Database, we describe a gene expression dataset used to illustrate the application of our proposed methods. In the section Improved Estimators for the *Mt* and *At* values, we show our improved estimation method for the *Mt* and *At* values that incorporates pixel-level variability. In the section Estimators for the Variances of the *Mt* and *At* Values, we discuss some criteria that can be used for proper setting of the parameters of the LOWESS within-slide normalization and we propose an algorithm for selecting the optimal value for the smoothing parameter (denoted by *f*).

### Intestinal Metaplasia Database

Due to a chronic inflammatory process, the normal squamous mucosa of the stomach may be replaced by columnar intestinaltype epithelium, characterizing a disease called intestinal metaplasia of the stomach. Since adenocarcinoma of the stomach and inflamed intestinal mucosa are strongly associated (Coussens and Werb, 2002), intestinal metaplasia may be a significant risk factor for gastric cancer.

We analyzed data from a two-color microarray experiment with tissues samples from 90 different subjects, being 35 from tissues representing type II intestinal metaplasia and 55 from tissues representing the normal condition, obtained from the Tumor Bank at A.C. Camargo Cancer Center/Antonio Prudente Foundation.

It was used the standard reference design (Churchill, 2002), in which each sample is hybridized against a pool of normal tissues using the same orientation of dye labeling. Gene expression levels were measured on Agilent Whole Human Genome Microarrays 4x44K G4112F (design ID 014850), each slide containing 41,093 unique probes. The scanned images of the microarray slides were processed by *Agilent Feature Extraction* software, version 9.5, where statistics (mean, standard deviation, and covariance) of the foreground and local background pixels were computed for each spot, in both test and reference channels. Each microarray spot contains about 60 foreground pixels.

This study was carried out in accordance with the recommendations of the international guidelines for investigations involving human beings with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the Ethics Institutional Committee of the A.C. Camargo Cancer Center (process number 1023/07).

#### Improved Estimators for the *Mt* and *At* Values

Usually, in microarray analysis, the test channel is denoted by (red), and the reference channel is denoted by *G* (green), following this usual notation, denoted by *Rtj* and by *Gtj*, the intensity value of the *j*th pixel within the th spot, respectively, in the test and reference channel. The relative expression of pixel *j* within spot is denoted by *Mtj* and defined as follows:

$$M\_{\circ j} \doteq \log\_2 \left( \frac{R\_{\circ j}}{G\_{\circ}} \right) = \log\_2(R\_{\circ j}) - \log\_2(G\_{\circ}).\tag{1}$$

 The average expression of pixel within spot is denoted by *Atj* and defined as follows:

$$A\_{tj} \doteq \frac{1}{2} \left( R\_{tj} G\_{tj} \right) = \frac{\log\_2(R\_{tj}) + \log\_2(G\_{tj})}{2}. \tag{2}$$

Usually, image analysis software does not provide all pixel intensity values within each spot. Nonetheless, it provides several descriptive statistics of the foreground and background pixel intensities, including sample estimates for the mean, median, variance, and covariance between the two channels.

To incorporate the pixel-level variability in the analysis, we derived an approximation of the expected values of *Mtj* and *Atj* by using the *multivariate delta method* (Casella and Berger, 1990). Assuming that the functions (1) and (2) are twice differentiable on an open interval which contains the point ( ( *R G* ), ) ( ) *tj tj* , we computed their second-order Taylor's expansions, around the point ( ( *R G* ), ) ( ) *tj tj* , and then derived their expected values. The derivation is presented in Appendix 4.

It is reasonable to assume that the variables *Rtj*, *Gtj*, *Mtj* and *Atj* have a distribution with well-defined mean and variance. Particularly, Hoyle et al. (Hoyle et al., 2002) empirically showed that the distribution of the pixels within a spot is heavy-tailed (a non-Gaussian distribution) and well-approximated by a log-normal distribution. Consequently, *Mtj* and *Atj* follow a distribution which is well-approximated by a Gaussian distribution and all the variables have at least the first and second moments finite.

Let *Rtc* and *Gtc* be non-zero estimates of, respectively, ( ) *Rtj* and ( ) *Gtj* , which represent average foreground signals after correction for removing the background influence. The subscript indicates dependence on the background correction. Also, let σˆ ( ) 2 *Rt* and σˆ ( ) 2 *Gt* be estimates of, respectively, Var (*Rtj*) and Var (*Gtj*), which are assumed to be independent of the background correction. Note that mean and variance estimates are calculated across observed foreground pixel intensities within the spot at the respective channel.

We can derive improved estimators for ( ) *Mtj* and ( ) *Atj* as follows:

$$\begin{split} \bar{\boldsymbol{M}}\_{t} \doteq \mathbb{E}(\boldsymbol{M}\_{t}) &= \log\_{2} (\overline{\mathcal{R}}\_{\scriptscriptstyle{\boldsymbol{\alpha}}}) - \log\_{2} (\overline{\mathcal{G}}\_{\scriptscriptstyle{\boldsymbol{\alpha}}}) \\ &+ \frac{1}{2\ln \left( 2 \right)} \bigg( -\frac{\hat{\sigma}^{2} \left( \boldsymbol{R}\_{t} \right)}{\overline{\mathcal{R}}\_{\scriptscriptstyle{\boldsymbol{\alpha}}}^{2}} + \frac{\hat{\sigma}^{2} \left( \boldsymbol{G}\_{t} \right)}{\overline{\mathcal{G}}\_{\scriptscriptstyle{\boldsymbol{\alpha}}}^{2}} \bigg), \end{split} \tag{3}$$

$$\begin{split} \tilde{A}\_{t} \doteq \mathbb{E}(A\_{t\circ}) &= \frac{1}{2} \Big( \log\_{2}(\overline{R}\_{\iota\iota}) + \log\_{2}(\overline{G}\_{\iota\iota}) \Big) \\ &- \frac{1}{4\ln\left(2\right)} \Big( \frac{\hat{\sigma}^{2}(\mathcal{R}\_{\iota})}{\overline{R}\_{\iota\iota}^{2}} + \frac{\hat{\sigma}^{2}(G\_{t})}{\overline{G}\_{\iota\iota}^{2}} \Big). \end{split} \tag{4}$$

Note that the conventional estimators for the *Mtj* and *Atj* values, given by

$$
\hat{M}\_t \doteq \log\_2(\overline{R}\_{\text{tc}}) - \log\_2(\overline{G}\_{\text{tc}}),
\tag{5}
$$

$$\hat{A}\_t \doteq \frac{\log\_2(\overline{R}\_{\mathfrak{t}\mathfrak{t}}) + \log\_2(\overline{G}\_{\mathfrak{t}\mathfrak{t}})}{2},\tag{6}$$

are approximations of, respectively, ( ) *Mtj* and ( ) *Atj* derived from only the zeroth-order Taylor's expansion of the functions that define *Mtj* and *Atj*. Thus, the conventional estimators ignore the known measures of pixel-variability, which represent uncertainties in the gene expression measurements.

**Figure 1** illustrates the differences between the estimators for the ( ) *Mtj* and ( ) *Atj* for a randomly chosen microarray slide of the database described in the section *Intestinal Metaplasia Database*. Since these estimators may suffer from numerical instability if the corrected foreground signals, *Rtc* and *Gtc*, are very close to zero, we removed the background influence by applying the *normexp* method (Ritchie et al., 2007) with offset equals to 50. The top 20 spots with the highest pixel-level variability are highlighted in red plus symbols. Several of these spots have low average intensity (small estimates for ( ) *Atj* ) and a small difference between the intensities of the two channels (estimates for ( ) *Mtj* close to zero), but they are not the majority. The differences between the proposed estimators, defined in Eq. (3) and (4), and the conventional estimators, defined in Eq. (5) and (6), are shown in **Figures 1C**, **D**. These differences are due to the distinct parts between their respective formulas. When computing the *M j* estimates, the ratio of the pixel-level variability to the squared expected value in the test channel appears in Eq. (3) with an opposite sign to the same term in the reference channel. Thus, positive and negative differences between the estimates for ( ) *Mtj* may occur if such terms do not cancel each other out. **Figure 1C** shows the *ilde M t* estimates were smaller than the *M*ˆ *t* estimates for the genes with highest pixel-level variance, indicating a larger variance in their test channels. **Figure 1D** shows some *A t* estimates were smaller than

influence was removed from the foreground signals by the *normexp* method with offset.

the ˆ *At* estimates. The reduction is explained by the fact that the additional terms in Eq. (4) are negative for any positive pixellevel variability in any channel.

#### Estimators for the Variances of the *Mt* and *At* Values

Since we have also available the sample covariance between *Rtj* and *Gtj*, denoted by σˆ( , *R G* ) *t t* , we applied the multivariate delta method for deriving estimators for the variances of the *Mtj* and *Atj*. We calculated the variance of the first order Taylor's expansion of the functions (1) and (2) that define, respectively, *Mtj* and *Atj*, as shown in Appendix 5. The variance estimators for *Mtj* and *Atj*, for pixels *j* within spot *t* are:

$$\hat{\boldsymbol{\sigma}}^{2}(\boldsymbol{M}\_{t}) \dot{=} \frac{1}{\ln^{2}(2)} \left( \frac{\hat{\boldsymbol{\sigma}}^{2}(\boldsymbol{R}\_{t})}{\overline{\boldsymbol{R}}\_{\boldsymbol{\kappa}}^{2}} + \frac{\hat{\boldsymbol{\sigma}}^{2}(\boldsymbol{G}\_{t})}{\overline{\boldsymbol{G}}\_{\boldsymbol{\kappa}}^{2}} - 2 \frac{\hat{\boldsymbol{\sigma}}(\boldsymbol{R}\_{t}, \boldsymbol{G}\_{t})}{\overline{\boldsymbol{R}}\_{\boldsymbol{\kappa}} \overline{\boldsymbol{G}}\_{\boldsymbol{\kappa}}} \right), \qquad (7)$$

$$\hat{\boldsymbol{\sigma}}^{2}(\boldsymbol{A}\_{\mathrm{t}}) \dot{\coloneqq} \frac{1}{4\ln^{2}(2)} \left( \frac{\hat{\boldsymbol{\sigma}}^{2}(\boldsymbol{R}\_{\mathrm{t}})}{\overline{\mathcal{R}}\_{\mathrm{u}}^{2}} + \frac{\hat{\boldsymbol{\sigma}}^{2}(\boldsymbol{G}\_{\mathrm{t}})}{\overline{\mathcal{G}}\_{\mathrm{u}}^{2}} - 2\frac{\hat{\boldsymbol{\sigma}}(\boldsymbol{R}\_{\mathrm{t}},\mathbf{G}\_{\mathrm{t}})}{\overline{\mathcal{R}}\_{\mathrm{u}}\overline{\mathcal{G}}\_{\mathrm{u}}} \right). \tag{8}$$

The variances of *Mtj* and represent pixel-level uncertainties of the th spot. They can be used, for instance, for assessing the quality of the th spot or for constructing confidence intervals for the parameters ( ) *Mtj* and ( ) *Atj* .

#### Optimal Selection of the LOWESS Parameters

To simplify the notation, we will denote the estimates for ( ) *Mtj* and ( ) *Atj* , independently of the estimation method used, by, respectively, *Mt* and *At* values.

It is necessary to remove from these *Mtj* intensity values the dependent dye-specific biases and other systematic errors by using some within-slide normalization method.

In the LOWESS within-slide normalization method, one estimates for each microarray slide a smoothing function µˆ that maps each *At* observed value to a smoothed *Mt* value, µˆ( ) *At* . Since µˆ( ) *At* is considered an estimate of a dye-dependent bias, it must be subtracted from the corresponding observed *Mt* value to obtain a residual value representing, presumably, the biologically relevant gene expression level.

An appropriate LOWESS estimation depends on the choice of its parameters. According to loader (Loader, 1999), the weight function and the number of iterations of the robustness algorithm are not critical parameters. Cleveland (Cleveland, 1979) comments that good choices for these parameters are, respectively, the tricube function and three iterations. However, the degree of the local polynomials and the smoothing parameter *f*, which, in the nearest neighbor method, is a number between and indicating the proportion of data used in each local fit, affects the bias and the variance of the fit.

Specifically, the higher the degree of the local polynomial (related to the complexity of the model), the lower the bias of the fit (probably, fitting the data very well). However, the additional parameters of this more complex model increase the variance of the fitted values, yielding a poor generalization ability (i.e., the model will have a large error). Thus, to avoid unstable LOWESS estimates, several references as (Loader, 1999; Yang et al., 2001; Dudoit et al., 2002; Smyth and Speed, 2003) recommend using local polynomials of degree one, mainly in the presence of sparsity, as is the case of microarray data.

The effects of the smoothing parameter *f* on the bias and variance of the fit are opposite to those of the degree of the local polynomials. Since the *f* parameter indicates the number of observations that will be used in the local polynomial estimation, when *f* value is large, a simple polynomial may not fit well to all observations in the neighborhood, distorting or ignoring essential features. In other words, the estimation of the smoothing function can be significantly biased. On the other hand, when a low *f* value is chosen, the number of observations may be insufficient to capture the general behavior of the data, resulting in a very noisy (large variance) fitness function.

In the next section, we propose a method for selecting a value for the *f* parameter, focusing on microarray data normalization. Our method takes into account the intrinsic characteristics of the bias and variance of the fit as well as of gene expression data.

#### Lowess Smoothing Parameter Selection

For microarray data normalization, the ideal LOWESS fitted curve captures only trends and effects from systematic errors, retaining all biological variation. However, it critically depends on the choice of the *f* parameter value.

**Figure 2** illustrates the MA plot of the microarray slide shown in **Figure 1B**, with different LOWESS fits yielded by *f* values varying from 0.05 to 0.9. The improved estimation method was used to obtain the *Mt* and *At* values, that is, the *M t* and *A t* estimates.

The quality of a LOWESS estimator can be assessed by the MSE, which measures how close the estimator µˆ is of the true mean function μ :

$$MSE(\hat{\mu}) = \mathbb{E}[(\mu - \hat{\mu})^2].$$

Since the real curve μ is unknown, we need a criterion to evaluate the MSE. Under the assumption of heteroskedasticity, Cleveland and Devlin (Cleveland and Devlin, 1988) propose the Mallows' Cp criterion for local fitting that can be used as as MSE estimator. In the presence of heteroskedasticity, as usual for microarray data, the heteroskedasticity-robust Cp (HRCp) criterion, proposed by Liu and Okui (Liu and Okui, 2013), may be a more appropriate MSE estimator. We detail this MSE estimator next.

Considering {(*A M*, )} *t t t T* =1 within-slide data points, the evaluation of the LOWESS smoothing function on any point is given by a linear combination of the observed points, whose weights {(*l A*( )} *t t T* =1 are assigned according to the distance of *A* to the *At* observed points:

$$\hat{\mu}(A) = \sum\_{t=1}^{T} l\_t(A)M\_t... $$

Consider the *T* × *T* matrix *L* which maps the observed to the fitted values:

$$\left( \begin{array}{c} \hat{\mu}(A\_1) \\ \vdots \\ \hat{\mu}(A\_T) \end{array} \right) = \mathbf{L}M = \left( \begin{array}{c} l\_1(A\_1)\dots l\_T(A\_1) \\ \vdots \\ l\_1(A\_T)\dots l\_T(A\_T) \end{array} \right) \left( \begin{array}{c} \mathcal{M}\_1 \\ \vdots \\ \mathcal{M}\_T \end{array} \right).$$

Two commons definitions of the effective degrees of freedom of µˆ are: (1) *v*1 tr (*L*) and (2) *v*2 tr (*L L*′ ), where tr stands for the trace operator.

Supposing that the variance of *Mt* , across *T* spots of a microarray slide, is constant and equals to σ2 , the Mallows' Cp for local fitting is defined as:

$$Cp(\hat{\mu}) = \frac{1}{\sigma^2} \sum\_{t=1}^{T} (M\_t - \hat{\mu}(A\_t))^2 - T + 2\nu\_1.$$

Cleveland et al. (1988) shows that σ2 can be estimated as follows:

$$
\hat{\sigma}^2 \doteq \frac{\Sigma\_{t=1}^T [M\_t - \hat{\mu}(A\_t)]^2}{n + \nu\_2 - 2\nu\_1}.
$$

When heteroskedasticity is present, Mallows' Cp criterion is not an appropriate MSE estimator. Considering the *T* × *T* diagonal matrix Σ, whose th diagonal element is given by a nonhomogeneous variance σ*t* 2 of *Mt* , a robust MSE estimation can be achieved by using the HRCp criterion, defined as:

$$HRCp(\hat{\mu}) = \sum\_{t=1}^{T} (M\_t - \hat{\mu}(A\_t))^2 + 2\operatorname{tr}(\Sigma L).$$

According to Loader (1999), σ*t* 2 can be estimated locally by calculating the error variance (the residual sum of squares divided by the corresponding degrees of freedom) of a nearly unbiased LOWESS fit, which can be yielded using a very small value for the smoothing parameter (e.g., *f* = 0.1. Since the local variance estimates can be very noisy, it may be appropriate to smooth them using a gamma kernel.

Several authors suggest to choose the *f* value which minimizes a measure of error of the LOWESS fit, such as the MSE criterion (Berger et al., 2004; Futschik and Crompton, 2004a; Lee et al., 2008). However, other authors (Mallows, 1973; Cleveland and Devlin, 1988; Loader, 1999) argue that a selection based only on minimizing the MSE criterion is a poor procedure since it ignores the intrinsic information of the bias and variance of the fit. Therefore, following their suggestion, we propose a method based on a graphical tool called M-plot. It is a graph of the MSE estimate as a function of the effective degrees of freedom of the fit.

M-plots illustrating the *f* parameter selection method for a typical microarray slide (ID<PHONE_NUMBER>95\_1.4) are shown in **Figure 3**. Dots show MSE estimates (by HRCp criterion) and respective degrees of freedom (by *v*2 definition) of LOWESS fits (on the *M*ˆ *t* and ˆ *At* estimates, in the first M-plot, and on the *M t* and *A t* , in the second M-plot) obtained with *f* parameter varying from to 0.2 We fixed the other LOWESS parameters (local polynomials of degree one, tricube weight function, and three iterations) so that the M-plot curve shows only the effect of the *f* parameter on the bias–variance compromise. Large *f* values tend to yield simple fits (with fewer degrees of freedom), which have a small variance, but a large bias. On the other hand, minimal *f* values tend to yield complex fits (with many degrees of freedom), which have a small bias, but a large variance.

For the microarray slide in **Figure 3**, a selection method based only on the minimization of the MSE curve would choose the smallest evaluated *f* value (0.2). However, any *f* value within the flattening region near to the minimum (the region with light-colored dots) is a good choice, in the sense that it yields a low-MSE fit (Cleveland and Devlin, 1988; Loader, 1999). Depending on the type of application, we can choose between one value which yields a low-bias fit (with more degrees of freedom) or a low-variance fit (with fewer degrees of freedom). Since we want to estimate a natural phenomenon behavior, we propose to select from the flattening region the *f* value which yields the simplest LOWESS fit (the one with fewest effective degrees of freedom). The biggest dot in each M-plot indicates the selected *f* value. The detection of the flattening region is made by searching points for which the derivative of the MSE curve is small. We check for each sequence of three points near the minimum whether the difference between the MSE values

251485069395\_1.4). The flattening region is represented by the light-colored dots and the selected *f* value by the biggest dot. The LOWESS fits were yielded using values ranging from 1 to 0.2 (from lowest to highest degree of freedom).

is small. If so, these points are considered as belonging to the flattening region.

The *f* parameter selection method can be summarized in the following discrete and constrained optimization problem. Consider a sequence of *l* different values for *f*, {*f*1, *f*2, ... , *fl* }, and denoted by µˆ *fk* , the LOWESS fit yielded by using the value *fk* for the *f* parameter. Also, let:

$$\begin{aligned} \mathcal{F} &= \{\hat{\mu}\_{f\_k}; f\_k \in \{f\_1, f\_2, \dots, f\_l\}, f\_{k+1} < f\_k, \text{ for } k = 1, \dots, l-1\}; \\ f\_{\text{min}} &= \operatorname\*{arg\,min}\_{f\_k} \text{HRCp } (\hat{\mu}\_{f\_k}), \text{ such that } \hat{\mu}\_{f\_k} \in \mathcal{F}; \\ f\_{\text{max}} &= \operatorname\*{arg\,max}\_{f\_k} \text{HRCp } (\hat{\mu}\_{f\_k}), \text{ such that } \hat{\mu}\_{f\_k} \in \mathcal{F}; \text{ and} \\ \Delta\_{\text{MSE}} &= 0.05(\text{HRCp } (\hat{\mu}\_{f\_{\text{max}}}) - \text{HRCp } (\hat{\mu}\_{f\_{\text{min}}})). \end{aligned}$$

Since *v*2 function provides the effective degrees of freedom of a given fit, the selected *f* value is the solution *f \** , if it exists, of the following problem:

$$\begin{aligned} f^{\star} &= \underset{f\_k}{\text{arg min }} \,\nu\_z(\hat{\mu}\_{f\_k}) \\ \text{subject to:} \\ &\hat{\mu}\_{f\_k} \in \mathcal{F}; \\ &\text{HRCp}\left(\hat{\mu}\_{f\_k}\right) \leq \text{ HRCp}\left(\hat{\mu}\_{f\_{\text{obs}}}\right) + \,\Delta\_{\text{MSE}}, \text{for } k = 1, 2; \\ &\text{HRCp}\left(\hat{\mu}\_{f\_{k-1}}\right) \leq \text{ HRCp}\left(\hat{\mu}\_{f\_{\text{obs}}}\right) + \,\Delta\_{\text{MSE}}, \text{for } k = 3, \dots, l; \\ &\left| \text{HRCp}\left(\hat{\mu}\_{f\_k}\right) \right| - \,\text{HRCp}\left(\hat{\mu}\_{f\_{k-1}}\right) \right| < \Delta\_{\text{MSE}}, \text{ for } k = 2, \dots, l; \text{and} \\ &\left| \text{HRCp}\left(\hat{\mu}\_{f\_k}\right) \right| - \,\text{HRCp}\left(\hat{\mu}\_{f\_{k-2}}\right) \right| < \Delta\_{\text{MSE}}, \text{ for } k = 3, \dots, l. \end{aligned}$$

If the minimum of the M-plot curve is far away of the point corresponding to the second lowest MSE estimate, the previous problem has no solution. In that case, the *f* value that yields the fit with lowest MSE estimate is selected. Specifically, the *f* parameter value is selected by solving the following problem:

$$f^\star \doteq \operatorname\*{arg\,min}\_{f\_k} \text{HRCp } (\hat{\mu}\_{f\_k}) \text{, such that } \hat{\mu}\_{f\_k} \in \mathcal{F}.$$

#### APPLICATION ON INTESTINAL METAPLASIA DATA

To investigate the effects of the proposed methods, we preprocessed the data described in the section *Intestinal Metaplasia Database* by using all discussed methods and compared the identified differentially expressed genes.

First, we applied the *normexp* method with offset value of for removing the background influence. Then, we compute the *Mt* and *At* values both by the conventional estimation methods, defined in Eq. (5) and (6), and by the proposed estimation methods, defined in Eq. (3) and (4). The LOWESS withinslide normalization was carried out as discussed in the section *Optimal Selection of the LOWESS Parameters*. For comparison purpose, the *f* smoothing parameter was selected both by the OLIN method (considered by us as a conventional approach) and by the proposed method, discussed in the section *LOWESS Smoothing Parameter Selection*. Since data from all microarray slides present overdispersion, we used the HRCp criterion as cost function of our selection method.

Therefore, the following four preprocessing procedures were applied separately to the original data:

**Figure 4** shows the distribution of the optimal values for the LOWESS *f* parameter, according to the proposed selection method with HRCp criterion, for the entire database, separated by normal and intestinal metaplasia conditions (both, hybridized against a pool of normal tissues). In the first plot, the LOWESS curve was fitted on the *M*ˆ *t* and ˆ *At* estimates and, in the second plot, on the *M t* and *A t* estimates. The average of the selected *f* values was close to 0.5.

As expected from a method that neither takes into account heteroskedasticity of the data nor attempts to make a good balance between bias and variance, the OLIN method selected the smallest evaluated value (0.2) for most of the slides. Same results were obtained when the *Mt* and *At* values were estimated by the conventional and by the proposed estimator. Such behavior has been reported in the literature, implying that the optimal *f* values according OLIN are usually close to the default one (Chiogna et al., 2009).

After preprocessing the data, a two-sample t-test assuming unequal variance was performed for each spotted gene to determine whether its expression is statistically different between gastric tissues in normal and intestinal metaplasia groups. However, since we are interested in directly assessing the impact of each proposed method on the t-statistics and p-values rather than making inference about differential expression, the comparative study was performed before applying a multiple testing correction.

#### Comparison of the Results

Results of a pairwise comparison among the p-values and t-statistics obtained by the four preprocessing methods are shown in **Figure 5**. In the left-column plots, we compare the p-values and, in the right-column plots, we show the changes in the difference between the group means (the absolute value of the t-statistic numerator) and in the within-group variability (the denominator of the t-statistic). Only genes with p-value less than 5% were considered.

The left-column plots show that most of the points are distributed around the 45-degree line. Thus, the p-values and, consequently, the total number of differentially expressed genes, even at a lower significance level, were similar among the four methods.

The first- and second-row plots show how p-values and t-statistics were affected by estimating the *Mt* and *At* values with the proposed method, which takes into account the pixellevel uncertainties. The genes represented by blue plus signs were identified as differentially expressed only when using the proposed estimator for the *Mt* and *At* values.

The genes represented by green crosses were identified as differentially expressed only when using the conventional estimator for the *Mt* and *At* values.

When the LOWESS *f* parameter is selected by OLIN (first-row plots), it is clear that the within-group variability decreases when using the proposed estimators for the *Mt* and *At* values. When the LOWESS parameter is selected by our method (second-row plots), there is still a reduction in the within-group variability. However, this impact is less clear because of the variability introduced when the LOWESS *f* parameter is selected by our method.

The third- and fourth-row plots compare p-values and t-statistics obtained by OLIN and the proposed approach for selecting the LOWESS *f* parameter. The genes represented by blue plus signs were identified as differentially expressed only when *f* was selected by the proposed method. The genes represented by green crosses were identified as differentially expressed only when selecting *f* by OLIN. It is clear that, for most genes, both within-group variabilities increased, implying that the normalization procedure was more conservative, and thus, more potentially relevant information is retained. In addition, for many genes, the increase in the within-group variability was counterbalanced by an increase in the distance between the groups. Such effect is even most pronounced when the proposed estimator for the *Mt* and *At* values are used. Thus, their respective p-values reduced enough to consider them as differentially expressed genes.

The diagrams in **Figure 6** show a comparison of the methods with respect to the total number of genes with p-value less than 5%. On the left, the p-values were not corrected for multiple tests, while on the right, the p-values were adjusted by the false discovery rate (FDR) correction (Benjamini and Hochberg, 1995).

Note that the four methodologies are quite different in terms of which genes were identified as differentially expressed. As a consequence of the more conservative (milder) noise reduction performed in the LOWESS withinslide normalization procedure with *f* parameter selected by our method, fewer genes are identified as differentially expressed. However, regardless of the normalization method, more genes could be identified as differentially expressed when the *Mt* and *At* values were estimated by the proposed estimation method that incorporates pixel-level variability. Given that both proposed methods make the analysis more robust by incorporating and preserving information neglected by the conventional methods, we can argue that they are contributing to the reduction of both false-positive and falsenegative rates.

#### Validation Analysis

To check the consistency of our analysis, we compared our results with those reported in the literature. Out of the genes which are associated with intestinal metaplasia according to the Gene Expression Omnibus platform (Edgar et al., 2002) of the NCBI (National Center for Biotechnology Information) website, 80 spotted genes (corresponding to 63 unique genes) have p-value (before FDR correction) less than 5%, and 35 spotted genes (corresponding to 29 unique genes) have p-value (after FDR correction) less than 5%. These findings are summarized respectively in **Tables 1**, **2**. In addition, **Figure 7** compares the total number of validated genes identified by each method with p-value less than 5% (before FDR correction).

Greater differences in inference were observed among the genes whose p-value is close to the significance level. These

column plots compare the difference between the absolute values of the numerators with the difference between the denominators of the t-test statistic.

genes have a more subtle differential expression, which can be easily damaged by measurement errors and poor estimation and normalization methods. Thus, the more accurate and careful analysis provided by the proposed methods is especially important for making decisions on the differential expression of these more sensitive genes.

Two replicates of the HSPB1 gene could not be identified as differentially expressed when using both the conventional estimators for the *Mt* and *At* values and our selection method for the LOWESS *f* parameter. Thus, the estimation of the *Mt* and *At*

values by the proposed estimators was crucial in determining the differential expression of the HSPB1 gene.

The genes PTEN, CTNNB1, MLH1, CXCR4, and CXCR1 could only be identified as differentially expressed when the LOWESS parameter was selected by our proposed method. Particularly, the gene CXCR4 only was determined as differentially expressed when the improved estimators for the *Mt* and *At* values were also used. In contrast, the gene KRT14 was no longer identified as differentially expressed when the LOWESS *f* parameter was selected by our proposed method.

In the following, we briefly describe the association of those genes with intestinal metaplasia of the stomach according to the literature data:

in esophageal adenocarcinoma when compared to normal esophagus (Lv et al., 2019).

### Genes Involved in Cancer

By performing a gene enrichment analysis, we identified, at a significance level of 5% (after FDR correction), 31 differentially expressed genes that are involved in cancer. Their respective p-values and fold changes are shown in **Table 3**. We remark that their association with intestinal metaplasia has not been clearly demonstrated yet. Thus, further investigation has to be done to confirm such conclusions.

Particularly, two replicates of the CCND1 gene and the LAMB2 gene were identified as differentially expressed only by the conventional approaches, suggesting that they may be false positives. Next, we briefly describe their association with cancer:

### DISCUSSIONS

Faced with the growing trend of multi-omics data integration in the midst of a replication crisis, improved microarray data analyses are crucial to identifying more reliable results (Ritchie et al., 2015a).

TABLE 1 | Genes reported in the literature as associated with intestinal metaplasia of the stomach that were identified as differentially expressed in our analysis at a significance level of 5% (after FDR correction).

*Each column shows the p-value (p), the FDR-corrected p-value (adj. p), and the fold change (FC) obtained in a variant of the database. P-values greater than 5% are shown in bold type.*

Ribeiro et al.

TABLE 2 | Other genes reported in the literature as associated with intestinal metaplasia of the stomach that were identified as differentially expressed in our analysis at a significance level of 5% (without FDR correction).

*Each column shows the p-value (p), the FDR-corrected p-value (adj. p), and the fold change (FC) obtained in a variant of the database. P-values greater than 5% are shown in bold type.*

Ribeiro et al.

Given that several pixel-level summary statistics are readily available in microarray databases, but are usually discarded in conventional approaches, we propose an improved estimation method for the *Mt* and *At* values, which takes into account the pixel-level variability. Specifically, we applied the multivariate delta method to derive estimators for the expected values of *Mt* and *At* , considering their Taylor's expansion up to the second-order terms. The conventional estimators, nonetheless, approximate the expected values considering only the zeroth-order term. Since the functions that define *Mt* and *At* are analytic (they are combinations of logarithmic function through addition or subtraction), the higher the number of terms of the Taylor expansion, the better the approximation of the function. Thus, we expect that the proposed estimators provide a better quantification of the hybridization signal. Also, by using these improved estimators, pixel-level dispersion can play an essential role in the analysis, increasing reliability.

To minimize the propagation of errors, the *Mt* and *At* values have to be properly normalized. Thus, we also propose a method for selecting the LOWESS smoothing parameter *f* that provides an optimal bias–variance compromise, considering some specific characteristics of microarray experiments, such as heteroskedasticity. This optimal normalization method leads to a more parsimonious correction of the systematic biases and, consequently, to greater preservation of the biological variation of interest.

By using the proposed methods, more variability information is considered and retained, improving inferences and preventing false conclusions. Thus, we expect to perform a more conservative analysis, where possibly fewer but more reliable differentially expressed genes are identified. In other words, we expect a reduction in both the false-positive and false-negative error rates.

Besides the theoretical support, relevant empirical observations could be drawn by a comparative study between the methods using real intestinal metaplasia microarray data. The results shows that inferences on differential gene expression were moderately affected by the incorporation of the pixel-level variability in the estimation of the *Mt* and *At* values and significantly affected by the LOWESS within-slide normalization using a smoothing parameter selected by the method. Both proposed methods tend to increase the within-group variability (the denominator of the t-statistic). However, for many genes, such increase occurred along with an increase in the difference between the group means (the absolute value of the t-statistic numerator), significantly reducing their respective p-values. Thus, many genes were identified as differentially expressed only when the proposed methods were used and some of them have been validated by other studies.

It is important to remark that most of the genes reported in the literature as differentially expressed in intestinal metaplasia were validated with a very strong association with the disease. Thus, these genes are probably more robust to difference approaches for estimating and normalizing the gene expression levels. On the other hand, genes sensitive to methods that address essential uncertainties in measurements are precisely those plagued with major reproducibility issues. Measurement error is one of the most damaging sources of error and has been neglected in many published analyses, thereby increasing uncertainty in parameter estimates and even inflating the estimates of effect sizes (Loken and Gelman, 2017). Thus, particularly for those sensitive genes, a more robust analysis is needed so that false conclusions are not made.

In this paper, we focused on gene expression from twocolor microarray data, but it is possible to use the same ideas to improve estimation and normalization of any fluorescent signal quantified by microarray image analysis. Also, the proposed methods could be adapted for oligonucleotide (one-color) microarray data. Particularly, the cyclic LOWESS normalization method (Bolstad et al., 2003) could be extended by just considering that the *Mt* and *At* values are defined by comparing pairs of arrays instead of pairs of channels and that the LOWESS normalization is applied to all distinct combination of two arrays. Although not so straightforward, it is also possible to adapt our methods to handle next-generation sequencing (NGS) data. Recently, Law et al. (Law et al., 2014) showed that RNA-Seq counts after log transformation and normalization by sequencing depth (log-counts per million, or log-cpm) can be properly analyzed by methods based on the normal distribution if a precision weight for each observation is taken into account. It was used to adapt all methods in the limma package (initially developed for microarrays) to also handle RNA-Seq and other sequence count data (Ritchie et al., 2015b). Therefore, considering the current need for accounting and propagating measurement uncertainties through analyses of NGS data (O'Rawe et al., 2015), a possible future work is to adapt our ideas to improve transcriptome profiling from RNA-Seq data.

 *Each column shows the p-value (p), the FDR-corrected p-value (adj. p), and the fold change (FC) obtained in a variant of the database. P-values greater than 5% are shown in bold type.*

LAMB2 5.18 × 10−3 6.76 × 10−2 −0.52 2.58 × 10−3 4.25 × 10−2 −0.49 4.42 × 10−3 6.10 × 10−2 −0.49 2.61 × 10−3 4.28 × 10−2 −0.49

to each base call, usually provided by base-calling algorithms, into the log-cpm estimator, leading to a more accurate gene expression quantification from RNA-Seq data.

### DATA AVAILABILITY

The omicsMA R package contains the source code of the proposed methods and part of the metaplasia dataset analyzed in this study. It was implemented using R, version 3.5.1, and depends on the locfit (Loader, 2013), maigesPack (Esteves et al., 2016), and limma (Ritchie et al., 2015b) R packages. The omicsMA R package is available at https://github.com/adele/ omicsMA, and the latest release is available at https://github. com/adele/omicsMA/releases/latest.

### ETHICS STATEMENT

This study was carried out in accordance with the recommendations of the international guidelines for investigations involving human beings with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the Ethics Institutional Committee of the A.C. Camargo Cancer Center (process number 1023/07).

### REFERENCES

### AUTHOR CONTRIBUTIONS

AR and RH conceived of the presented ideas. AR derived the models, implemented the methods, and analyzed the data. AR wrote the manuscript with support from RH and JS. All authors discussed the results and contributed to the final manuscript. RH and JS supervised the project.

### FUNDING

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001; National Council of Technological and Scientific Development (CNPq); and NAP eScience - PRP - USP. It was also supported by the Foundation for Research Support of the State of São Paulo (FAPESP) [grants 06/03227-2, 2011/50761-2 and 2015/01587-0].

### ACKNOWLEDGMENTS

We are greatly thankful to Luiz F. L. Reis, Director of the Research and Education, Hospital Sírio-Libanês, São Paulo, for making the intestinal metaplasia database available to us, and to Professors Eduardo Jordão Neves and Luís Gustavo Estevesfor helping us with data preprocessing and analysis.

**Conflict of Interest Statement:** 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 © 2019 Ribeiro, Soler and Hirata. 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.*

### APPENDIX

### Estimation of *E*(*Mtj*) and *E*(*Atj*) by the Delta Method

Let *f* (*Rtj*, *Gtj*) be a twice differentiable function of two random variables, *Rtj* and *Gtj*.The second-order Taylor's expansion of at ( ( *R G* ), ( )) *tj tj* is:

$$\begin{split} &f(R\_{\eta},G\_{\eta}) \sim f(\mathbb{E}(R\_{\eta}),\mathbb{E}(G\_{\eta})) + \frac{\partial f}{\partial R\_{\eta}}(\mathbb{E}(R\_{\eta}),\mathbb{E}(G\_{\eta}))(R\_{\eta}-\mathbb{E}(R\_{\eta})) + \\ &\frac{\partial f}{\partial G\_{\eta}}(\mathbb{E}(R\_{\eta}),\mathbb{E}(G\_{\eta}))(G\_{\eta}-\mathbb{E}(G\_{\eta})) + \\ &\frac{1}{2}\Big{(}\frac{\partial^{2}f}{\partial R\_{\eta}^{2}}(\mathbb{E}(R\_{\eta}),\mathbb{E}(G\_{\eta}))(R\_{\eta}-\mathbb{E}(R\_{\eta}))^{2} + 2\frac{\partial^{2}f}{\partial R\_{\eta}\partial G\_{\eta}}(\mathbb{E}(R\_{\eta}),\mathbb{E}(G\_{\eta})) \\ &\Big{[}(R\_{\eta}-\mathbb{E}(R\_{\eta}))(G\_{\eta}-\mathbb{E}(G\_{\eta}))\Big{]} + \frac{\partial^{2}f}{\partial G\_{\eta}^{2}}(\mathbb{E}(R\_{\eta}),\mathbb{E}(G\_{\eta}))(G\_{\eta}-\mathbb{E}(G\_{\eta}))^{2}\Big{]}. \end{split}$$

An approximation of ( ( *f R*( ) *G* ), *tj tj* can be determined by the expected value of the second-order Taylor's expansion of *f*:

$$\mathbb{E}(f(R\_{\eta}, G\_{\eta})) \sim \mathbb{E}[f(\mathbb{E}(R\_{\eta}), \mathbb{E}(G\_{\eta}))] + \frac{\partial f}{\partial R\_{\eta}}(\mathbb{E}(R\_{\eta}), \mathbb{E}(G\_{\eta})) \mathbb{E}(R\_{\eta} - \mathbb{E}(R\_{\eta})) + \varepsilon$$

$$\frac{\partial f}{\partial G\_{\eta}}(\mathbb{E}(R\_{\eta}), \mathbb{E}(G\_{\eta})) \mathbb{E}(\mathbf{G}\_{\eta} - \mathbb{E}(G\_{\eta})) +$$

$$\frac{1}{2} \left( \frac{\partial^{2} f}{\partial R\_{\eta}^{2}}(\mathbb{E}(R\_{\eta}), \mathbb{E}(G\_{\eta})) \mathbb{E}[(R\_{\eta} - \mathbb{E}(R\_{\eta}))^{2}] +$$

$$2 \frac{\partial^{2} f}{\partial R\_{\eta} \partial G\_{\eta}}(\mathbb{E}(R\_{\eta}), \mathbb{E}(G\_{\eta})) \mathbb{E}((R\_{\eta} - \mathbb{E}(R\_{\eta}))(G\_{\eta} - \mathbb{E}(G\_{\eta}))) +$$

$$\frac{\partial^{2} f}{\partial G\_{\eta}^{2}}(\mathbb{E}(R\_{\eta}), \mathbb{E}(G\_{\eta})) \mathbb{E}((\mathbf{G}\_{\eta} - \mathbb{E}(G\_{\eta}))^{2}) \right).$$

Considering that

$$\begin{aligned} \text{Var}\left(R\_{\circ}\right) &= \mathbb{E}[\left(R\_{\circ} - \mathbb{E}(R\_{\circ})\right)^{2}], \\ \text{Var}\left(G\_{\circ}\right) &= \mathbb{E}[\left(G\_{\circ} - \mathbb{E}(G\_{\circ})\right)^{2}], \text{and} \\ \text{Cov}\left(R\_{\circ}, G\_{\circ}\right) &= \mathbb{E}[\left(R\_{\circ} - \mathbb{E}(R\_{\circ})\right)(G\_{\circ} - \mathbb{E}(G\_{\circ}))], \end{aligned}$$

the following simplified expression for the expected value of *f* (*Rtj*, *Gtj*) is obtained:

$$\mathbb{E}(f(R\_{\eta}, G\_{\eta})) = f(\mathbb{E}(R\_{\eta}), \mathbb{E}(G\_{\eta})) + \frac{1}{2} \left( \frac{\partial^2 f}{\partial R\_{\eta}^2} (\mathbb{E}(R\_{\eta}), \mathbb{E}(G\_{\eta})) \text{Var} \,(R\_{\eta}) + 1 \right)$$

$$2 \frac{\partial^2 f}{\partial R\_{\eta} \partial G\_{\eta}} (\mathbb{E}(R\_{\eta}), \mathbb{E}(G\_{\eta})) \text{Cov } (R\_{\eta}, G\_{\eta}) +$$

$$\frac{\partial^2 f}{\partial G\_{\eta}^2} (\mathbb{E}(R\_{\eta}), \mathbb{E}(G\_{\eta})) \text{Var} \,(G\_{\eta}) \text{.}$$

Since

$$M\_{\mathfrak{q}} = f(R\_{\mathfrak{q}}, G\_{\mathfrak{q}}) \doteq \log\_2(R\_{\mathfrak{q}}) - \log\_2(G\_{\mathfrak{q}}),$$

the first and second derivatives of the function that defines *Mtj* are:

$$\begin{aligned} \frac{\partial f}{\partial R\_{lj}} &= \frac{1}{R\_{lj} \ln \text{ (2)}}; \qquad \frac{\partial f}{\partial G\_{lj}} = -\frac{1}{G\_{lj} \ln \text{ (2)}};\\ \frac{\partial^2 f}{\partial R\_{lj}^2} &= -\frac{1}{R\_{lj}^2 \ln \text{ (2)}}; \qquad \frac{\partial^2 f}{\partial G\_{lj}^2} = \frac{1}{G\_{lj}^2 \ln \text{ (2)}}; \qquad \frac{\partial^2 f}{\partial R\_{lj} \partial G\_{lj}} = 0. \end{aligned}$$

Assuming that ( ) *Rtj* and ( ) *Gtj* are non-zero, an approximation of ( ) *M R* (log ( ) log (*G* )) *tj* = − 2 2 *tj tj* can be obtained by using its second-order Taylor's expansion:

$$\begin{split} \mathbb{E}(M\_{\eta}) &= \mathbb{E}(\log\_{2}(R\_{\eta}) - \log\_{2}(\mathbb{G}\_{\eta})) \\ &\approx \log\_{2}(\mathbb{E}(R\_{\eta})) - \log\_{2}(\mathbb{E}(\mathcal{G}\_{\eta})) + \frac{1}{2} \left( -\frac{Var\,\,(\mathcal{R}\_{\eta})}{\ln\left(2\right)\mathbb{E}^{2}\left(R\_{\eta}\right)} + \frac{Var\,\,(\mathcal{G}\_{\eta})}{\ln\left(2\right)\mathbb{E}^{2}\left(\mathcal{G}\_{\eta}\right)} \right) \\ &= \log\_{2}(\mathbb{E}(R\_{\eta})) - \log\_{2}(\mathbb{E}(\mathcal{G}\_{\eta})) + \frac{1}{2\ln\left(2\right)} \left( -\frac{Var\,\,(R\_{\eta})}{\mathbb{E}^{2}\left(R\_{\eta}\right)} + \frac{Var\,\,(\mathcal{G}\_{\eta})}{\mathbb{E}^{2}\left(G\_{\eta}\right)} \right). \end{split}$$

Let the non-zero background-corrected signals be estimators for the expected values of the foreground signals, i.e.,

$$
\hat{\mathbb{E}}(R\_{\sharp j}) = \overline{R}\_{\omega} \text{, with } \overline{R}\_{\iota \varepsilon} \neq 0,
$$

$$
\hat{\mathbb{E}}(G\_{\iota j}) = \overline{G}\_{\iota \varepsilon} \text{, with } \overline{G}\_{\iota \varepsilon} \neq 0.
$$

Denote the sample variance estimators, obtained across the pixel intensities within each spot, as σˆ ( ) 2 *Rt* (for the test channel) and σˆ ( ) 2 *Gt* (for the control channel). Also, assume that these estimators do not depend on thebackground correction. We can derive an estimator for ( ) *Mtj* as follows:

$$
\bar{\mathcal{M}}\_t \doteq \log\_2(\overline{R}\_{\underline{n}}) - \log\_2(\overline{G}\_{\underline{n}}) + \frac{1}{2\ln \text{ (2)}} \left( -\frac{\hat{\sigma}^2(R\_t)}{\overline{R}\_{\underline{n}}^2} + \frac{\hat{\sigma}^2 \left( G\_t \right)}{\overline{G}\_{\underline{n}}^2} \right).
$$

Since

$$A\_{\mathfrak{t}\mathfrak{f}} = f(R\_{\mathfrak{t}\mathfrak{f}}, G\_{\mathfrak{t}\mathfrak{f}}) \doteq \frac{\log\_2(R\_{\mathfrak{t}\mathfrak{f}}) + \log\_2(G\_{\mathfrak{t}\mathfrak{f}})}{2},$$

we can estimate ( ) *Atj* in a similar way to ( ) *Mtj* . The first and second derivatives of *Atj* are:

$$\begin{aligned} \frac{\partial f}{\partial R\_{\circ j}} &= \frac{1}{2\ln\left(2\right)R\_{\circ j}}; \qquad \qquad \frac{\partial f}{\partial G\_{\circ j}} = \frac{1}{2\ln\left(2\right)G\_{\circ j}};\\ \frac{\partial^2 f}{\partial R\_{\circ j}^2} &= -\frac{1}{2\ln\left(2\right)R\_{\circ j}^2}; \qquad \qquad \frac{\partial^2 f}{\partial G\_{\circ j}^2} = -\frac{1}{2\ln\left(2\right)G\_{\circ j}^2}; \qquad \qquad \frac{\partial^2 f}{\partial R\_{\circ j}\partial G\_{\circ j}} = 0, \end{aligned}$$

An approximation of ( ) *Atj* is obtained by using its secondorder Taylor's expansion:

$$\begin{split} \mathbb{E}(A\_{\eta}) &\approx \frac{1}{2} (\log\_{2} \left( \mathbb{E}(R\_{\eta}) \right) + \log\_{2} \left( \mathbb{E}(G\_{\eta}) \right)) + \frac{1}{2} \bigg( -\frac{Var \left( R\_{\eta} \right)}{2 \ln \left( 2 \right) \mathbb{E}^{2} \left( R\_{\eta} \right)} - \frac{Var \left( G\_{\eta} \right)}{2 \ln \left( 2 \right) \mathbb{E}^{2} \left( G\_{\eta} \right)} \bigg) \\ &= \frac{1}{2} \Big( \log\_{2} \left( \mathbb{E}(R\_{\eta}) \right) + \log\_{2} \left( \mathbb{E}(G\_{\eta}) \right) \bigg) - \frac{1}{4 \ln \left( 2 \right)} \left( \frac{Var \left( R\_{\eta} \right)}{\mathbb{E}^{2} \left( R\_{\eta} \right)} + \frac{Var \left( G\_{\eta} \right)}{\mathbb{E}^{2} \left( G\_{\eta} \right)} \right). \end{split}$$

Considering the sample estimators of the expected values and variances of *Rtj* and *Gtj*, we can derive the following estimator for ( ) *Atj* :

$$\tilde{A}\_t \doteq \frac{1}{2} \Big( \log\_2(\overline{R}\_{\scriptscriptstyle\rm tr}) + \log\_2(\overline{G}\_{\scriptscriptstyle\rm tr}) \Big) - \frac{1}{4\ln \text{ (2)}} \Big( \frac{\hat{\mathfrak{G}}^2(R\_t)}{\overline{R}\_{\scriptscriptstyle\rm tr}^2} + \frac{\hat{\mathfrak{G}}^2(G\_t)}{\overline{G}\_{\scriptscriptstyle\rm tr}^2} \Big).$$

#### A.2.Estimation of *Var* (*Mtj*) and *Var* (*Atj*) by the Delta Method

We can derive an estimator for Var (*f* (R*tj*, G*tj*)) by computing the variance of the first-order Taylor's expansion of *f* (R*tj*, G*tj*) at ( ( *R G* ), ( )) *tj tj* :

$$\begin{split} \text{Var}\left(f(R\_{\eta},G\_{\eta})\right) &\approx \left(\frac{\partial f}{\partial R\_{\eta}}(\mathbb{E}(R\_{\eta}),\mathbb{E}(G\_{\eta}))\right)^{2} \text{Var}\left(R\_{\eta}\right) + \\ & \quad \left(\frac{\partial f}{\partial G\_{\eta}}(\mathbb{E}(R\_{\eta}),\mathbb{E}(G\_{\eta}))\right)^{2} \text{Var}\left(G\_{\eta}\right) + \\ & \quad 2\left(\frac{\partial f}{\partial R\_{\eta}}(\mathbb{E}(R\_{\eta}),\mathbb{E}(G\_{\eta}))\right) \left(\frac{\partial f}{\partial G\_{\eta}}(\mathbb{E}(R\_{\eta}),\mathbb{E}(G\_{\eta}))\right) \text{Cov}\left(R\_{\eta},G\_{\eta}\right). \end{split}$$

The second-order term was not considered because Var ( ) *Rtj* 2 and Var ( ) *Gtj* 2 cannot be usually estimated.

Since *M f tj tj tj* = − *R G R G tj tj* ( , ) l og ( ) log ( ) 2 2 , with the first and second derivative showed in Appendix 5, we can obtain an approximation of Var (*Mtj*) as follows:

$$\begin{split} \mathrm{Var}\left(M\_{\eta}\right) &\approx \left(\frac{1}{\ln\left(2\right)\mathrm{E}\left(R\_{\eta}\right)}\right)^{2} \mathrm{Var}\left(R\_{\eta}\right) + \left(-\frac{1}{\ln\left(2\right)\mathrm{E}\left(G\_{\eta}\right)}\right)^{2} \mathrm{Var}\left(G\_{\eta}\right) + 1 \\ &\quad 2\left(\frac{1}{\ln\left(2\right)\mathrm{E}\left(R\_{\eta}\right)}\right) \left(-\frac{1}{\ln\left(2\right)\mathrm{E}\left(G\_{\eta}\right)}\right) \mathrm{Cov}\left(R\_{\eta}, G\_{\eta}\right) \\ &= \frac{1}{\ln^{2}\left(2\right)} \left(\frac{\mathrm{Var}\left(R\_{\eta}\right)}{\mathrm{E}^{2}\left(R\_{\eta}\right)} + \frac{\mathrm{Var}\left(G\_{\eta}\right)}{\mathrm{E}^{2}\left(G\_{\eta}\right)} - 2\frac{\mathrm{Cov}\left(R\_{\eta}, G\_{\eta}\right)}{\mathrm{E}\left(R\_{\eta}\right)\mathrm{E}\left(G\_{\eta}\right)}\right). \end{split}$$

Consider the sample estimators of the expected values of *Rtj* and *Gtj*, denoted by, respectively, *Rtc* and *Gtc* , and assume that they are non-zero. Also, consider their variance and covariance sample estimators, denoted by, respectively, σˆ ( ) 2 *Rt* , σˆ ( ) 2 *Gt* , and ˆ( , ) σ *Rt G*t , and assume that they are independent of the background correction. We can derive the following estimator for Var (*Mtj*) :

$$
\hat{\boldsymbol{\sigma}}^{2}(M\_{\boldsymbol{t}}) \dot{=} \frac{1}{\ln^{2}(2)} \left( \frac{\hat{\boldsymbol{\sigma}}^{2}(R\_{\boldsymbol{t}})}{\overline{R}\_{\boldsymbol{t}\boldsymbol{c}}^{2}} + \frac{\hat{\boldsymbol{\sigma}}^{2}(G\_{\boldsymbol{t}})}{\overline{G}\_{\boldsymbol{t}\boldsymbol{c}}^{2}} - 2 \frac{\hat{\boldsymbol{\sigma}}(R\_{\boldsymbol{t}}, G\_{\boldsymbol{t}})}{\overline{R}\_{\boldsymbol{t}\boldsymbol{c}} \overline{G}\_{\boldsymbol{t}\boldsymbol{c}}} \right).
$$

Considering that *Atj* is defined by the function

$$f(R\_{\mathfrak{g}}, G\_{\mathfrak{g}}) \doteq \frac{\log\_2(R\_{\mathfrak{g}}) + \log\_2(G\_{\mathfrak{g}})}{2},$$

we can estimate Var (*Atj*) in a similar way to Var (*Mtj*).

By using the first and second derivatives of A*tj*, which are showed in Appendix (Barrett et al., 2012), we obtain the following approximation of Var (*Atj*):

$$\begin{split} \mathrm{Var}\,(A\_{\eta}) &\approx \left(\frac{1}{2\ln\left(2\right)\mathbb{E}(R\_{\eta})}\right)^{2} \mathrm{Var}\,(\mathcal{R}\_{\eta}) + \left(-\frac{1}{2\ln\left(2\right)\mathbb{E}(\mathcal{G}\_{\eta})}\right)^{2} \mathrm{Var}\,(\mathcal{G}\_{\eta}) + 1 \\ &\quad 2\left(\frac{1}{2\ln\left(2\right)\mathbb{E}(R\_{\eta})}\right) \left(-\frac{1}{2\ln\left(2\right)\mathbb{E}(\mathcal{G}\_{\eta})}\right) \mathrm{Cov}\,(\mathcal{R}\_{\eta}, \mathcal{G}t\_{\eta}) \\ &= \frac{1}{4\ln^{2}\left(2\right)} \left(\frac{\mathrm{Var}\,(\mathcal{R}\_{\eta})}{\mathbb{E}^{2}(\mathcal{R}\_{\eta})} + \frac{\mathrm{Var}\,(\mathcal{G}\_{\eta})}{\mathbb{E}^{2}(\mathcal{G}\_{\eta})} + 2\frac{\mathrm{Cov}\,(\mathcal{R}\_{\eta}, \mathcal{G}\_{\eta})}{\mathbb{E}(\mathcal{R}\_{\eta})\mathbb{E}(\mathcal{G}\_{\eta})}\right). \end{split}$$

Rewriting the above expression using the sample estimators for the expected value, variance and covariance of *Rtj* and *Gtj*, we derive the following estimator for Var (*Atj*) :

$$
\hat{\boldsymbol{\sigma}}^{2}(\boldsymbol{A}\_{\boldsymbol{t}}) \dot{\boldsymbol{\epsilon}} \doteq \frac{1}{4\ln^{2}(2)} \left( \frac{\hat{\boldsymbol{\sigma}}^{2}(\boldsymbol{R}\_{\boldsymbol{t}})}{\overline{\boldsymbol{R}}\_{\boldsymbol{u}}^{2}} + \frac{\hat{\boldsymbol{\sigma}}^{2}(\boldsymbol{G}\_{\boldsymbol{t}})}{\overline{\boldsymbol{G}}\_{\boldsymbol{t}\boldsymbol{c}}^{2}} + 2\frac{\hat{\boldsymbol{\sigma}}\left(\boldsymbol{R}\_{\boldsymbol{t}},\boldsymbol{G}\_{\boldsymbol{t}}\right)}{\overline{\boldsymbol{R}}\_{\boldsymbol{t}\boldsymbol{c}}\overline{\boldsymbol{G}}\_{\boldsymbol{u}\boldsymbol{c}}} \right).
$$

# The Construction and Comprehensive Analysis of ceRNA Networks and Tumor-Infiltrating Immune Cells in Bone Metastatic Melanoma

#### *Edited by:*

*Liang Cheng, Harbin Medical University, China*

#### *Reviewed by:*

*Chuan-xing Li, Karolinska Institute (KI), Sweden Dapeng Hao, Baylor College of Medicine, United States*

#### *\*Correspondence:*

*Jie Zhang<EMAIL_ADDRESS>Tong Meng<EMAIL_ADDRESS>Zongqiang Huang<EMAIL_ADDRESS>
*†These authors have contributed equally to this work*

#### *Specialty section:*

*This article was submitted to Statistical Genetics and Methodology, a section of the journal Frontiers in Genetics*

> *Received: 26 June 2019 Accepted: 12 August 2019 Published: 25 September 2019*

#### *Citation:*

*Huang R, Zeng Z, Li G, Song D, Yan P, Yin H, Hu P, Zhu X, Chang R, Zhang X, Zhang J, Meng T and Huang Z (2019) The Construction and Comprehensive Analysis of ceRNA Networks and Tumor-Infiltrating Immune Cells in Bone Metastatic Melanoma. Front. Genet. 10:828. doi: 10.3389/fgene.2019.00828*

*Tong Meng2,3,4\* and Zongqiang Huang1\* 1 Department of Orthopaedics, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China, 2 Division of Spine, Department of Orthopedics, Tongji Hospital affiliated to Tongji University School of Medicine, Shanghai, China, 3 Tongji University School of Medicine, Tongji University, Shanghai, China, 4 Department of Orthopedics, Shanghai General Hospital, School of Medicine, Shanghai Jiaotong University, Shanghai, China, 5 Shanghai East Hospital, Key Laboratory of*

*Arrhythmias, Ministry of Education, Tongji University School of Medicine, Shanghai, China*

*Runzhi Huang1,2,3†, Zhiwei Zeng1†, Guangyu Li1†, Dianwen Song4, Penghui Yan1, Huabin Yin4, Peng Hu1, Xiaolong Zhu1, Ruizhi Chang1, Xu Zhang1, Jie Zhang5\*,*