A punctured HCC tumor and paired adjacent tissue was obtained from a patient (57 years, male) at the Youan Hospital, Capital Medical University of China (Beijing, China) and complied with the principles of The Declaration of Helsinki. The patient was infected with HBV and received no radiation and chemotherapy prior to radiofrequency ablation.

The DNA was extracted using an E.Z.N.A.^® Tissue DNA Kit (Omega Bio-Tek, Inc., Norcross, GA, USA) and the extracted DNA was captured using Agilent Human All Exon 50 M kit (Agilent Technologies, Inc., Santa Clara, CA, USA) following the protocols recommended by the manufacturer. Sequencing machines generated a large volume of data at a rapid speed by sequencing paired-end DNA fragments in parallel using Illumina His-seq2,000 (Illumina, Inc., San Diego, CA, USA) (17,18). Following a series of library construction and actual sequencing, a large quantity of raw data was produced.

Raw reads generated by a sequenator are usually affected by adverse factors, including adaptor contamination, poor base sequence quality and guanine-cytosine (GC) bias (19). Once the raw data was obtained, the quality of raw reads was assessed and the adaptor was clipped using fastq-mcf (version 1.04.636; www.github.com/ExpressionAnalysis/ea-utils/blob/wiki/FastqMcf.md). The sequencing data was then processed using the FastQC tool (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to analyze the distribution of base GC content and sequence quality scores.

Following the quality control analyses, the processed reads were aligned to an established reference genome (version hg19), which was provided by the University of California Santa Cruz (Santa Cruz, CA, USA) (20). Millions of short reads were aligned efficiently to the reference genome using Burrows-Wheeler Aligner (BWA) software with default parameters, which were based on the Burrows-Wheeler transform (21). The aligned reads were then stored in BAM file (.bam) using samtools software (22), which was able to sort and index the BAM file to save space and help subsequent process. For the assembled genome data, the picard tool (http://picard.sourceforge.net/index.shtml) was combined with bamtools to filter out the mismatching and inappropriate reads. In addition, picard removed the read duplicates derived from library PCR. The data distribution and reads coverage were then evaluated using the CalculateHsMetrics package. Recalibration and realignment were performed using GATK (version 2.8; Broad Institute, Cambridge, MA, USA; www.broadinstitute.org/gatk/). Finally, the resulting data were used for further variation identification.

A key step in the analysis of cancer exome sequencing data is the identification of variants. The depth of sequence coverage determines the choice of somatic mutation algorithms used for identification of variants mutation. The different identification abilities in different allele frequencies of GATK (version 2.8.1), MuTect (version 1.1.4; Broad Institute; http://www.broadinstitute.org/cancer/cga/mutect), and VarScan (version 2.3.6; http://varscan.sourceforge.net/), and the joint analysis strategy by combining the three softwares (the three-caller pipeline approach), were taken into consideration when identifying somatic mutations.

Oncotator (http://portals.broadinstitute.org/oncotator/) was used to annotate the screened variations (23). All of the candidate mutations were validated visually using the Integrated Genomics Viewer (IGV) (24) and were confirmed using Sanger sequencing in paired samples. The tools, Polyphen-2 (www.genetics.bwh.harvard.edu/pph2/index.shtml) and scale-invariant feature transform (SIFT; www.sift.jcvi.org/), were integrated to predict whether mutations affected protein function based on the structure and function of the protein, and the conservation of amino acid residues in different species sequences.

The gene sets screened were used for functional annotation analysis by the Database for Annotation, Visualization and Integrate Discovery software (25), which consists of the Kyoto Encyclopedia of Genes and Genomes and Gene Ontology database. The significance of gene groups enrichment was defined by a modified Fisher's exact test and P<0.05 was considered to indicate a statistically significant difference.

WES was analyzed in one HCC tumor and paired adjacent tissues with the three-caller approach. The present study acquired 96.30X and 79.18X coverages for the tumor and paired adjacent tissues, respectively, in all of the targeted exonic regions, with 93.4% of the base targeted at 20-fold and ≥99.1% bases by a depth of at least two times. To identify the somatic mutations, a flow chart was created with the following steps: i) Quality evaluation of the raw reads; ii) reads map to a reference genome; iii) somatic mutation identification with the three-caller approach; iv) variant annotation; v) data visualization; and vi) pathway analysis (Fig. 1).

Variant filtering was performed by GATK with the following filter parameters: Low coverage (DP <5), low quality (QUAL >30.0 and QUAL <5.0), very low quality (QUAL <30), hard to validate [MQ0 ≥4 and MQ0/(1.0*DP)>0.1)] and quality-by-depth (QD <1.5). The exome data from the samples were calculated by running these parameters and reserved in a VCF file. GATK was primarily used for identifying somatic mutations in the sequencing data, including SNVs and indels.

In order to identify the low allelic-fraction mutations, MuTect was used to generate more performance in low coverage (12). To illustrate how high the sensitivity was based on allele fraction and sequencing depth, a strategy was established based on the published data to analyze the data (14). As shown in Fig. 2, the sensitivity of mutation was detected by MuTect approaching >90% at allele frequency 10% with >80X sequencing depth and 80% at allele frequency 5% with >80X.

The calling of SNVs by MuTect software was executed through Java (version 1.6.0_45; www.oracle.com/technetwork/java/javase/downloads/java-archive-downloads-javase6-419409.html). The default parameters of MuTect were kept to identify mutations. Input database texts, including reference sequence hg19, dbsnp v.135 and cosmic v54, were used for the MuTect algorithm. Somatic point mutations were only identified by MuTect; GATK (version 1.5) was used to analyze indels. SNVs located in exome regions were screened with ≥20 coverage in the tumor, which was coupled with ≥4 alternate alleles and ≥4 allelic fraction of the altered base. The paired normal sample also had 10X coverage at least in a certain base. As many low coverage or low allelic fraction SNVs were characterized by MuTect, SNVs with variants from low purity samples not blindly rejected.

VarScan outperformed the other tools at higher allelic fraction. A threshold of 6X for tumor and 8X for normal was set, with ≥20% variation frequency. Subsequently, the present study preferentially analyzed 20X coverage in the tumor, including alternated variation accounting for 10X coverage, to eliminate false positives.

The present study proposed 75 candidate somatic variants through the three-algorithm strategy (Fig. 3), including 50 nonsynonymous mutations, 2 nonsense mutations, 20 synonymous mutations and 3 indels. The nonsynonymous to synonymous somatic SNV ratio was 2.5.

The predictive impact of amino acid substitution on functional evidence was analyzed using PolyPhen-2/SIFT (Table I). The P94Q mutation was predicted to affect the protein function of cell division cycle 7 protein, which may be associated with neoplastic transformation of some tumors and affect protein serine/threonine kinase activity. All of the putative somatic mutations were validated manually using IGV. The T>C transversion at position_9056725 in mucin 16 (MUC16) was identified (Fig. 4), which was then validated by Sanger sequencing.

The 75 genes with tumor-specific mutations demonstrated significant functional enrichment of cell adhesion and regulation of Ras GTPase activity (P<0.05; Table II). Notably, the genes encoding cell adhesion demonstrated the most prevalent enrichment (P=0.0089), indicating that the enriched mutations of cell adhesion genes may serve pivotal roles in HCC development.