OTU assignment

Each 16s sequence was assigned to an operational taxonomic unit (OTU). To count the reads number for each OTU (microbe), 16s sequences obtained from MiSeq were aligned to the reference Greengene OTU builds. The Qiime script assigne_taxonomy.py (see more at http://qiime.org/scripts/assign_taxonomy.html) was performed in the data processing. Reference sequences are pre-assigned with OTU described in the id_to_taxonomy file. Any sequence alignment tools, such as uclust, SortMeRNA, blast, RDP, Mothur etc, can be called by the assign_taxonomy script for the sequence alignment between the 16s sequences and reference sequences. For example, the script will assign taxonomy with the uclust consensus taxonomy assigner by default using the following command, assign_taxonomy.py -i repr_set_seqs.fasta -r ref_seq_set.fna -t id_to_taxonomy.txt. OTU redundancy matrix was normalized from the sequence number of each sample. Since these less abundant microbes are unlikely to be a destroying force for host immune system, we selected the top-259 most abundant OTUs for further studying.

Comparison with the controls

To study the remarkable difference of microbiota between cancer cases and the controls, we compared our 16s raw data to those data from 300 healthy human subjects released by Human Microbiome Project (HMP) [22] (http://www.hmpdacc.org). To find a map between OTUs from our data and OTUs from healthy data, a commonly used alignment tool blastn was performed to compare their representative sequences. These pairs with evalue <1e-5 and pident >80% were used for establishing the map. These OTUs matched with a same OTU in HMP were collapsed into one OTU. The Qiime scripts were performed to analyze the 16s raw data [23].

Calculation of correlation

Abundant microbes and DEGs were selected for reconstructing host-microbial networks. DEGs in cervical cancer were collected from published data [9], which were verified in five cohorts of tumor and normal samples. Hence, the DEGs are more reliable than these obtained from only one cohort. The spearman rank correlation method was employed to calculate the correlation between each pair of nodes. Note that, the gene expression data and 16s rRNA were tested on the same sample. Therefore, the spearman correlation in the network makes sense. In contrast to pearson correlation, spearman correlation coefficient can efficiently avoid the environmental noise and experimental errors caused from the non-uniform samples.

Error correction

To improve the confidence of the host-microbial network calculated by spearman correlation, we removed these edges that are less likely to be a true one (false positive errors) and added some new edges that are very likely to correlate with each other (false negative errors). The false positive edges include two scenes: 1) these negatively correlated edges that connected two interactors with a same type of regulation (i.e. both of them are up regulated or down regulated); 2) these positive correlated edges that connected two interactors with different types of regulation (i.e. one is up regulated, the other down regulated); 3) self-loops; 4) multiple loops. All these false positive edges are removed in our network. These false negative edges are these pairs of nodes between OTUs and DEGs which satisfying two conditions: 1) the OTU was collapsed from a set of sub-nodes; 2) all these sub-nodes strongly correlated with the DEG. All these false negative edges were added in the host-microbial network. False positive and false negative edges were detected and corrected according to the coherence of regulation and correlation relationships. A workflow of the reconstruction of host-microbial network was illustrated in Fig. 1.

Bi-partite betweenness centrality

To search for risk factors from host-microbial network, bi-partite betweenness centrality (BpBC) [24], adapted from betweenness centrality, was used to quantify the risk of a given node, written as g(v). The definition can be formatted as \(g(v)=\sum _{s,t}\delta _{st}(v)/\delta _{st}\). Here, s and t are two nodes from two separate sub-networks. And δ _st represents the number of shortest paths from s to t, δ _st (v) the number of shortest paths going through node v from s to t. Given a node v, g(v) reflects the probability of how likely a shortest path could go through v from one sub-network to another.

Composition of the microbiota

To study the microbial community in cervical cancer, we examined the 16s raw data of cancer cases and assigned taxonomy to each sequence. The definition of operational taxonomic unit (OTU) was used to classify groups of closely related microbiome based on sequence similarity. Reference data sets and id-to-OTU maps for 16s rRNA sequence was downloaded from the Greengenes reference OTU builds [25]. All these sequences were grouped into different categories based on their family-level OTU labels. As shown in Fig. 2, prevotellaceade followed by tissierellaceae appears to be the most abundant microbes, accounting for 13.7% of the microbiota community. There are four other groups accounting for more than 5% of the microbiota, which are fusobacteriaceae, porphyromonadaceae, planococcaceae and bacteroidaceae. Totally, twenty-six family-level OTU groups make up more than 87% of the whole community. To examine the diversity of cervical microbiota, the PCoA analysis was performed to analyze the microbial community in cervical carcinoma, skin, mouth and vagina. As shown in Fig. 3, microbiota in cervical carcinoma (red dots) is less diverse than microbiota in any other body sites. Hence, we indeed found remarkable changes of microbial composition in the cancer cases.

Reconstruction of host-microbial network

A host-microbial network was reconstructed from the 16s raw data and gene expression data. Nodes in the network refer to microbes or DEGs, edges the regulation relationships between each pair of microbes. Two nodes were connected if and only if they are strongly correlated (i.e. |γ|>0.4 and p-value < 0.05). As show in Fig. 4, a network with 997 nodes was connected by 4262 edges. Nodes in the network consist of 259 microbes and 738 DEGs. We grouped all the DEGs into four categories, named as cell cycle, antiviral response, epithelial cell differential and the other DEGs, according to their function in the development of cervical cancer. The three functional DEGs groups (excluding the other DEGS) are three major densely connected sub-networks in the host-microbial networks. They are functionally enriched by GO terms cell cycle, response to virus, epithelial cell differentiation respectively. They don’t have any overlap between each pair of groups. In the whole network, 403 edges are negatively correlated, 3859 positively correlated. Negative correlation indicates inhibition between two biological subjects. In a negative correlation, one variable increases as the other decreases. Positive correlation indicates activation or co-existence between two subjects of interest. In a positive correlation, one variable increases as the other increase, or one variable decreases while the other decreases. This network integrates all the regulation relationships between host genes and microbiota.

Risk factors in cervical cancer

The risk factors in cancer may activate oncogene expression and cause a series of functional disorder in metabolic and immune systems. In the development of cancer, the most remarkable differences between tumor and normal samples are: 1) the up-regulation of viral responses; 2) the speed-up in the progression of cell cycle; 3) the inhibition of epithelial cell differentiation. To study how microbiota regulates the viral response, cell cycle and epithelial cell differentiation, we searched for key risk factors using BpBC. These key factors are thought to be cancer-causing agents that can drive the progress of cervical carcinogenesis. Nodes that organizing communication between two cancer-related groups are more likely to be key factors. Since BpBC is such a measure to evaluate the importance of a node in the network topology, we choose these nodes in the top list of BpBC as candidates of key factors. These key factors with high BpBC value may play crucial roles in the communication between two different sub-networks.

The results show that Anaerococcus (labeled as OTU_97.18428) and proteasome subunit beta 10 (PSMB10) are significantly higher than the others (see in Fig. 5 left) between the sub-networks of microbe and antiviral response genes. PSMB10 was an up-regulated gene in cervical cancer. Between the sub-networks of microbe and cell cycle, KCNIP1 and Hydrogenophilaceae (labeled as OTU_97.2777) are the most important regulators (see in Fig. 5 middle). Eubacterium (labeled as OTU_97.10051) and KRT13 are the most important regulators between the sub-networks of microbes and epithelial cell differentiation (see in Fig. 5 right). It proves that the interplay between microbiota and differentially expressed genes might be the driving force that regulates the progress of cell cycle, epithelial cell differentiation and viral response.

Query and visualization

In order to fast and easily query and visualize the host-microbial networks, we developed a web-based tool KF-finder. Multiple web programming languages were used in the development, which includes PHP, mysql and javascripts. Each node and its neighborhood in the network can be searched by a query term in the panel of Search. And the induced sub-network will be visualized in the panel of View. For example, one can input a gene symbol CYP2A7 as a query term in the Search panel. A list of nodes associated with CYP2A7 will show out in a user-friendly panel, as well as a graphic view of the induced sub-network (see in Fig. 6). Except for visualization and query, KF-finder can also sort microbes and DEGs in a decreasing order by the value of BpBC in microbe-antivirus, microbe-cell cycle or microbe-epithelial cell differentiation. Download and advanced search have been enabled on the web server. All our test datasets and results of users’ personal jobs can be downloaded. Advanced search allows us search for genes and microbes based on string patterns or value constriction. KF-finder enables us to query and visualize the knowledge of host-microbial network in a fast-and-easy way. It can be accessed at http://www.nwpu-bioinformatics.com/KF-finder/.

A case study of PSMB10 in cervical cancer

Most vertebrates express immunoproteasomes (IPs) that possess three IFN- γ-inducible homologues: PSMB8, PSMB9 and PSMB10. Many studies show that expression of IP genes including PSMB10 is up-regulated in most cancer types [26]. IP genes can be expressed by non-immune cell, and that differential cleavage of transcription factors by IPs has pleiotropic effects on cell function. Indeed, IPs modulate the abundance of transcription factors that regulate signaling pathways with prominent roles in cell differentiation, inflammation and neoplastic transformation (e.g., NF-kB, IFNs, STATs and Wnt) [27]. Therefore, PSMB10 is indeed a risk factor involved in the antiviral response of cervical caner.

A case study of KRT13 in cervical cancer

KRT13’s full name is keratin 13 in human, also known as K13 and CK13, located in a region of chromosome 17q21.2. It is a down-regulated gene in cervical carcinoma, and a risk factor that involves in the progress of uncontrolled epithelial cell differentiat,ion. Previous work suggests that the loss of K13 or low K13 mRNA expression is associated with invasive oral squamous cell carcinoma (OSCC) [28, 29]. Epigenetic alteration of K13 is one major reason resulting the inhibition of K13 in OSCC. Besides, K13 was also reported that it played a directive role in prostate cancer bone, brain and soft tissue metastases [30]. More than 1000 single nucleotide polymorphisms of K13 were found in the dbSNP database. Totally, 51 variations mentioned K13 in ClinVar, seven out of which are pathogenic. All these evidences suggest KRT13 is very likely to be a key risk factor involved in cervical cancer.