A total of 21 PCa patients who underwent radical prostatectomy between 2011 and 2014 at the Department of Urology, Yokohama City University Graduate School of Medicine, were included in the present study. Several parts of each resected prostate that had been indicated to contain cancer tissues by preoperative examinations were embedded in OCT compound (Sakura Finetek Japan) and immediately stored at −80°C. The patient clinical information is summarized in Table I. All the patients were Japanese, with a mean age of 67 years (range, 51–76 years), and serum prostate-specific antigen (PSA) values ranging from 4.4 to 31.0 ng/ml (mean, 10.4±7.36 ng/ml). All tumors were diagnosed as non-metastatic adenocarcinomas and assigned a Gleason score of 6–9 at the Department of Pathology. When multiple Gleason scores had been assigned to one patient, the highest score was used, as shown in Table I.

In order to profile the genetic alterations in Japanese patients with PCa in an efficient as well as highly sensitive manner, two original panels were prepared for the targeted sequencing of genes that were reported to be altered in previous whole-genome or whole-exome sequencing studies in Western countries. The KCC71 panel was for DNA samples designed to detect single-nucleotide variations (SNVs) and small insertions and deletions (indels) in 71 PCa-related genes and driver genes reported by whole-exome and whole-genome sequencing (4–11,14,15) (Table SI). The PCaFusion panel was for RNA samples designed to detect transcripts derived from 38 previously reported fusion transcripts, together with 8 control transcripts (Table SI). In addition, the PCaFusion panel was designed to detect fusion transcripts with different exonic junctions from the same fusion partners (1–5,14–18). The multiplex-PCR primer sets for the KCC71 and PCaFusion panels are provided in Tables SIIA and SIIB.

To obtain DNA and RNA samples, a thin-sliced section was prepared from each stored frozen specimen, embedded in OCT compound and stained with hematoxylin and eosin (HE). Based on information on the area of tumor tissues in the HE-stained section, PCa tissues were obtained directly from the remaining OCT-embedded specimen, from which DNA and RNA were extracted using ZR-Duet DNA/RNA miniprep (Zymo Research), following the manufacturer's protocol. DNA and RNA were quantified with Qubit 2 (Thermo Fisher Scientific, Inc.). To assess the DNA and RNA quality, ratios of optical densities, A260/A280 and A260/A230, were further evaluated by NanoPhotometer (Implen). A total of 10 ng of genomic DNA or total RNA was used to create panel libraries for each specimen. Library amplification was performed using Ion Torrent AmpliSeq™ technology, along with sequencing with the Ion PGM next-generation sequencer (Thermo Fisher Scientific, Inc.).

Torrent Suite v4.0.2 and Ion Reporter version 4.4 (Thermo Fisher Scientific, Inc.) softwares were used to process and analyze the sequenced data from the Ion PGM. Quality control reports were obtained from the Torrent Suite. To identify somatic variants, the SNVs and indels with a coverage rate of ≥20 and with coding amino acid sequence substitutions when compared with the UCSC hg19 reference genome sequence were first selected. Next, single-nucleotide polymorphisms (SNPs) were excluded by using the sequences as queries against the data in the databases COSMIC (http://cancer.sanger.ac.uk/cosmic), dbSNPs (NCBI, NIH; http://www.ncbi.nlm.nih.gov/projects/SNP/), the 1000 Genomes Project (http://www.1000genomes.org), and other publicly accessible databases. For SNVs for which it remained unclear whether they were SNPs or somatic variants after database analysis, Sanger sequencing on DNA from the non-neoplastic counterpart of each specimen was performed to obtain definitive results. Finally, sequence alterations with an allele frequency of >5% were defined as somatic variants in the present study.

Torrent Suite v4.0.2 and Ion Reporter version 4.4 were used to process and analyze the sequenced data from the PCaFusion panel. Quality Check reports were obtained from the Torrent Suite server. The unclear fusion transcripts identified by the PCaFusion panel were further verified by reverse-transcription (RT)-PCR followed by Sanger sequencing of the products.

A summary of the sequencing statistics is presented as a representative case of KCC71 panel analysis in Fig. S1A. The means of the obtained reads and coverage were 3,902,663 and 1,850, respectively. The data on average alignment ratios revealed that 97.6% of the total reads were aligned properly to the hg19 human genome reference sequence. A summary of the sequencing statistics as a representative case of the PCaFusion panel analysis is shown in Fig. S1B. The mean number of obtained reads was 539,860. The data on the average alignment ratios revealed that 88.1% of the total reads were aligned properly to the hg19 human genome reference sequence.

As indicated in Materials and methods, confirmed somatic nucleotide sequence alterations of non-synonymous SNVs and indels, with or without frameshifts, were considered as somatic variants in the present study. Somatic variants were detected in 17 of 21 patients by the present panel analyses. No variants were detected in 4 patients (cases 18–21). A total of 33 somatic variants were identified in 24 of 71 PCa-related genes in the KCC71 panel. The results are summarized in Fig. 1 (detailed information is provided in Table SIII).

A total of 7 probable pathogenic variants in 5 genes were identified in the present KCC71 panel analyses. Evident driver gene variants in the literature were found in BRAF (p.K601E) (22) and SPOP (p.F102V and p.F133L) (5,6). Although not well characterized as driver genes in the literature, somatic variants with a high pathogenic score predicted by FATMM (23–25) were also identified in CDH1 (p.E880K) and RB1 (p.R621S). Two variants found in BRCA2, namely p.I1859fs and p.R2318ter, which may result in premature termination and truncation of the BRCA2 polypeptide, were considered as pathogenic, although the identical alterations did not appear in COSMIC.

The remaining 26 variants were considered as variants of uncertain/unknown significance (VUSs), including AR (p.K610E), CDH1 (p.G62V), FOXA1 (p.R265-K267 del) and TP53 (p.V31I). These variants were found in the COSMIC v82 database with labels of ‘n/a’ or ‘neutral’ based on FATHMM score. The CDH1 (p.G62V) variant is not registered in COSMIC; however, the identical mutation was reported as a germline mutation detected in families with hereditary diffuse gastric cancer (26). This non-synonymous mutation was in a region encoding a pro-domain and is generally considered to be non-pathogenic (23–25). The remaining 21 somatic mutations did not appear in COSMIC.

The existence of gene fusion transcripts was analyzed to identify the presence of fusion genes with the original PCaFusion panel. All 8 non-fusion transcripts evaluated as positive controls were detected in all specimens. In the present study, fusion transcripts were designated as follows: The 5′ gene symbol (number of the exon located at the fusion site)-the 3′ partner gene symbol (exon number). For SLC45A3-ELK4 fusion transcripts, SLC45A3(1)-ELK4(2) and SLC45A3(1)-ELK4(4) were detected in all cases. By contrast, SLC45A3(2)-ELK4(2) was identified in only 1 case (case 5). USP9Y-TTTY15 fusion transcripts were detected in all examined cases. Among the TMPRSS2-ERG fusion gene transcripts, TMPRSS2(1)-ERG(4) was identified in only one case (case 11). No other fusion transcripts were identified in the PCaFusion panel analysis. The results are summarized in Fig. 1.

A comparative analysis of the results of the KCC71 panel with TCGA and other big data registered in cBioPortal (http://www.cbioportal.org) was performed. Briefly, the frequency of somatic variants in the aforementioned public databases were examined, including the databases of TCGA, Broad Institute, Freed Hutchinson Cancer Research Center, and Memorial Sloan Kettering Cancer Center (hereafter referred to as ‘cBioPortal databases’) for the 71 genes in the KCC71 panel, and this was compared with the frequency obtained in the present study. The total frequency of somatic variants in the 71 genes, calculated as the total variant number identified per examined case, was higher compared with that in the cBioPortal databases (present study, 33 different variants in 21 cases; summary of cBioPortal databases, 849 variants in 1,656 cases) (Fig. S2, and Tables SIV and SV). This difference was particularly notable for the frequencies of variants in ATBF1, BRCA2 and LRP1B, which were all ≥10% compared with those in the cBioPortal databases. By contrast, the frequency of variants of TP53 was low (1 in 21 cases, 4.8%; Fig. 1). To compare those mutation frequencies in terms of pathways, the ratio between the number of patients with and without mutations in the genes in a particular pathway was calculated. Then, the ratios from our database and TCGA databases were compared using Fisher's exact test. We observed that those ratios in genes belonging to the androgen receptor (AR; -log10 Fisher's P-value=2.67) and DNA Repair (-log10 Fisher's P-value=4.40) signaling pathways were particularly high compared with those in TCGA database (Tables II, SIV and SV).

Our network analysis consisted of four steps as explained below (Fig. 2A-D). In the first step, genes with mutations from the KCC71 panel analysis were selected as ‘initial seed genes’ (Fig. 2A). In the second step, the 7 gene networks of PCa in TCNG were selected as graphical representations of the regulatory relationships between genes. The Gene Expression Omnibus (GEO) ID and information from each selected gene network for PCa are shown in Fig. 2B and Table SVI. In the third step, the ‘initial seed genes’ were mapped on the 7 PCa gene networks, and a subnetwork around the ‘initial seed genes’ was extracted from each of the gene networks; each subnetwork consisted of the ‘initial seed genes’ and downstream genes within the two passes around the initial seed genes (Fig. 2C). Genes with one path from ‘initial seed genes’ are referred to as ‘child genes’ and genes with a path from the child genes are referred to as ‘grandchild genes’. The obtained subnetworks show the regulation around the seed genes. Next, we attempted to identify ‘extended common seed genes’ that are shared among ≥6 subnetworks. In the last step, a putative ‘core network’ of PCa was estimated by integrating subnetworks around ‘the extended common seed genes,’ referred to as ‘extended subnetworks’ (Fig. 2D). The above network operations were conducted by using the functions of igraph, a package of R version 3.5.0. The subnetworks (relationships) of gene regulation were demonstrated by graphical visualization using Cytoscape software version 3.5.1 (19), as shown in Fig. S3, and the core network was presented using igraph.

Two publicly available tools, Database for Annotation, Visualization, and Integrated Discovery (DAVID; Laboratory of Human Retrovirology and Immunoinformatics, http://david.ncifcrf.gov/home.jsp) (20) and Reduce + Visualize Gene Ontology (REVIGO; Redjer Boskivic Institute; http://revigo.irb.hr) (21), were then used to investigate the biological functions associated with the gene groups involved in the core network. DAVID v6.8, which mainly provides typical batch annotation and Gene Ontology (GO) term enrichment analysis, was used to highlight the most relevant GO terms associated with a given gene list, in order to elucidate the biological meaning behind a large list of genes. Enrichment analysis was performed using DAVID for the core network genes listed in Table SVII, and a functional annotation chart including a list of GO IDs and P-values for the enrichment tests was obtained. The functional annotation chart report of DAVID shows categories, enriched terms associated with a gene list of interest, related term search, genes involved in the term, and percentages or modified Fisher's exact P-values. REVIGO was used to summarize the results obtained from DAVID (21). REVIGO provides a functional interpretation of genes defined by GO with statistical methods. A list of GO IDs and P-values was entered from the functional annotation chart report of DAVID. The REVIGO GO tree map shows a two-level hierarchy of GO terms.

The workflow for the reconstruction of the core network of PCa was schematically summarized with TCNG (Fig. 2A-D). The 24 genes with somatic variations from the KCC71 panel analysis (Fig. 1) were selected as ‘initial seed genes.’ Although only 5 well-characterized pathogenic driver gene mutations were identified in the present analysis and the remaining 19 genes were considered as VUSs, all ‘initial seed genes’ were reported to be involved in PCa in the literature and databases (Materials and methods, Sample preparation and target sequencing with the original panels). Subnetworks were extracted from 7 public PCa gene networks in TCNG (Table SVI), which included initial seed genes (parent nodes), and parent-child and child-grandchild genes in the network (Fig. S3). ‘Extended common seed genes’, TNK2, SOX9, CDH1, FOXA1 and TP53, that were commonly included in the subnetworks, were extracted. To identify the core network of PCa s examined, extended subnetworks around the ‘extended common seed genes’ were further extracted from each of the original 7 PCa networks, integrated, and the core network was finally reconstructed.

The core network around the ‘extended common seed genes’ was further analyzed. The 3 extended common seed genes, SOX9, CDH1 and FOXA1, were connected via edges with each other through the genes EMX2, NKX3-1 and TFAP2A, and formed a closed loop (Fig. 3, red arrows). EMX, NKX3-1 and TFAP2A were the only genes in the network located between the extended common genes. The top 15 genes with high connectivity (hub genes), are listed in Table SVIII. All 5 extended common genes are listed in Table SIV, but none of the initial seed genes appears in it. Only AR and SPOP as the initial seed genes appear in the final network, with few edges.

The enrichment analysis using DAVID and REVIGO found 50 GO terms with adjusted P-values (Table SIX). REVIGO generates tree maps of the GO terms, as shown in Fig. 4 and Fig. S4. The GO terms are joined into ‘superclusters’ of loosely related terms and depicted with different colors. The most significant GO term found was ‘positive regulation of transcription from RNA polymerase II promoter’, followed by ‘epithelial cell differentiation’, ‘response to water deprivation’, ‘tissue homeostasis’, and ‘amino acid transport’.