extfieldsofstudy:
Biology
Computer Science
Medicine

text:
Integrative network analysis of TCGA data for ovarian cancer

Background Over the past years, tremendous efforts have been made to elucidate the molecular basis of the initiation and progression of ovarian cancer. However, most existing studies have been focused on individual genes or a single type of data, which may lack the power to detect the complex mechanisms of cancer formation by overlooking the interactions of different genetic and epigenetic factors. Results We propose an integrative framework to identify genetic and epigenetic features related to ovarian cancer and to quantify the causal relationships among these features using a probabilistic graphical model based on the Cancer Genome Atlas (TCGA) data. In the feature selection, we first defined a set of seed genes by including 48 candidate tumor suppressors or oncogenes and an additional 20 ovarian cancer related genes reported in the literature. The seed genes were then fed into a stepwise correlation-based selector to identify 271 additional features including 177 genes, 82 copy number variation sites, 11 methylation sites and 1 somatic mutation (at gene TP53). We built a Bayesian network model with a logit link function to quantify the causal relationships among these features and discovered a set of 13 hub genes including ARID1A, C19orf53, CSKN2A1 and COL5A2. The directed graph revealed many potential genetic pathways, some of which confirmed the existing results in the literature. Clustering analysis further suggested four gene clusters, three of which correspond to well-defined cellular processes including cell division, tumor invasion and mitochondrial system. In addition, two genes related to glycoprotein synthesis, PSG11 and GALNT10, were found highly predictive for the overall survival time of ovarian cancer patients. Conclusions The proposed framework is effective in identifying possible important genetic and epigenetic features that are related to complex cancer diseases. The constructed Bayesian network has identified some new genetic/epigenetic pathways, which may shed new light into the molecular mechanisms of ovarian cancer. Electronic supplementary material The online version of this article (doi:10.1186/s12918-014-0136-9) contains supplementary material, which is available to authorized users.

the cancer status. Thus a systems biology approach that combines multiple genetic and epigenetic profiles for an integrative analysis provides a new direction to study the regulatory network associated with ovarian cancer.
The rapid advances in next-generation sequencing technology now allow genome-wide analysis of genetic and epigenetic features simultaneously. The timely advent of TCGA project has provided the most comprehensive genomic data resource from over 20 types of cancers (http://cancergenome.nih.gov/). For example, the TCGA ovarian cancer data contain both clinical and molecular profiles from 572 tumor samples and 8 normal controls. The molecular profile includes gene expression (microarray), genotype (SNP), exon expression, MicroRNA expression (microarray), copy number variation (CNV), DNA methylation, somatic mutation, gene expression (RNAseq), MicroRNA-seq and protein expression. The clinical information includes records on recurrence, survival, and treatment resistance. These massive complex data sets have driven enthusiasm to study the molecular mechanism of cancers through computational approaches [1,[6][7][8]. Among the developed methods, Bayesian Network (BN) is one of the most frequently used multivariate models. The BN approach is more appealing than graphs constructed based on correlation or mutual information metrics for it allows rigorous statistical inference of causality between genetic and epigenetic features. However most of the existing studies have been focused on one type of data either continuous or discrete [9][10][11][12][13]. How to combine different types of complex data for causal inference in BN poses a big challenge. In addition, deducing the complex network structure from data remains an open problem partially due to the lack of prior information, relatively smaller sample size and the high dimensionality of data (number of possible nodes) [13,14].
A necessary and important step to construct a BN from tens of thousands of features is feature selection, i.e., to identify a subset of the most-relevant features. Removing irrelevant or redundant features helps improve computing efficiency and estimation accuracy in the causal network. Existing feature selection methods can be roughly classified into two categories: wrapper approach [15,16] and filter approach [17][18][19]. For large data sets, the filter approach using significance test for difference between the cancer and control samples is more commonly used due to its simplicity. As some features could be causal to other features while having no direct association with the cancer phenotypes, the independent test can filter out many related features (see a simulation study in the Methods section). One innovation of this paper is a novel stepwise correlation-based selector (SCBS) that mimics the hierarchy of the BN for feature selection. The selected features from the TCGA data are a mixture of continuous and categorical variables. To integrate them into the same BN, we discretize the continuous variables and use a logit link function for casual inference. The proposed approach is applied to the TCGA ovarian cancer data and leads to a series of interesting findings that shed light into the genetic/epigenetic mechanisms of ovarian cancer.

Preprocessing of TCGA ovarian cancer data
In this paper, we only consider four types of molecular data including gene expression, DNA copy number variation, promoter methylation and somatic mutation (summarized in Table 1). This data set contains the expression values of 17,812 genes, out of which, 12,831 had methylation level measured for each CpG island located in their promoter regions. If multiple CpG islands exist for a given gene, we took the average as the overall methylation level. The copy number was measured for each chromosomal segment, recorded as a seg.mean value, with the segment length varying from hundred up to tens of millions base pairs by the Circular Binary Segmentation (CBS) performed by TCGA. Out of 17,812 genes, 15,352 have well defined locations on the genome provided by UCSC Genome Browser (http://genome.ucsc.edu/) and each of them was assigned a value as the copy number. If a gene entirely falls within a chromosome segment, we assigned it the corresponding seg.mean value (236 out of 15,352 genes that spans two chromosomal segments were not assigned any seg.mean value). For the somatic mutation data, we defined a binary variable where "1" stands for all the non-silent mutations (coding different amino acid) and "0" for silent mutations or not being mutated, resulting in 9,895 genes with somatic mutations (mutation occurred in at least one sample). For the gene expression and methylation data, we applied an existing method by [8] to remove the effects due to different age groups and batches. Figure 1 illustrates an example for BRCA1 gene where the boxplots showed this procedure was effective in removing the age and batch effects.
The second step in data preprocessing is to discretize the continuous random variables. We classified the gene expression level into three groups using k-means clustering algorithm, namely low, midium, and high. Likewise the promoter methylation level was classified  Boxplots of BRCA1 expression before (on the left) and after (on the right) removal of age effect, x-axis is age group (<50 yrs old, 50-70 yrs old and >70 yrs old) and y-axis is the expression level. In the preprocessing step, we removed batch and age effects of expression level and methylation level for every single gene.
into either hyper or hypo state, and the copy number status into two states: gain or loss. The discretized variables were the input of the Bayesian network to be discussed below.

Feature selection
The pipeline we propose in this paper assumes that cancer phenotype is directly associated with gene expression, which can be potentially driven by genetic and epigenetic changes (this assumption was made from the biological point of view. It can be dropped without affecting the modeling and computing). Figure 2 illustrates the workflow of our proposed framework. We first identify a set of tumor suppressors and oncogenes by differential expression analysis between the cancer and control groups. This set of genes, together with those that were previously reported in the literature as cancer-relevant, form the set of seed genes for further feature selection. These seed genes are then fed into our proposed stepwise correlation-based selector (SCBS) to select other features. The SCBS is motivated by the hierarchy of causality in the Bayesian network. For example, suppose there is a causal relationship A→B→cancer. Though A to B or B to cancer has strong directed association, the association between A and cancer could decay greatly so that it cannot be detected. The SCBS procedure starts with detection of features significantly correlated with the cancer and then progressively selects subsequent features that correlate with the selected features. Our simulation study presented in the Methods section shows that the proposed SCBS is more effective to select important features that are involved in the phenotype-related pathways but indirectly associated with the cancer phenotype. The first step in feature selection is to define a set of seed genes out of 12,000 genes that have records of expression level and at least one of the three (epi)genetic factors (i.e., CNV, methylation and mutation). We conducted a nonparametric test to identify the most differentially expressed genes between case and control groups, as well as the most differentially methylated genes, most differentially mutated genes and genes with most differentiate copy numbers. A gene was defined as a candidate oncogene or tumor suppressor if it satisfies all the following three conditions: (1) the p-value of gene expression was statistically significant under Benjamini-Hochberg (BH) procedure with false discovery rate (FDR) ≤ 0.05; (2) the p-value of somatic mutation/ promoter methylation/copy number was also significant under BH procedure with FDR ≤ 0.05; and (3) the absolute value of correlation coefficient between gene expression and somatic mutation/promoter methylation/copy number was greater than 0.4. The p-values were calculated using Wilcoxon rank-sum test and the correlation coefficients were calculated using Pearson's method. This procedure resulted in 48 potential tumor suppressors or oncogenes ( Table 2), 16 of which were well-studied tumor suppressors and oncogenes for ovarian cancer from literature (out of a total of 36 presented in [2], Table 3). The union of the two sets of 68 genes were defined as seed genes and generated an additional 271 nodes out of more than 50,000 candidate features by the stepwise correlation based selection (SCBS) procedure, which include 177 gene expressions, 82 copy number variation sites, 11 methylation sites and one somatic mutation site at gene TP53.

Bayesian network prediction
The 339 nodes (discretized if continuous) were fit into the Bayesian network through a logit link function using the blockwise coordinate descent algorithm for penalized maximum likelihood estimation procedure [20] (Methods).The predicted network contains 698 edges ( Figure 3, details are tabulated in Additional file 1: Table  S1), where the direction of the edge indicates the downstream feature is regulated by the upstream one. We found the CNVs are the major factor that accounts for differential gene expression. In addition most of the 82 genes were CNV-amplified in cancer samples, suggesting that many amplified genes may act as cancer drivers, confirming findings from a breast cancer study [21]. This network also confirmed many previously reported gene-gene interactions. To name a few, for example, the edge from TPX2 to AURKA could be explained by the fact that the protein encoded by the TPX2 gene activates AURKA by inducing autophosphorylation [22]. The connection between BRCA1 and NBR2 could be due to the shared bi-directional promoter between the two genes [23]. The connection between STAT3 and ETV6 was suggested previously that ETV6 is a negative regulator of STAT3 activity [24]. The edge from CDKN2A to CCNE1 is a known gene-gene regulation in the RB signaling pathway [1]. The edge from MYC to IMPHD2 confirms that MYC depletion results in repression of IMPHD2 (a gene coding rate-limiting enzyme) [25]. These results suggested that  the proposed pipeline is capable of revealing important genetic or epigenetic pathways that underlie the complex cancer phenotype. The average degree (indegree plus outdegree) of the graph is 4.124 indicating that the inferred network is sparse with moderate complexity. Due to the directionality of BN, one could also obtain the indegree and outdegree separately for every single node. Figure 4 shows the observed distribution of outdegree, where the mean and standard deviation are 2.15 and 2.31 respectively. We identified 13 nodes with significantly larger outdegree (greater than mean + 2 × SD, ≥ 7 edges) in the network including ARID1A, C19orf53, CSNK2A1, DERL1, TRMT6, COL5A2, TCF21, LUM, TPX2, UBE2C,  Table S3. DPM1, NDUFB7, and NDUFB9 (Table 4). These hub genes all have known functions and have causal effect on at least seven neighboring genes, suggesting that they may play important roles in driving corresponding local subnetworks. Some of the hub genes have been reported in the literature that are related to ovarian cancer. For instance, ARID1A is known to promote the formation of SWI/SNF chromatin remodeling complexes containing BRG1 or BRM, and is a candidate tumor suppressor not only in clear cell ovarian cancer but also in endometrioid cancers and uterine endometrioid carcinomas [22,26,27]. C19orf53 is known to be associated with Leydig cell tumors which are a member of the sex cord-stromal tumor group of ovarian and testicular cancers and it has a potential role in hypercalcemia of malignancy [22]. CSKN2A1 is a well-known oncogene that can phosphorylate a number of key intracellular signaling proteins implicated in tumor suppression (P53 and PTEN) and oncogenesis (MYC, JUN, NF-kappaB). This gene also influences Wnt signaling via beta-catenin phosphorylation and the PI3K signaling pathway via the phosphorylation of AKT [22]. Interestingly these 13 hub genes can clearly distinguish the cancer samples from the normal samples as revealed by a multi-dimensional scaling plot (MDS, [28]) based on the correlation dissimilarity metric (comparable clustering effect was observed based on the entire set of 245 genes, Figure 5a,b). This suggests that the thirteen hub genes may present the major difference between the cancer and normal samples. The early-stage and high-grade tumor samples however are not well distinguished.

Gene clusters
The 245 genes (listed in Additional file 2: Table S2) were identified to fall into four major clusters corresponding to distinct functions by k-means clustering method (Figure 6b). Cluster 1 (black in Figure 7a) contains 18 genes, mainly related to cell division, mitosis, spindle formation etc. Cluster 2 (red) contains 23 genes, most of which are functionally related to growth factor, cell shape, cell motility, tumor invasion etc. Cluster 3 (green)  Table 4.
contains 20 genes, mostly related to mitochondrial system, membrane process etc. Cluster 4 (blue in Figure 7a) is the largest and most complicated cluster harboring 184 genes. This large cluster communicates between the other three clusters (as shown in Figure 7b), which are nearly independent from each other. This can be seen from the summary statistics of within and between cluster causal edges in Table 5. These findings could be implicative of some important molecular pathways, which may or may not have been identified, that drive the development of ovarian cancer. We also looked into the subnetwork within each cluster. Figure 8 shows the local subnetwork corresponding to the first cluster which is involved mainly in cell division processes. Two hub genes, TPX2 and UBE2C, are in the central positions of this network. Our finding that there are four gene clusters may suggest that the development of ovarian cancer could be partially driven by the cell cycle regulation, as well as the pathways related to cell shape and motility, and mitochondrial system. Intervention (activation or suppression) on the hub genes or other important genes may alter the entire network, therefore may control key aspects of disease development.

Survival-centric network prediction
We applied SCBS and Bayesian network model to construct a survival-centric network. In feature selection, we used continuous data to select genetic/epigenetic features that were most relevant to the overall survival time of ovarian cancer patients. Starting from the first node (overall survival time in days), the SCBS procedure selected 66 genes ( Table 6) and two methylation sites as the nodes in the Bayesian Network. Interestingly only 6 genes including CCDC19, MMP1, SLC family, TEKT2, WDR family and ZMYND10 had been reported relevant to cancer survival in a separate study (where a total of 88 genes were reported, [29]). The death risk (binary) is then used as the phenotype node to build a phenotypeinduced network. We defined the overall survival time for less than 2 years as the "high-risk" (134 patients) and the survival time greater than 4 years as the "lowrisk" (101 patients). In our predicted graph, there are 9 hub genes (similarly defined as above) that drive corresponding local subnetworks including C2orf39, FAP, SLC2A2, LAPTM5, CD53, THBS2, CCDC63, SLC17A, and LCT. Gene FAP has a known function to control fibroblast growth or epithelial-mesenchymal interactions during development, tissue repair, and epithelial carcinogenesis [22]. The inferred Bayesian network (Figure 9, details are tabulated in Additional file 3: Table  S3) identified two genes, namely PSG11 and GALNT10, that may be directly associated with the overall survival time of ovarian cancer patients ( Figure 10). Both genes are functionally related to glycoprotein synthesis, as well as many other genes in the network such as SLC2A2, SLC17A, CD53, THBS2, LCT, GYPA [22]. This indicates the biological pathway related to glycoprotein synthesis may be implicative of death risk of ovarian cancer patients. As reported in literature [30,31] several tumor-associated glycoproteins were found on the surface of many cancer cells including ovary, breast, colon, and pancreatic cells and they may play potential roles in early detection of caners. One well-known such glycoprotein is CA-125 (encoded by gene MUC16), which is the primary protein used to measure serous cancer tumor load, especially during recurrence, and it is heavily glycosylated [31]. Pregnancy-specific glycoproteins (PSG) are mainly produced by the placental syncytiotrophoblasts during pregnancy and these proteins comprise a subgroup of the carcinoembryonic antigen family [22]. The protein encoded by gene GALNT10 may have increased catalytic activity toward glycosylated peptides compared to activity toward non-glycosylated peptides [22]. As pointed out by several research groups [32][33][34], some certain glycoproteins are closely associated with women cancers such as ovarian cancer and breast cancer, affecting the death risk, chemotherapy resistance and prognosis of ovarian cancer patients. The network also involves genes of other important functions including microtubules (RSHL3, TEKT1 and OLR1), extracelluar processes (ECM1, THBS2), hematopoietic (MS4A4, SRGN, LAPTM5), and human immune system (LILRA1, SIGLEC7, LAT2, LAIR1). This network suggested many causal relationship between different features, some of which have been known, for example, the edge HOTAIR→ HOXC10 may be due to the fact that HOTAIR (a noncoding RNA gene) is located within the Homeobox C (HOXC)

MDS plot for 245 cancer related genes
First direction gene cluster and it regulates the expression of HOX genes such as HOXC and HOXD [22,35].

Discussion
In this paper we proposed an integrative approach in the Bayesian network framework for causal inference between genetic and epigenetic features in complex cancer data. It presents novelty in two aspects. First, we demonstrated that the stepwise correlation-based selection approach is more effective than simple single-round selection method in identifying important features in the genetic/epigenetic pathways, particularly those that are indirectly associated with the underlying phenotype. The method we proposed relies on the correlation strength among connected nodes and may fail when the connections are weak, especially for high-dimension data where assigning significance is challenging. Using literature-verified genes as seeds helps to better select features relevant to the phenotype. The SCBS procedure is model-free and computationally efficient and it can be applied to other graphical models such as Markov Random Field (MRF, undirected graph) and gene-gene or protein-protein interaction (PPI) network problems. Second, we proposed a method for modeling causal relationships between features of different types (continuous or discrete) in a Bayesian network through a logit link function. The block-wise coordinate descent (BCD) algorithm accompanied with the Bayesian network model provides a simple and efficient way to estimate the parameters in the model. With a moderate sample size, this method achieves reasonable accuracy even for a moderate-scale network containing 200 nodes. This pipeline can be readily applied to other complex cancer data for pathway analysis or to find the common pathways between different but correlated phenotypes. For example, the TCGA now has accumulated more than 1,000 samples from breast cancer patients. It would be highly intriguing to know whether these two common diseases share any common molecular basis, especially for the basal-like breast cancer. As previously reported in literature [36], the basal-like subtype of breast cancer is the most distinct of the four subtypes (Luminal A, Luminal B, HER2-enriched and Basal-like) and it is similar to high-grade serous ovarian cancer at the mRNA expression level. The Bayesian network model allows strict probabilistic inference, but also has limitations. First, it lacks flexibility to model cyclic causal relationships. The BN cannot model any cyclic pathways, for example, A→B→C→A, which though may exist in gene regulation in vivo. Second, most of the existing BN learning algorithms assume sparsity for computational feasibility. If the true network is dense or locally dense (dense subnetworks), the weak causations may fail to be detected. Our proposed BN with logit link function can model categorical variables and discretized continuous variables simultaneously, while discretization procedure may cause loss of information. A more desirable way could be to model the discrete variable nodes with logit distribution, but the continuous variable nodes with Gaussian or other continuous distributions. This direction demands future research. The coordinate descent is successful in solving the LASSO-type problem, especially in the sparse BN problem. Due to the super-exponentially increased number of possible BNs, the traditional sampling-based methods such as Metropolis-Hastings algorithm ( [13,14]) is computational infeasible to estimate the network with moderate number of nodes, e.g., 100 nodes. In our problem involving 339 nodes, the BCD algorithm took about 20 minutes on a single CPU to complete the estimation for a given penalizing parameter λ (Materials and Methods).
We illustrated the flexibility of this pipeline with two phenotypes, namely the ovarian cancer phenotype itself and the survival time of the cancer patient. The induced network by the cancer phenotype contains a set of 245 genes forming into four major clusters of distinct major functions, coordinated by 13 hub genes. Some of the hub genes (e.g., ARID1A) have been reported by other researchers for their important roles in genetic pathways, while other genes (NDUFB7 and NDUFB9) are newly discovered in this study. Their functions in ovarian cancer need to be further investigated. Our discovery that the mitochondrial systems are regulated in serous tumors is consistent with the hypothesis that the Warburg effect impacts tumor progression as suggested in published studies ( [37]). We also found that pathways related to glycoprotein synthesis, hematopoietic and immune systems correlate with the survival rate of ovarian cancer patients. In particular, we discovered that the two genes related to glycoprotein synthesis, PSG11 and GALNT10, can significantly affect the overall survival time of ovarian cancer patients.

Conclusions
Understanding the biological mechanism of ovarian cancer has significant practical importance for clinical   Table S3.
diagnosis and treatment. Timely advent of TCGA project has provided the most comprehensive genomic data resource for cancer study at the molecular and system level. Nevertheless, how to utilize these complex data for discovery of molecular mechanism of cancers remains one of the biggest challenges in this field. To this end, we propose a new integrative approach in this paper, which presents two innovations: a stepwise feature selection procedure and a Bayesian network model that incorporates both continuous and discrete features for causal a b inference. The predicted graph for the ovarian cancer data confirmed numerous genetic pathways reported in the literature, as well as many new ones that may provide new clues to guide future research. The graph suggested 13 hub genes that may drive certain subnetworks therefore play important roles in ovarian cancer. Clustering analysis suggested four gene clusters corresponding to distinct biological processes including cell division, tumor invasion and mitochondrial system. In addition, we found that genes related to glycoprotein synthesis, hematopoietic, immune system could be highly predictive of overall survival time of ovarian cancer patients.

Data
Datasets in ovarian cancer were downloaded from the Cancer Genome Atlas (TCGA) data portal (http://tcgadata.nci.nih.gov). We extracted four types of molecular data including gene expression, promoter methylation, somatic mutation and DNA copy number variation, using the "data matrix" tool provided by TCGA data portal.

Bayesian network with a logit link function
Bayesian Network can be used to model a set of random variables (nodes) and their conditional dependencies (directed edges) [11]. In general, the joint likelihood function of nodes X 1 , . . . , X p in a BN can be expressed as: where graph G = (V , E) represents the topological structure of the Bayesian network, V = X 1 , . . . , X p denotes the set of nodes and E = X j → X i , X j ∈ G i denotes the set of edges, and G i ⊆ X 1 , . . . , X p \ {X i } stands for the parent set of X i . We say X j ∈ G i if X j causes X i , written as X j → X i . A BN is called Gaussian Bayesian Network (GBN) if X i is normally distributed with the mean equal to a linear combination of X j ∈ G i . The GBN is the most popular BN model and its structure learning problem has been discussed by several researchers [20,38]. In our motivating example, however the network involves both continuous and discrete random variables. Here we discretize the continuous random variable and consider a multinomial logistic model. Let X i take values from {1, . . . , K i } with probabilities π i1 , . . . , π iK i s.t.
where k = 1, . . . , K i − 1 and i = 1, . . . , p. Here we transform the network structure to a coefficient matrix where β ikjl = 0 for all k and l means X j → X i , and otherwise X j → X i . Therefore estimating the structure of G is equivalent to estimating matrix {β ikjl }. For simplicity, we illustrate the parameter estimation using binomial logistic model where all X i 's only take values 0 or 1. Define π i ≡ prob(X i = 1), then Suppose we observe data from N subjects. Let X n = (X n0 , X n1 , . . . , X np ) T , where X n0 = 1 is the dummy variable, and X ni = 0 or 1 for n = 1, . . . , To achieve the sparsity, we apply the L 1 penalty to the loglikelihood ( [20,38]) as follows: We aim to optimize the objective function (5) under the constraint of acyclicity. Finding the global maximizer is typically difficult in such a high dimensional space. Here we consider the coordinate descent (CD) algorithm, which has been successfully used to solve lasso regression problems [20,39]. The CD algorithm is based on single-parameter updating strategy to minimize the objective function coordinate-by-coordinate. For our model, the single-parameter updating can be done as follows, in particular, we seek the maximizerβ ij of the following objective function given all the other parameters, denoted β −ij : After excluding the constant part, we have: log 1 + C 2n exp C 3n β ij , do the following to find the maximizer of L * i β ij |β −ij : The acyclicity constraint brings a major difficulty in BN learning problem, especially when the topological order of nodes is unknown. One immediate result of this constraint is that β ij and β ji cannot be both nonzero. To take advantage of this implication, Fu and Zhou (2013) proposed a blockwise coordinate decent (BCD) algorithm where the p(p − 1) parameters are partitioned in to p(p − 1)/2 blocks. Each block consists of β ij and β ji . The BCD algorithm [20] can be implemented as follows (starting from an empty network where all β ij =0): • Step 1: For each pair of β ij and β ji ,β ij ⇐ 0 stands for thatβ ij has to be 0 under the constraint of acyclicity: a. Ifβ ji ⇐ 0, find the maximizerβ ij of L * i w.r.t. β ij . b. Ifβ ij ⇐ 0, find the maximizerβ ji of L * j w.r.t. β ji . c. If eitherβ ij orβ ji can be nonzero, then compare the two sums: . Find maximizer of max(S 1 , S 2 ).
• Step 2: Repeat step 1 until the maximum difference between two successive cycles is below some threshold.
To check the acyclicity of the candidate graph after edge i −→ j is added, we use a simple breadth-first algorithm detailed as follows. The time complexity of this algorithm is O(|V | + |E|): • Step 1: Remove all the edges coming into j and identify the children set of j, denoted by CS j .
• Step 2: If i ∈ CS j , then stop. Otherwise remove all the edges coming into CS j and find the children set of CS j . • Step 3: Repeat until i is found or all the edges were removed.
Notice that in multinomial setting, the acyclicity constraint forces multiple β's to be zero. For instance, in ternary case, X i → X j indicates that four β's are zeroes simultaneously. Another important issue in Lasso regression is the choice of the tuning parameter λ. Cross-validation is the most commonly used method for selecting λ, which however, tends to select a too small λ resulting in high false positive rate ( [20]). To overcome this difficulty, an empirical λ selection method was proposed in [20] that can guarantee significant increase of the maximized likelihood value as a function of the graph complexity (number of edges). We employed this method for the selection of λ and the sequence of candidates is set to be {1/16, 1/8, 1/4, 1/2, 1, 2, 4, 8, 16}.

Simulation I: BCD algorithm for BN with discrete nodes
In the first simulation, we evaluated the performance of BCD algorithm in the proposed Bayesian network with logit link function as follows. We first simulated a random graph (i.e., G = (V , E)) with p nodes and 2p edges respectively. The simulated graph is a weakly connected directed graph generated using R package bnlearn ( [40]). Given each G, we simulated N independent observations, i.e., X n = (X n1 , . . . , X np ), for n = 1, .., N. If X nj is caused by other nodes, then we simulated it based on Bernoulli distribution with success probability following the logistic regression model. For simplicity, for each observation the causal effect β ij (if X i → X j in G) in the logistic regression model was set as a constant. For those nodes in the network but not caused by any other nodes, we simulated them from binomial distribution independently with success probability randomly generated from uniform distribution from 0.1 to 0.9. In this simulation we consider a crossed design of p = 50, 100, 200, β ij = 0.5, 1 and N = 500, 1000, 2000. For each setting of (p, β ij , N), 10 replicated samples were generated. We evaluated the estimated network under two different criteria: the directed version and undirected version of network (skeleton). In the former, we count an edge as a true positive only if it has the correct link and direction. For the skeleton comparison, an edge is counted as a true positive as long as it has the correct link. Table 7 presents the average number of predicted edges (P) from the 10 replicates for each setting, true positive rate (TPR) and false discovery rate (FDR) for both directed and undirected edges. Unsurprisingly the estimation of network structure is affected by the complexity of network, magnitude of causal effect and sample size. For a sparse network with 200 nodes, and sample size N = 2000, and β = 1, the BCD algorithm achieves an average TPR of 0.7 and 0.96 for directed and undirected (skeleton) graphs respectively. This simulation demonstrates that the BCD algorithm performs reasonably well when applied to categorical data in a moderate complex network when the sample size is relatively large.

Comparison of three BN models on real data
To benchmark our logistic BN model, we compared it with two other BN models, namely the Gaussian BN model and multinomial BN model, on a popular data set [41] where the true causal network is known and experimentally-validated. This data set contains the abundance measurement of 11 proteins in 5400 samples, and has been used to elucidate the signaling pathway structure. Both the continuous and discrete versions of data are available online. The known protein-protein network is a Bayesian Network containing 11 nodes and 20 directed edges. Because this data set is based on experimental interventions, we adapt our model by deleting the intervention terms from the likelihood function. Figure 11 shows the true graph and estimated graphs by three different models and Table 8 summarizes the true positive rate and false discovery rate by three models. In terms of TPR and FDR, the logistic a b c d  model appears to perform slightly better than other two models.

Simulation II: stepwise correlation-based feature selector
In the second simulation, we compare the proposed stepwise correlation-based feature selector with other existing methods. The feature selection step is to define a set of nodes of smallest possible size but include most possible nodes that are involved in the casual structure in the true phenotype-induced Bayesian network. Suppose we aim to select p variables from S (S >> p) candidates as the nodes in Bayesian network based on a random sample of N observations. Assume that the phenotype-induced Bayesian network truly has p + 1 nodes (including the phenotype itself ). The proposed stepwise feature selection method starts with the phenotype node and selects the features that are most correlated with the current nodes in a stepwise fashion based on a correlation or mutual information metric (in this paper, we use correlation). This procedure is a natural mimic of network structure and can identify those nodes indirectly associated with the phenotype. In practice, the method can be implemented as follows: • Step 1: Calculate the correlation coefficients between the current node X i and all the other nodes, denoted by ρ ij , j = i. Keep k most correlated nodes with X i based on ρ ij for further filtering. • Step 2: Calculate the p-value of correlation coefficient for each of the k nodes from step 1, select the node if the p-value is significant under Benjamini-Hochberg (BH) procedure with FDR≤ 0.05. • Step 3: Repeat step 1 and 2 until p nodes are selected.
In practice we need to pre-define the value of p and k based on the complexity of the network. The choice of p is subject to the feature pools size S and the scale of the network to build. The computing time is sub linear to p. We recommend to choose a k of 4, 5 or 6 to attain moderate complexity or sparsity of the network (see a simulation study below for the choice of k).
To test the SCBS method, we conducted a simulation study with S = 10, 000 features, among which only p = 49 are truly involved in the phenotype related network. We first generated 50 random sparse graphs (i.e., G = (V , E)) consisting of 50 nodes and 100 directed edges (one node will be randomly chosen to be the phenotype node) using bnlearn ([40]). For each graph, we simulated eight samples according to the binomial logistic model, four with constant β ij = 1 and the rest with β ij = 2 at four different sample sizes N = 500, 1000, 2000, 5000. By using the topological order, each node can be simulated conditioning on the outcome of its parent nodes. For those nodes in the directed network but not caused by any other nodes, or the 9950 features outside the network, we simulated them from binomial distribution independently with success probability randomly generated from uniform distribution from 0.1 to 0.9. The SCBS with k = 4 was applied to each data set and 49 features are identified from the sea of 10,000 candidates. Compared to the Pearson's Chi-square test (single-round test for features between two phenotypic groups), the proposed SCBS performs uniformly better in all situations in terms of true positive rate (as shown in Figure 12, total positives were controlled at 49 for both methods). In particular, a small β poses more challenges in estimation or testing, while the SCBS can outperform the Chi-square test by 3-4 folds in true positive rate of identified phenotype-related features.
A further simulation study on the choice of k was carried out using the same strategy as above with S = 10, 000, p = 49, β ij = 1 and N = 1000. We generated 50 random graphs with average degree 4 (100 edges), 50 graphs with  average degree 8 (200 edges) and 50 graphs with average degree 12 (300 edges), then apply SCBS with different k (k = 1, 2, . . . , 10) to each data set and select 49 features. Figure 13 shows the true positive rates of the selected features under different choices of k for networks with different average degrees. It is shown that increasing network density leads to increasing true positive rates, and a k of 4, 5 or 6 performs better in all situations. Therefore for networks with moderate complexity, we recommend to choose k = 4, 5, 6 since smaller k tends to miss weakly connected nodes and larger k tends to catch more false positives.


extfieldsofstudy: Biology Medicine

text:
Protein-protein Interaction Network Prediction by Using Rigid-Body Docking Tools: Application to Bacterial Chemotaxis

Core elements of cell regulation are made up of protein-protein interaction (PPI) networks. However, many parts of the cell regulatory systems include unknown PPIs. To approach this problem, we have developed a computational method of high-throughput PPI network prediction based on all-to-all rigid-body docking of protein tertiary structures. The prediction system accepts a set of data comprising protein tertiary structures as input and generates a list of possible interacting pairs from all the combinations as output. A crucial advantage of this docking based method is in providing predictions of protein pairs that increases our understanding of biological pathways by analyzing the structures of candidate complex structures, which gives insight into novel interaction mechanisms. Although such exhaustive docking calculation requires massive computational resources, recent advancements in the computational sciences have made such large-scale calculations feasible. different rigid-body docking tools with different scoring models. We found that the predicted interactions were different between the results from the two tools. When the positive predictions from both of the docking tools were combined, all the core signaling interactions were correctly predicted with the exception of interactions activated by protein phosphorylation. Large-scale PPI prediction using tertiary structures is an effective approach that has a wide range of potential applications. This method is especially useful for identifying novel PPIs of new pathways that control cellular behavior.

INTRODUCTION
Identifying protein-protein interaction (PPI) networks is useful in understanding molecular biological phenomena such as signal transduction and the regulation of gene expression [1]. The most powerful methods used for PPI network prediction includes sequence based searching [2][3][4] and evolutional information based methods [5]. In addition, there are approaches that also employ protein tertiary structure to investigate PPIs. Ogmen et al. developed a PPI prediction method based on the structural alignments of query proteins to known interface surfaces in the database [6]. Acuner Ozbabacan et al. applied this method to predict PPI in the human apoptosis pathway and identified several candidates that were not listed in the corresponding database [7]. These methods work well for most of the cases in which the interaction databases have similar sequences or interaction surface structures to those of the query proteins. However, there are examples, such as the epidermal growth *Address correspondence to this author at the Graduate School of Information Science and Engineering, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8550, Japan; ‡ Present address: Department of Physics, Chuo University, 1-13-27, Kasuga, Bunkyo-ku, Tokyo 112-8551 Japan; Tel: +81 3-5734-3645; Fax: +813-5734-3646; E-mail: akiyama@cs.titech.ac.jp † Yuri Matsuzaki and Masahito Ohue has equally contributed to this work. factor receptor (EGFR) signaling pathway, where a complete set of interactions is unknown. In such cases it is difficult to predict novel PPIs using existing methods because they depend on known interaction regions. In order to discover novel interactions, the candidate interaction sites are searched throughout the entire surface of each protein structure.
Identifying large-scale PPI networks based on the full tertiary structure is computationally expensive. However, recent progress in computer science has facilitated the development of supercomputers that are able to perform such large-scale analyses. In the PPI research field, some pioneering studies have been done by taking advantage of this huge computational power. Mosca et al. presented a docking prediction database of the yeast interactome by using existing and modeled structures [8]. This idea was also later employed in the field of structural systems biology [9].
Recent studies have demonstrated novel methods of PPI network prediction by employing exhaustive rigid-body protein-protein docking and post-docking analysis [10][11][12]. Protein-protein docking is usually used to predict possible binding conformation between two proteins that are already known to undergo direct interaction. However, these docking calculations have now been extended to include all combinations of proteins, including protein pairs that were not previ-ously known to directly interact with one another. In these works, the affinities of given protein pairs are evaluated using groups of models with high interaction energy scores. In our previous study, we examined up to ten thousand possible structure pairs and showed prediction performance that was better than random [12]. The same procedure was subsequently applied to RNA molecules [13]. A similar procedure has been used by Wass et al. with a larger dataset (about 50,000 pairs), using a different rigid-body docking tool Hex with shape complementarity score function and post-docking analysis [14]. Rigid-body docking software, where the relative coordinates of atoms comprising the query protein structures remain unchanged, is commonly used in this kind of approach.
The rigid-body docking method is the most feasible procedure for analyzing the interactome based on protein structures using current computational technology [14]. However, flexibility and dynamics are known to play an important role in PPIs such as in the induced fit. Soft docking software (e.g., RosettaDock [15] and FiberDock [16]), which consider backbone or side-chain flexibility, are advantageous in terms of obtaining structures of protein complexes with high resolution accuracy. One possible strategy for analyzing largenumbers of PPIs, such as the interactome, is to initially perform rigid-body docking based PPI predictions and then use particular protein pairs of interest for flexible docking or structural refinements.
To explore how the current rigid-body docking based method performs on real biological data, we applied a method used in our previous studies [12][13] to reconstruct the well-known bacterial chemotaxis signaling pathway (Fig.  1). The bacterial chemotaxis pathway has been studied for several decades and most of the functional relationships among the proteins involved in this signal process have been identified, especially those involving the core part of the signaling system. However there are still uncertainties concerning how flagellar motor proteins are assembled and operate (reviewed in [17]). Also in the existing databases there are some interactions not listed in conventional network descriptions (Fig. 2).
Here, we conducted PPI network prediction by exhaustive docking using two different docking engines: ZDOCK 3.0.1 and MEGADOCK 2.5. ZDOCK, which is currently one of the most popular rigid-body docking engines, uses a scoring function that includes shape complementarity, electrostatics and a heuristic potential called atomic contact energy [18]. MEGADOCK is a similar system to ZDOCK that searches probable docking structures in a grid-based 3D space using fast Fourier transform (FFT). MEGADOCK employs a much simpler score function in which only shape complementarity and electrostatics is considered and thus makes the calculations 8.8 times faster than ZDOCK. The method is set up to perform massive numbers of calculations that are run on parallel computing systems. There are also other types of useful rigid-body docking tools without a 3D grid-base searching system, such as Hex [19,20], LZerD [21] and PatchDock [22]. Here we used FFT-based docking software because it thoroughly searches whole 3D space whilst at the same time being computationally efficient, which maritima. The dotted circle highlights proteins that make up the flagellar motor. CheC, D and X (circled in bold) are only found in T. maritima, while CheZ (circled by a double line) is only found in E. coli and S. typhimurium. When bacteria sense external signals, such as the presence of chemicals or changes in temperature, chemotactic receptors in the form of a complex with histidine kinase CheA linked to CheW transmit a signal to the cytoplasm by changing the autophosphorylation rate of CheA.In this study, a serine receptor Tsr structure was used. In the unstimulated state, the phosphorylated CheA transfers its phosphate group to CheB and CheY. The phosphorylated form of CheY displays increased affinity for the flagellar motor switch protein FliM. This enhances the clockwise rotation of the flagellar motor, which causes the cell body to change direction in a random fashion (i.e., tumbling motion). When the external signal is of a favorable type, such as the presence of nutrients, it works to reduce the CheA phosphorylation rate. The net effect of this is to lower the CheY phosphorylation level, thus reducing the frequency of changing swimming direction. Consequently, the cell quickly enters a suitable environment for growth. The CheY phosphorylation level is also regulated by phosphatases (CheC, D, X, Z). CheR and CheB play a role in adaptation. Specifically, CheR methylates the activated receptor and CheB demethylates the receptor when phosphorylated. Methylated receptors have a reduced ability to suppress the rate of autophosphorylation of CheA. makes it suitable for large-scale exhaustive docking on parallel computing environments. We also considered using Hex for these experiments because it is one of the fastest docking engines for shape complementarity score function. Nonetheless, we used MEGADOCK as our docking engine with a simple score function that was mainly based on shape complementarity. Specifically, MEGADOCK is designed to be highly efficient in the exhaustive docking calculations and is implemented using parallelizing techniques. MEGADOCK computes each of the docking jobs in parallel using threads by OpenMP. The docking jobs are distributed to thousands of nodes by MPI. In this way, MEGADOCK usually completes the exhaustive docking calculations in just a few minutes, depending on the number of dockings needed and the number of available nodes. MEGADOCK has a further advantage in that the calculation time is almost unaffected by adding an electrostatic term to the score function [23].

Protein Structure Data
We selected the well characterized species E. coli, S. typhimurium and T. maritima as the targets for this study of bacterial chemotaxis. The following procedures were performed to retrieve the data. Firstly, we retrieved the PDB ID list that corresponded to the proteins in the bacterial chemotaxis pathway. Pathway data were obtained from KEGG [24] (KEGG pathway ID: eco02030, stm020230, tma02030 for each species). PDB IDs were obtained through LinkDB [25].
The collected protein structure files were prescreened according to the following criteria, as in the recently published protein-protein docking benchmark version 3.0 [26]: (i) experimental method: X-ray diffraction, resolution better than 3.25 Å; (ii) polypeptides consisting of more than 30 residues. In principle, mutant and synthetic data were excluded with the one exception of CheZ, for which only mutant data were available. Structural data for ligand binding domains of membrane proteins, which are located in the periplasm, were also excluded. Each PDB file was divided into data for one polypeptide chain. Using this procedure we obtained 101 protein structures, which corresponded to 13 protein species (CheA, CheB, CheC, CheD, CheR, CheW, CheX, CheY, CheZ, Tsr, FliM, FliG, FliN). Where possible we used multiple structural data of each protein species. The protein structure data used in this study are shown in Table 1.

Protein-protein Interaction Prediction
The PPI prediction scheme employed in this study consists of two stages. First, we conducted rigid-body docking simulations on all the possible binary combinations of target proteins. Using this process we obtained a group of high scoring docking complexes for each pair of proteins. Next, we analyzed each group of interactions by clustering based on the distance matrix calculated from each central coordinate of the complex. The sufficiently populated clusters were then taken and the docking energy score of the highest ranked complex involved in the selected clusters were used  [34] but it was not included because the structure was unavailable. A total of seven interactions that are not coloured dark grey were found in the STRING database [35] by (i) searching interactions associated with experimental reports or (ii) those annotated in databases (KEGG, BioCyc). The interactions are: CheY-FliG, CheY-CheW, CheB-CheW, Tsr-CheZ, Tsr-CheA, CheR-FliN, CheR-CheZ. These interactions were not considered as "correct" in this study because they have not been characterized.
Tsr Tsr for PPI prediction. The deviation of the selected docking energy score from the score distribution of a high scoring group as Z-score was used to assess possible interactions. This prediction protocol has been described in detail elsewhere [12,13]. Potential complexes that have no other high scoring interactions nearby were rejected. Thus, we consider likely binding pairs have at least one populated area of high scoring structures, one of which may be the true binding site.

Exhaustive Protein-protein Docking
We applied two FFT-based rigid-body protein-protein docking engines, ZDOCK 3.0.1 [27] and MEGADOCK 2.5 [13,23]. Both docking engines exhaustively search 3-D grid space around one protein (receptor) for probable docking interactions by rotating and translating another protein (ligand). The 3D-grid was constructed by a pitch of 1.2 Å. The rotation angle of ligand protein was set as every 15 degrees. The highest docking scores were recorded for each rotation angle throughout all the translation patterns. Using this procedure we obtained 3600 high scoring interactions. In this study we conducted dockings of 101 as receptors 101 as ligands = 10201 structure data combinations.

Clustering of High-scoring Docking Conformations
We used the top 1000 ranked conformations for the postdocking analysis. Before reducing the number of potential hits, results of MEGADOCK were re-ranked by using ZRANK [28] to eliminate the energetically unlikely interactions.
Targeting the 1000 top ranked interactions, we conducted clustering based on the positional differences between the various models. The group average method was used for clustering. Through this method we created at least 20 clusters per docking. The parameters used for the clustering as well as the clustering method itself were carefully chosen in order to obtain the best possible prediction performance by employing a general benchmark dataset as described in our previous study [12]. The distance matrix of interactions was calculated as described previously [12]. Specifically, we take the sum of the two Euclid distances as follows: (i) 3-D coordinates of the ligand proteins when receptor protein positions are superimposed; (ii) 3-D coordinates of the receptor proteins when ligand protein positions are superimposed.

Prediction of Binary Interactions
After clustering analysis we calculated a predicted affinity index E for each protein pair as follows. First we defined sufficiently populated clusters C' from {C | C i to C number of clusters } with the cutoff value of m * , where for the cluster C i , the m i was calculated as the Z-score of the population of C i throughout all the clusters. Second, the interaction with the highest docking score among all the data included in C' was selected as the representative interaction. Finally we calculated the deviation of the docking score of the representative interaction, "E", as the Z-score using scores of 1000 models as the population.
Using the value E we made predictions as to whether a pair of proteins interacts by means of the cutoff value E * . The value of E * was determined independently for each docking engine as the value that yielded the best F-measure when applied to a general benchmark dataset as described previously [12].
In many cases, there were multiple X-ray protein structures available for the docking procedure. For example when we investigate CheA (21 structures) and CheY (51), our docking results provided 21 51=1071 satisfactory results. We then calculate the affinity score (E) for each of these docking results. If the E values are larger than the cutoff (E * ) among any of the 1071 docking results, we consider that CheA and CheY interact. Figure 2 shows the predicted PPIs by ZDOCK and ME-GADOCK. The best F-measure value was 0.52 (ZDOCK, TP=12, TN=57, FP=11, FN=11, recall=0.52, precision=0.52, when E * =7.9 and m * =2.0) and 0.48 (MEGADOCK, TP=14, TN=47, FP=21, FN=9, recall=0.61, precision=0.40, when E * =5.5 and m * =0.0). For both the ZDOCK and MEGA-DOCK predictions, parameter values E * and m * were set as the same values that yielded the best F-measure value when applied to general benchmark data used in a previous study [12].

Predicted PPIs
Previously known PPIs are colored gray in (Fig. 2). Relevant PPIs are defined based on published data [29][30][31][32][33]. The interactions of short form CheA [34] were not consid-ered because its structure was unavailable. In addition, interactions based on genetic observations alone were excluded. FliG, FliM and FliN were considered as binding to the protein species because they make solid flagellar motor machinery. For in vitro studies, large numbers of interactions are listed in public databases such as the STRING database [35]. However, these data sets were not included in this study because the physical interactions for those PPIs are not characterized.
By restricting target proteins using localization information, both ZDOCK and MEGADOCK yielded better Fmeasure values, with both precision and recall values higher than that of the whole dataset. These results show that more accurate PPI predictions are made if protein localization is taken into consideration. Figure 3 shows a comparison of the results obtained by using ZDOCK and MEGADOCK. In total, 17 out of 23 relevant PPIs were detected when at least one of the docking software programs are used. Among 17 true positives, 9 were predicted by both of the software packages. Among the 28 false positives, 4 were common for both software packages and 7 were specific to ZDOCK, while 17 were specific to MEGADOCK. Thus, a lower precision value was obtained using MEGADOCK when compared to ZDOCK.

Performance of PPI Prediction
On the prediction of binary interactions, both ZDOCK and MEGADOCK yielded F-measure values of more than 0.4. When localizations of proteins were considered, ZDOCK performed better than MEGADOCK. It should be noted that MEGADOCK employs shape complementarity and an electrostatic score function whereas ZDOCK also takes into account heuristic score function based on atomic contact energy (IFACE) [18].
As shown in (Fig. 3), eight interactions were detected only by one of the docking software programs. Tsr-CheR, Tsr-CheW and CheA-CheA interactions were detected only by ZDOCK. Tsr-CheD, CheA-CheB, CheD-CheC, CheY-FliN and FliG-FliG interactions were detected only by ME-GADOCK. Two out of three of the interactions detected by ZDOCK, Tsr-CheW and CheA-CheA, are tight binding interactions that constitute the receptor complex. In the case of MEGADOCK, with the exception of FliG-FliG, all the five detected interactions are transient. These results suggest there are differences in the type of interactions detected when different score functions are applied. It is very interesting to see the difference of the predicted PPIs by using different score functions.
Applying other score functions for docking or conducting re-ranking calculations with more sophisticated score functions to the generated decoys would be useful for analyzing the effects of score function on PPI prediction. To investigate this further we require a more thorough dataset such as that used by Kastritis and Bonvin [36] to evaluate any correlation between score function types and known proteinprotein binding affinities.
One such example of applying different PPI prediction procedures is given in (Fig. 4), which shows PPI prediction results for a chemotaxis dataset using the PRISM protocol [37]. PRISM uses a template dataset of known proteinprotein binding interfaces extracted from PDB. The surface of the target monomer protein, for which we want to identify binding partners, is analyzed against all the interface templates by structural alignment. Target protein pairs whose surface structures are aligned to any known interface pairs in the template dataset are then refined and scored by Fiber-Dock [16].
PRISM identified four candidates that potentially interact with most of the proteins in the dataset. Specifically, of the 13 proteins involved in the chemotaxis dataset, CheA and CheZ were predicted to interact with 11 while Tsr and CheY were predicted to interact with 10. Interestingly, the dockingbased procedure applied to CheZ gave fewer predicted binding partners i.e., ZDOCK (three partners) and MEGADOCK (four partners). In this chemotaxis dataset only two interactions are confirmed for CheZ i.e., oligomerization of CheZ and interaction with CheY-p. However, CheZ is also known to interact with the short form of CheA [34] and localize in the cell pole area where receptor complexes are located [38]. CheY, the main target for CheZ, moves between the receptor area and flagellar motor area. The template-based PPI prediction suggests that CheZ may undergo non-specific interactions. Thus, it would be of interest to further analyze the role of CheZ in the receptor complex area.

Protein Localization
Both docking tools yielded better performances when flagellar motor related proteins are excluded from the target, while random prediction with a recall value of 0.5 yielded similar F-measure values (0.34 for the whole dataset and 0.35 for the restricted targets). It should be noted that direct binary interactions among flagellar motor proteins are still unclear and the true interacting pairs might be different from the "correct" interactions used here. Combining protein localization prediction methods such as PSORT [39] and SO-SUI [40], especially for forecasting whether a given protein is membrane associated or soluble, to our PPI prediction would be useful when applying our method to large numbers of target proteins.

False Negative Interactions
When using both ZDOCK and MEGADOCK predictions, 7 interactions were not detected; FliG-FliG, FliM-FliM, FliG-FliM, FliG-FliN, FliM-FliN, Tsr-CheB and CheY-FliM (Fig. 3). When we removed interactions among flagellar motors to consider protein localization, only the interaction between Tsr-CheB and CheY-FliM were not detected by both ZDOCK and MEGADOCK. It is known that for these interactions to occur CheB and CheY both need to be phosphorylated [41,42]. In our dataset there were no protein structures of the phosphorylated forms. However, CheY structures used here include the mimicked activated state by using BeFe 3-, such as PDB ID: 1ZDM [43]. This activated structure showed only modest differences with the native structure in terms of backbone geometry [43]. Based on these results, we cannot assess whether a rigid-body docking method is capable of distinguishing the phosphorylated from the nonphosphorylated state. Nonetheless, our findings are understandable because we did not use a flexible docking tool that considers phosphorylation mediated conformational changes of CheY and CheB.
One possible mean of obtaining increased sensitivity in our PPI prediction model is to construct likely structural variations of target proteins and then use the ensemble set as a docking target.

False Positive Interactions
There are four common false positive PPIs (CheD-CheY, CheA-CheD, CheR-CheY, CheY-CheY) predicted both by ZDOCK and MEGADOCK, three of which include CheY. In total, there are 51 structures of CheY in the dataset and 7 interactions out of 13 target proteins were predicted for CheY by both docking tools. Positive predictions were ob-tained using various structures of single protein species. This, however, does not mean that only specific protein structure pairs generate positive interactions. CheA has 21 structures in the dataset. Both of the docking tools predicted 5 interactions out of 13. The availability of more structural data for a given protein enriches structural variation and serves to increase sensitivity. In such cases we can consider using higher E * value to get better precision.
This result is also understandable from the fact that CheY has multiple binding partners. Bacterial chemotaxis is a twocomponent signal transduction system consisting of a histidine kinase (CheA) and response regulators (CheB, CheY). CheA operates in the form of a complex with receptors and CheW. Phosphorylated CheB works as a modifier of receptor proteins, which accumulate at the cell pole [44]. While CheB operates in the local area around the receptor complex, CheY accepts signals from CheA and transmits them to the flagellar motors, which are evenly distributed around the entire surface of the cell. There are several processes that modify CheY activity during transmission of the signal; CheC, CheX (T. maritima) and CheZ (E. coli, S. typhimurium) have activity that dephosphorylates CheY [33,45,46]. Thus, CheY undergoes transient interactions with several different proteins during the signal transduction process. Indeed, our conclusion that CheY undergoes non-specific binding with many types of proteins is in agreement with our findings given in (Fig. 2). It should also be noted that both docking tools predict CheY undergoes dimerization. Moreover, sequence homology based interlog search using PiSite [2] also suggests that dimerization of CheY is likely.

CONCLUSION
We conducted a reconstruction of the protein-protein interaction network using two distinct physical docking tools. The predicted interactions generated from the two tools were slightly different. However, when the positive predictions from both tools were combined, the vast majority of relevant interactions were represented. Indeed, there were only two exceptions, both of which required phosphorylation to activate the corresponding interaction.
Large-scale PPI prediction using tertiary structures is an effective approach that has a wide range of applications. This approach is especially useful for identifying PPIs associated with novel pathways that control cellular behavior. We are currently attempting to use a PPI network prediction method to find novel interactions in the signaling pathway involving epidermal growth factor receptor (EGFR) with related proteins whose expression levels change in response to specific drugs. In such an application, it is realistic to perform the prediction with a higher threshold (E * ) value where a high level of precision is expected. Using this approach, we can evaluate the likelihood of novel interactions that were previously unconfirmed. In such cases, it would be helpful to employ multiple docking tools to identify as many probable interactions as possible.


s2fieldsofstudy:
Biology
Agricultural And Food Sciences

text:
Comparison of innate immune agonists for induction of tracheal antimicrobial peptide gene expression in tracheal epithelial cells of cattle

Bovine respiratory disease is a complex of bacterial and viral infections of economic and welfare importance to the beef industry. Although tracheal antimicrobial peptide (TAP) has microbicidal activity against bacterial pathogens causing bovine respiratory disease, risk factors for bovine respiratory disease including BVDV and stress (glucocorticoids) have been shown to inhibit the induced expression of this gene. Lipopolysaccharide is known to stimulate TAP gene expression, but the maximum effect is only observed after 16 h of stimulation. The present study investigated other agonists of TAP gene expression in primary cultures of bovine tracheal epithelial cells. PCR analysis of unstimulated tracheal epithelial cells, tracheal tissue and lung tissue each showed mRNA expression for Toll-like receptors (TLRs) 1–10. Quantitative RT-PCR analysis showed that Pam3CSK4 (an agonist of TLR1/2) and interleukin (IL)-17A significantly induced TAP gene expression in tracheal epithelial cells after only 4–8 h of stimulation. Flagellin (a TLR5 agonist), lipopolysaccharide and interferon-α also had stimulatory effects, but little or no response was found with class B CpG ODN 2007 (TLR9 agonist) or lipoteichoic acid (TLR2 agonist). The use of combined agonists had little or no enhancing effect above that of single agonists. Thus, Pam3CSK4, IL-17A and lipopolysaccharide rapidly and significantly induce TAP gene expression, suggesting that these stimulatory pathways may be of value for enhancing innate immunity in feedlot cattle at times of susceptibility to disease.

Introduction
Tracheal antimicrobial peptide is an inducible β-defensin that is highly expressed in the bovine respiratory tract and has demonstrated antimicrobial effects on pathogens that cause bovine respiratory disease, including Mannheimia haemolytica, Histophilus somni and Pasteurella multocida [1][2][3]. TAP expression is induced by several pro-inflammatory cytokines as well as bacterial components such as lipopolysaccharide (LPS) [4][5][6][7]. However, infection with certain viruses and bacteria, or exposure to glucocorticoids, inhibits this response to LPS, implying that this may be one mechanism by which stress and viral infection increase the risk of bacterial pneumonia in cattle [8][9][10][11]. Importantly, the susceptibility of bovine bacterial pathogens to TAP suggests that restoring this suppressed response may be a viable method to prevent disease in stressed or virus-infected cattle [11].
The induction of TAP gene expression by LPS is delayed and does not reach zenith until 16 h of stimulation, perhaps suggesting indirect signal transduction. Thus, the purpose of this study was to determine which TLR receptors are constitutively expressed by bovine tracheal epithelial cells, and to identify agonists that induce high levels of TAP gene expression more rapidly than LPS. Toll-like receptor (TLR) ligands were a focus of the study, including lipoteichoic acid (TLR2/2), Pam3CSK4 (TLR2/1), FSL-1 (TLR2/6), LPS (TLR4), flagellin (TLR5), and CpG ODN 2007 (TLR9). We also examined IL-17A because of its emerging recognition as a critical stimulus for mucosal defence against bacterial infection, and interferon (IFN)-α because of the commercial availability of an immunostimulant (inactivated Parapoxvirus ovis, Zylexis) known to induce its expression [12].

Toll-like receptor and IL-17A receptor subunit A gene expression
Tracheal and lung tissue from three clinically and pathologically normal market-weight calves were obtained at the time of slaughter. From each calf, 30 mg samples of tracheal mucosal tissue and lung tissue were washed in PBS and preserved in RNAlater® at 4°C for 24 h prior to storage at −80°C. Adjacent samples of the same tissues were fixed in formalin for 48 h and histologic sections were prepared routinely. Histologically, the tracheal and lung samples from all calves were within normal limits and had only low numbers of lymphocytes, plasma cells, and few neutrophils.
Fresh tracheas were obtained as above and primary cultures of bovine tracheal epithelial cells (bTEC) were prepared as previously described [10]. Briefly, cells were cultured for 4-6 days until 80-90% confluency in a 1:1 mixture of Dulbecco's modified Eagle's medium and Ham's F-12 medium containing 10% fetal bovine serum, 0.1 mg/mL penicillin-streptomycin, 5 μg/mL amphotericin B, 0.5 mg/mL cell culture grade gentamycin, 25 μg/ mL bovine pituitary extract, 25 ng/mL epidermal growth factor and 1% insulin-transferrin-selenium. To confirm the epithelial nature of the cells, immunohistochemistry was performed on cells fixed in 70% ethanol and suspended in HistoGel™ (Thermo Fisher Scientific Inc.) for pancytokeratin (AE1/AE3, M3515, Dako, ON, Canada) and vimentin (M0725, Dako); negative control sections were processed in the same way with omission of the primary antibody. RNA was extracted from tissues preserved in RNA later® and from cultured bTEC using RNAeasy® Mini Kit capture columns (QIAGEN Inc, Mississauga, ON, Canada) and QIAshredder™ (QIAGEN Inc) according to manufacturer's instructions. Genomic DNA was removed by treatment with DNase I (QIAGEN Inc). The concentration and quality of extracted RNA was measured using a Nanodrop 2000 spectrophotometer (Thermo Scientific; Wilmington, DE, USA). Complementary singlestranded DNA (cDNA) was reverse transcribed from 100 ng total RNA using Super Script® III reverse transcriptase (Thermo Scientific) according to the manufacturer's protocol using oligo dTs as the primer. Once cDNA synthesis was completed, samples were incubated at 37°C for 20 min with ribonuclease H (Invitrogen) to remove any remaining RNA.
PCR amplification of cDNA was completed using primers for TLRs 1-10 [13,14]. For IL-17RA, sequences were obtained from NCBI (XM_002683626.2) and primers were designed using NCBI primer designer software: forward primer 5'-GTGCCCGACTGCAAGGACCC-3' , reverse 5'-CTGGGTGGCTTCCTGTGCGG-3' , 250 bp product size . The PCR profile was as follows: 95°C for  5 min, 95°C for 30 s, 58°C for 30 s for TLR primers   and 65°C for 30 s for IL-17RA primers, 72°C for 30 s,  repeat from step 2 for 35 cycles, 72°C for 5 min. The PCR products were analyzed by gel electrophoresis using 0.01% SYBR® Safe gel dye on a 5% agarose gel (Sigma, St. Louis, MO, USA) in 0.5× Tris/Borate/EDTA buffer (Sigma). Product sizes were compared to a 100 bp ladder, and products of each reaction were confirmed once by sequencing. Pooled cDNA from mesenteric lymph node and ileum obtained from a healthy Holstein calf were used as positive controls for the TLR gene expression study, because these tissues were previously shown to express all 10 bovine TLRs [13]. Negative controls consisted of PCR reactions that included receptor-specific primer pairs but contained no template cDNA.
To quantify relative gene expression, real-time reverse transcription quantitative PCR was conducted using Light-Cycler® 480 technology (Roche Diagnostics, Mannheim, Germany). Constant amounts of cDNA synthesized from 100 ng of total RNA were used in each reaction to measure gene expression of TAP relative to the reference gene glyceraldehyde-3-phosphate dehydrogenase (GAPDH). The GAPDH Ct values remained within 1 threshold cycle across various treatments, confirming constant gene expression under the conditions tested. The primer sequences and product sizes for TAP and GAPH were as follows: TAP forward 5'-TCTTCCTGGTCTGTCTGCT-3' , reverse 5'-GCTGTGTCTTGGCCTTCTTT-3' , 183 bp product size; GAPDH forward 5'-GGCGTGAACCACG AGAAGTAT AA-3' , reverse 5'-CCCTCCACGATGCCAAAGT-3' , 120 bp product size [10]. Amplification efficiencies of 1.91 for TAP and 1.88 for GAP were determined from the slope of a standard curve constructed from a series of 1:5 serial dilutions of cDNA synthesized from LPSstimulated bTEC. A middle dilution was used in each qPCR run as a calibrator to control for plate-to-plate variation. Each qPCR reaction consisted of 10 μL of 2× LightCycler®480 SYBR Green 1 master mix (Roche Diagnostics, Mannheim, Germany), 0.5 μL each of forward and reverse primers (10 μM), 7 μL of PCR grade water and 2 μL of 1:10 diluted template cDNA. The reaction was started by a pre-incubation step at 95°C for 7 min followed by 45 amplification cycles consisting of 95°C for 15 s, 63°C for 15 s, and 72°C for 20 s. Once completed, a melting curve was produced at 95°C for 5 s then 45°C for 1 min. The reaction was concluded by a cooling step at 40°C for 1 min. The samples were run in technical duplicates for both the gene of interest (TAP) and the reference gene (GAPDH) then averaged for advanced relative quantification analysis using LightCycler®480 software (Roche Diagnostics, Mannheim, Germany). Two types of negative controls were used in each trial: a non-stimulation control that differed from the other wells only in the absence of an agonist, and an additional negative control containing master mix but no template cDNA.

Statistical analysis
A general linear mixed model was used to fit the normalized ratios using Proc MIXED (SAS 9.2). Fixed effects were a combination of agonist and dose by time and the combination by time interaction were also included in the model. Calf identification was a random blocking variable (a random effect), and the random interaction of calf identification by doseagonist combination by time was used to handle subsampling (biologic triplicates). To assess the ANOVA assumptions, residual analyses were conducted. The residuals were formally tested for normality using the four tests offered by SAS (Shapiro-Wilk, Cramér-von Mises, Kolmogorov-Smirnov and Anderson-Darling tests). In addition, the residuals were plotted against the predicted values and explanatory variables included in the model. Log transformation and exclusion of 2 outliers produced a fundamentally normal distribution of data. P < 0.05 was considered significant. The data are shown graphically as mean ± standard error of the mean (SEM).

Toll-like receptor and IL-17A receptor subunit A gene expression in bovine respiratory cells and tissue
Using immunohistochemistry, primary cultures of bTEC had positive cytoplasmic labeling for pancytokeratin and about 5% had positive cytoplasmic labeling for vimentin ( Figure 1). The negative control, prepared by omitting the primary antibody, did not show background staining. Unstimulated bTEC, tracheal tissue, lung tissue, and the positive control (pooled RNA from mesenteric lymph node and ileum) each expressed mRNA for all 10 bovine TLRs examined (Figure 2). Findings in cell cultures, tracheal tissue, and lung tissue were similar for all 3 calves examined. The bTEC and the positive control each showed mRNA expression of IL-17RA (data not shown).

Effect of agonists on TAP gene expression
TAP gene expression was measured in bTEC stimulated with TLR agonists and cytokines. Pam3CSK4, an agonist of the TLR2/1 heterodimer, induced TAP gene expression in a dose-dependent manner from 31.6 ng/mL to a maximal effect at the highest dose tested of 3.16 μg/mL. This effect was significantly different from the non-stimulated control and from LPS-stimulated cells ( Figure 3A). At this dose, the magnitude of effect was significantly greater than induction by LPS. The findings were similar in both calves tested, although the magnitude of the response differed. Induction of TAP gene expression was detected as early as 4 h with 1 μg/mL Pam3CSK4 stimulation ( Figure 3B).
FSL-1, a TLR2/6 agonist, also showed significant induction of TAP gene expression. A dose of 31.6 ng/mL induced the highest levels of TAP gene expression ( Figure 3C). Using this concentration in the time course study, FSL-1 induced significantly greater TAP gene expression compared to the negative control and LPS as early as 4 h of stimulation ( Figure 3D). In contrast, lipoteichoic acid, an agonist of the TLR2/2 homodimer, did not significantly induce TAP gene expression in doses ranging from 10-1000 ng/mL ( Figure 3E), although the bTEC used in this experiment did respond as expected to LPS (confirming the viability and responsiveness of these cells). In an earlier time-course study using a lipoteichoic acid dose of 20 μg/mL and bTEC sourced from a different calf, significant induction of TAP gene expression was found at 8 and 16 h, albeit at much lower magnitude than for stimulation with LPS ( Figure 3F).
In bTEC stimulated with flagellin, a TLR5 agonist, a dose of 316 ng/mL (the highest dose tested) significantly induced TAP gene expression compared to the nonstimulated control. This effect was dose-dependent at both 8 and 16 h of stimulation ( Figure 3G). Using this dose, significant induction of TAP gene expression was first detected at 8 h of stimulation ( Figure 3F). However, flagellin did not induce significantly higher gene expression than LPS at any of the doses or times tested.
Stimulation of bTEC with various doses of CpG ODN 2007 (a bovine/porcine TLR9 agonist) and its non-CpG control did not induce significant TAP gene expression at 8 or 16 h ( Figure 4A). A minor non-significant increase was noted for both CpG ODN 2007 and the non-CpG control in the time course study, and a similar nonsignificant trend was found in the dose-response study. However, CpG did not induce gene expression significantly higher than its non-CPG control at any time ( Figure 4B).
Stimulation of bTEC with recombinant bovine IL-17A showed dose-dependent induction of TAP gene expression, to a maximum of 316 ng/mL (the highest dose tested) ( Figure 4C). At this dose, IL-17A significantly induced TAP expression after 8 h of stimulation, with a non-significant increase in TAP gene expression as early as 4 h ( Figure 4D).
Stimulation of bTEC with recombinant bovine IFN-α showed that a dose of 1000 ng/mL (the highest dose tested) significantly induced TAP gene expression in comparison to the non-stimulated control ( Figure 4E). Using this dose of IFN-α, significant induction of TAP gene expression occurred after 4 h of stimulation and remained significantly elevated at later time points. Even though there was a rapid induction caused by IFN-α, it was not higher than LPS ( Figure 4F).
A summative experiment was conducted to compare the effects of Pam3CSK4, IL-17A and LPS on induction of TAP gene expression in bTEC from 4 calves ( Figure 5). As above, all agonists had significantly greater effects than no stimulation. Pam3CSK4 induced significantly higher TAP gene expression than did IL-17A or LPS, at both 8 and 16 h of stimulation for all 4 calves examined. IL-17A also induced significantly higher gene expression than LPS after 8 h, but not at 16 h. Differences among the 4 animals tested were significant (P < 0.0001).
Three experiments (conducted on different days using cells from different calves) examined the combined effects of LPS, Pam3CSK4 and IL-17A on TAP gene expression ( Figure 6). All combinations of agonists tested had a greater effect than IL-17A alone (1-way ANOVA, P < 0.05). In one experiment, the combination of LPS and IL-17A had a greater effect than any single agonist ( Figure 6A), but the difference was small and not found in the second experiment ( Figure 6B). Other combinations of agonists did not have significantly different effects than that of Pam3CSK4 alone or of LPS alone.

Discussion
The present study investigated factors that induce antimicrobial peptide gene expression in bovine tracheal epithelial cells, with the intent of developing methods to stimulate this response to prevent respiratory disease. The morphology and cytokeratin expression of the primary cell cultures confirmed an epithelial phenotype; occasional cells co-expressed cytokeratin and vimentin, which was consistent with prior reports of tracheobronchial epithelial cells [15]. Initial findings identified gene expression of TLRs 1-10 as well as IL-17RA in bTEC as well as in tracheal and lung tissue, consistent with other findings using primary cultured airway epithelial cells [16]. Histologically, tracheal and lung samples from all calves contained lymphocytes, plasma cells, and neutrophils, albeit in low number. This study examined the  presence rather than the level of gene expression, and it is acknowledged that post-translational regulation of these expressed genes would determine functional activity of these cell-surface receptors. However, the subsequent finding that agonists of these receptors induced TAP gene expression implies that functional receptors were expressed. As a result of these findings, agonists for TLRs 1, 2, 5, 6 and 9 and IL-17RA were studied to determine their effect on induction of TAP gene expression.
TLR2 functions either as a homodimer or as heterodimers with TLRs 1 or 6. TLR 2/2 recognizes cell wall components of Gram-positive bacteria such as peptidoglycan and lipoteichoic acid; the TLR 2/2 agonist lipoteichoic acid is a lipid bound to a glycerol phosphate chain. TLR 2/1 recognizes triacylated lipopeptides found in the membrane of Gram-positive and Gram-negative bacteria; Pam3CSK4, a synthetic tripalmitoylated lipopeptide resembling the acylated terminus of bacterial lipoproteins, is an agonist of TLR 2/1. TLR 2/6 recognizes bacterial lipoproteins including FSL-1, an N-terminal fragment of the diacylated LP44 lipoprotein of Mycoplasma salivarium [17,18].
Pam3CSK4, the TLR 2/1 agonist, significantly induced TAP gene expression in comparison to non-stimulated controls and to LPS, and this stimulatory effect was observed in a dose-and time-dependent manner. These results are consistent with studies in humans that have shown up-regulation of human β-defensin-2 (hBD-2)- an inducible β-defensin homologous to TAP-in human airway epithelial cells [19]. Stimulation of human tracheobronchial epithelial cells with Pam3CSK4 induced antimicrobial activity against both Gram-positive and Gram-negative bacteria including Pseudomonas aeruginosa and Listeria monocytogenes [19]. The rapid effect suggests that Pam3CSK4 interacts with TLR 2/1 on the cell surface to trigger intracellular signalling that directly activates transcription of β-defensin genes.
As for Pam3CSK4, stimulation of bTEC with the TLR 2/6 agonist FSL-1 resulted in robust induction of TAP gene expression, in comparison to a non-stimulated control or LPS. In contrast, lipoteichoic acid, the TLR 2/2 agonist, did not induce greater TAP gene expression in the dose-response study than for LPS, and only produced a minimal induction during the time course study. In previous studies, LTA significantly upregulated gene expression of both TAP and lingual antimicrobial peptide in primary bovine mammary epithelial cells. However, when compared to LPS, the effect of LTA in stimulating β-defensin as well as pro-inflammatory cytokine gene expression was significantly lower [5]. Similarly, in a study using human lung epithelial cells, stimulation with LTA showed lower induction of IL-8 gene expression compared to stimulation with LPS [16,20]. Thus, the findings suggest that heterodimerization of TLR2 with TLR1 or TLR6 may augment its function in stimulating this innate immune response.
Flagellin, the major protein of flagella of Gram-negative and Gram-positive bacteria, activates TLR5. In this study, flagellin induced TAP gene expression in bTEC in a doseand time-dependent fashion. However, it never produced effects as high as or earlier than those induced by LPS. With regards to the in vivo relevance, the major Pasteurellaceae bacterial pathogens causing bovine respiratory disease do not possess flagella, suggesting that TLR5 stimulation is not important in pathogenesis of this disease.
Bacterial DNA contains 20 times more unmethylated CG dinucleotides than mammalian DNA, and these CpG motifs are immunostimulatory and recognized by mammalian TLR9 [21,22]. This study used CpG ODN 2007, a bovine/porcine-specific class B ODN, which had no significant effect on TAP gene expression compared to the non-CpG control. This result was consistent with other studies that have shown that CpG ODN had no effect on pro-inflammatory cytokine production in human bronchial epithelial cells [23], although another study found that microbial DNA induced gene expression of hBD-2 in airway epithelial cells by interaction with TLR9 and through activation of NF-κB [24]. CpG may need to be combined with other inflammatory cytokines in order to induce an innate response in epithelial cells; for example, CpG augmented IL-8 production in bronchial epithelial cells co-stimulated with IL-1β and tumour necrosis factor-α, and co-stimulation of mice with class C CpG-ODN2395 and Pam2CSK4 had greater resistance to bacterial challenge than either agonist alone [25][26][27]. Finally, it is noted that other classes of CpG ODNs were not tested in the present study and may have differing effects.
IL-17A is a potent pro-inflammatory and immunostimulatory cytokine produced by T helper 17 lymphocytes, type 3 innate lymphoid cells and bovine γδ-T cells [28,29]. IL-17 is important for defence of the respiratory tract against bacterial pathogens including Staphylococcus aureus, Pseudomonas aeruginosa and Legionella pneumophila infections in mice [30][31][32]. In the present study, IL-17A receptor subunit A gene expression was identified in bTEC, consistent with prior findings in human airway epithelial cells [33]. Exposure of these cells to IL-17A stimulated TAP gene expression in bTEC in a dose-and time-dependent fashion, to a greater degree than did stimulation with LPS. This result is consistent with studies of IL-17A inducing human β-defensin-2 expression in human airway epithelial cells and TAP in bovine mammary epithelial cells [33,34]. These findings suggest that IL-17A induces innate mucosal antimicrobial defences in the bovine respiratory tract.
IFN-α is an antiviral type I interferon. The presence of viral or bacterial nucleic acid within cells stimulates pattern recognition receptors including NOD-like receptors, RIG-I-like receptors, and TLRs 3, 4, 7, 8 and 9, activating interferon regulatory factors to induce IFN-α expression [35]. Secreted IFN-α then interacts with the IFNαβ receptor (IFNAR1/2) to activate the JAK/STAT signalling pathway and induce expression of over 300 immune response genes [36]. IFNAR is expressed in all  cell types including tracheal epithelial cells [37]. In the present study, IFN-α induced TAP gene expression in bTEC only at the highest dose tested and not consistently greater than LPS, although there was significant induction of TAP gene expression as early as 4 h.
Significant differences in agonist-induced TAP gene expression were identified when tracheal epithelial cells from different animals were tested in the same experiment. Similarly, replication of experimental findings using tracheal epithelial cells from a different animal yielded qualitatively similar but quantitatively different results. These findings suggest differences among animals in the magnitude of response to the various agonists. Thus, comparison between agonists is only valid for those experiments in which different agonists were compared in the same experiment.
Caution is warranted because in vitro nature of this study would not address the complex interactions between different cell types that occur in the respiratory tract. However, this is advantageous because it allows a specific focus on epithelial cell responses. Further, it is possible that the surface mucus layer present in vivo might alter the effect of agonists, and this in vitro work avoids this complication. Finally, additional studies are needed to determine whether in vivo stimulation of this immunoinflammatory response would protect calves against pneumonia or increase the risk of disease due to a heightened response to bacteria.
The findings of this study are important for understanding host responses to pathogens and for development of methods to stimulate innate immune responses in cattle. Tracheal epithelium is a dynamic sensor of threats entering into the lung and responds to pathogens by providing a physical barrier, a mucociliary escalator, and by secreting various antimicrobial peptides, immununomodulatory proteins, antioxidants and chemokines that attract other cells of both the innate and adaptive immune responses. Tracheal antimicrobial peptide is secreted from the respiratory epithelium in an inducible manner and has bactericidal activity against the major pathogens that cause pneumonia in cattle [1]. Other functions include antiviral effects, CCR6-dependent chemotaxis of immature dendritic cells and T lymphocytes, and epithelial cell proliferation and wound repair [38]. Prior studies have shown that LPS induces TAP gene expression, but the present findings indicate more rapid effects by stimulation of TLR2 and IL-17RA pathways. Although previous studies have shown that using combinations of agonists produced synergistic induction of β-defensins in airway epithelium [34,39,40], the use of combined agonists had minimal or no synergistic effect in this study. Further, combined stimulation may increase the risks of adverse effects, so identifying specific effective pathways is important and therefore was the focus of this work. TAP has antimicrobial activity against the pathogens that cause pneumonia in cattle, but induction of TAP gene expression in bTEC is blocked by stress (modeled by glucocorticoids) or bovine viral diarrhea virus infection. The present study suggests that TLR2 agonists or IL-17A may be of value to stimulate this aspect of innate immunity in immunosuppressed feedlot cattle.