A father-son pair were both diagnosed with AML at The General Hospital of Western Theater Command (Chengdu, China). The son was diagnosed in April 2014 while the father was diagnosed in May 2016. The father was aged 66 and the son aged 44 years. According to National Comprehensive Cancer Network Guideline for AML (8), both patients were diagnosed with AML and vein blood samples were collected, and DNAs were extracted for whole-genome sequencing (WGS) or WES using the phenol-chloroform method. The Hospital Ethics Committee approved the study and the two patients signed informed consent before donating samples. For transcriptome sequencing, the RNA sample was extracted from the blood of the son.

The DNA sample from the father was used for genome-wide sequencing. DNA was extracted from peripheral white blood cells using a DNA extraction kit (Qiagen GmbH) according to the manufacturer's instructions. A total of 2 µl purified DNA was used to measure the quality and quantity of the DNA (ng/µl) (A260/A280) using an UV visible spectrophotometer (BioTek Instruments, Inc.). Finally, 20 ng/µl DNA samples were stored at 20°C, until further use. DNA degradation and the molecular weight was analyzed using 0.2% agarose gel electrophoresis (loading 5 µl extracted DNA). DNA libraries were prepared following the manufacturer's protocol (Illumina DNA Prep with Enrichment; cat. no. 20025523; Illumina, Inc.). Briefly, for each sample, DNA was firstly sheared into fragments of ~350 bp. The fragments were then end-repaired, A-tailed, ligated to paired-end adaptors and PCR amplified (according to the manufacturer's instructions) for library construction. The final library concentration used was 6.2 pM. The resulting DNA libraries were subjected to 150 bp pair-end sequencing on the Illumina HiSequencing PE150 platform (Illumina, Inc.). The reference genome was downloaded from the University of California, Santa Cruz (UCSC) database (GRCh38/hg38; http://genome.ucsc.edu). Low-quality reads were discarded using the PRINSEQ software (v0.20.4) and the resulting clean data were aligned to the human reference sequence (hg19) with the Burrow-Wheeler Aligner (BWA) software (0.7.12-r1044) (9). Duplicate reads were removed using the Picard software (http://sourceforge.net/projects/picard/). The aligned reads were sorted with the Genome Analysis Toolkit (GATK) software (v2.8.1; http://github.com/RRafiee/GenomeAnalysisToolkit). SIFT (sift.jcvi.org) and PolyPhen2 (genetics.bwh.harvaed.edu) software were used to predict the effects of amino acid mutations on protein structure and function. All variants were annotated with ANNOVAR (10).

The DNA sample from the son was used for WES. The protocol of DNA extraction was identical to that of the father's sample. DNA was randomly digested into 150–200 bp fragments using a Bioruptor Pico ultrasound. Fragmented DNA was end-repaired, A-tailed and then connected. Sample labeling and enrichment of DNA were conducted using PCR amplification according to the manufacturer's instructions. The DNA library with the specific index was subjected to liquid phase hybridization with the biotin-labeled RNA probe, and the target gene exon was obtained using the streptavidin-labeled magnetic beads. Then, target exon genes were enriched by PCR amplification. Sequencing was performed using an Ilumina platform (AmpliSeq™ Focus Panel for Illumina; cat. no. 20019164; 10 ng) and the HiSeq 3000/4000 PE (150 paired-end) Cluster kit (both from Illumina, Inc.), according to the manufacturer's instrustions. Clean reads were mapped to the reference genome, GRCh37 using BWA. After removing duplications, single nucleotide polymorphisms, and insertions and deletions (Indel) were assigned and annotated using the (GATK) based on dbSNP build 150.

The total RNA from blood samples of the son was extracted using TRIzol^® reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. RNA purity was determined using a NanoDrop 2000 spectrophotometer (NanoDrop Technologies; Thermo Fisher Scientific, Inc.), and RNA integrity was verified using 1.5% agarose gels. Magnetic beads conjugated with oligo (dT) were used to isolate mRNA, which was fragmented into short lengths of ~200 bp using an RNA Fragmentation kit (Illumina, Inc.). The short RNA fragments were used as templates for first-strand DNA synthesis with random primers according to the manufacturer's instruction. Next, second-strand DNA was generated, purified using magnetic beads, end repaired, and a single adenine (A) nucleotide was added to the 3′ ends. DNA degradation and the molecular weight were analyzed using 0.2% agarose gel electrophoresis (loading 5 µl extracted DNA). RNA-Sequencing (RNA-Seq) libraries were prepared using the VAHTS Total RNA-seq (H/M/R) Library Prep kit for Illumina (Vazyme Biotech Co., Ltd.) following the manufacturer's instructions. The raw reads were processed to remove the adapter and primer sequences using SeqPrep software (v1.0; http://github.com/jstjohn/SeqPrep). Low-quality bases at the 3′ ends were deleted using Sickle software (v1.33; http://github.com/najoshi/sickle). Fragments <30 bp in length were excluded from further analyses and reads containing N >10% were removed. The obtained high-quality sequencing sequence after quality control was aligned with the designated Ensembl genome using TopHat (http://tophat.cbcb.umd.edu/). Saturation, duplicate reads and coverage analysis were conducted software from Majorbio, including SeQC-2.3.2 (http://code.google.com/p/rseqc/) and RSeQC-2.3.2 (http://code.google.com/p/rseqc/). Predictionof novel transcripts in annotated genomes was performed using Cufflink (http://cufflinks.cbcb.umd.edu/) (11). To identify the known long non-coding (lnc)RNAs, the predicted novel transcript was aligned with a known lncRNA using collection of databases, such as NONCODE (www.noncode.org), Ensembl (ensembl.org), NCBI (http://www.ncbi.nlm.nih.gov/), and UCSC (genome.ucsc.edu), LncRNAdb (http://www.lncrnadb.org/), GENCODE (https://www.gencodegenes.org/) and LncRNA Disease. The gene expression levels and differentially expressed genes (DEGs) were estimated using cufflinks (http://cole-trapnell-lab.github.io/cufflinks/) (12–14). The default parameters for DEGs were false discovery rate (FDR) ≤0.05 and log2 fold-change (FC) ≥1 or ≤-1.

Gene Ontology (GO) enrichment analysis was performed using the software Goatools (https://github.com/tanghaibao/GOatools). The P-values were corrected using four multiple test methods (Bonferroni's, Holm, Sidak and FDR) to control the calculated false positive rate. When the corrected P-value (or P-FDR) ≤0.05, it was considered significant (15–17). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using KOBAS (http://kobas.cbi.pku.edu.cn/home.do). A corrected P≤0.05 was considered to indicate a significant difference (18). Samtools (http://samtools.sourceforge.net/) and VarScan version 2.2.7 (http://varscan.sourceforge.net/) software were used to find candidate SNPs.

The WGS results for the father's blood sample were as follows. A total of 638,386,172 raw reads were generated as shown in Table I, achieving WGS mean depth of 25. The coverage of the whole genome was 96.9%. All properly pair-end mapped sequences were used for subsequent variants detection. A total of 3,897,164 SNPs and 780,834 indels were identified (Fig. 1). The distribution of each SNP locus was analyzed (Fig. 2A) (Table II). The results showed that C>T and G>A were the most common variations. The distribution of indel lengths is presented in Fig. 2B. All SNP and indel variants were annotated using the ANNOVAR program in conjunction with the reference gene annotation information of the UCSC Genome Browser. A total of 24,150 SNPs located in the exon region were probed and 727 indel variants in the exon region were found. The distribution of SNPs and indel variants in gene functional regions shown in Table II.

Regarding amino acid translation, there were 11,316 non-synonymous, 12 stop-loss, 12,033 synonymous, 92 stop-gain for SNPs, 63 frameshift insertions, 73 frameshift deletions, 242 non-frameshift insertions, 248 non-frameshift deletions, 4 stop-gain, 2 stop-loss for indel variants (Table III).

Non-synonymous SNPs can cause changes in the encoded amino acids. Together, 2,436 mutations were destructive to protein structure and function according to SIFT, and 1,154 mutations were probably harmful according to PolyPhen2. Comparing with 121 AML-related genes identified in previous AML studies, 14 genes were found in the father's exon region, including the following genes: HPSE2 (19), NT5C2 (20), AKR1C3 (21), MBL2 (22), ARID5B (23), TET1 (24), WT1 (25), CD44 (26) and FLT3 (27). However, only 10 of these had non-synonymous mutations changing encoded amino acids.

The WES results from the son achieved 96,396.31 Mb raw data, generating a mean depth of 28.61. The coverage rate was 99.91%. All properly pair-end mapped sequences were used for subsequent variants detection, as shown in Table II. Overall, 96,639 SNPs were identified, including 10,504 non-synonymous SNPs. For the types of base mutations in SNPs, C>T and G>A were the most commonly observed in this sample. The indel results showed 15,033 mutations, including 102 frameshift insertions, 117 frameshift deletions, 156 non-frameshift insertions, 190 non-frameshift deletions, three stop-gain and one stop-loss (Table III). Seven mutant genes were found in son's exon region compared with 121 AML-related genes, including FLT3, GATA2 (28,29), MCM7 (30), PTEN (31), and RUNX1T1 (32) (Table IV). Among them, six genes showed non-synonymous SNPs (Table IV).

Then the mutant genes in the exon region were compared between the father and son results, and some common ones were observed. Notably, two missense mutations in the FLT3 gene were identified. Internal tandem duplications in the FLT3 gene (FLT3-ITD) were the earliest described molecular alterations in AML (33). In the present results, the alterations of FLT3 were as follows. In WGS data from the father: Exon 6, c.C680T: p.T227M (NM_004119). In WES data from the son: Exon 1, c.A20G: p.D7G (NM_004119) and exon 12, c.A1423C: p.T475P (NM_004119). Another common gene was PTEN. As Tables III and IV show, in WGS data from the father, the following mutations were identified: c.G10A: p.G4R (NM_001304717) In exon 1, c.G194C: p.C65S(NM_001304717) in exon 2 and c.C511G: p.L171V (NM_001304717) in exon 2. In WES data from the son the c.G194C: p.C65S (NM_001304717) mutation was identified in exon 2.

Transcriptomic sequencing generated 69,537,564 raw reads containing 10,430,634,600 nucleotides for the son. Following clean-up and quality filtering, 68,019,440 and 75,105,808 clean reads were obtained, with a Q20 percentage (proportion of nucleotides with quality value >20) >98.51%, which indicated the RNA sequencing results were of high quality and suitable for use in further analysis.

As Fig. 3 shows, there were 54 differentially expressed mRNAs from transcriptomic analysis, of which 31 were significantly upregulated in AML, and 23 were significantly downregulated. However, no common differentially expressed mRNAs between WGS and WES were identified. As aforementioned, FLT3 and PTEN might play an important role in AML development. Compared with patients with AML, expression of FLT3 was decreased (P=0.035; FDR, 0.919; log2FC, 3.17). A lower expression of PTEN were found in the patients with AML without statistical significance (P=0.4135; FDR, 0.993; log2FC, −0.71).

GO annotations revealed that the aforementioned genes were enriched in ‘single-organism signaling’, ‘response to stimulus’, ‘biological regulation’ and ‘metabolic process’. The top 30 terms were all among the biological process class and are presented in Fig. 4. KEGG pathway analysis revealed 56 related pathways, while according to corrected P-value, only three pathways were significantly enriched: ‘Influenza A’, ‘hematopoietic cell lineage’ and ‘cell adhesion molecules’ (Fig. 5). Alternative splicing is an important mechanism for regulating gene expression and producing protein diversity, which can cause large differences in genes and proteins (12). There were 4,370 3′ end alternative splicing events, 2,575 3′ untranslated region (UTR) alternative splicing events, 4,927 5′ end alternative splicing events, 1,783 5′ UTR alternative splicing events, 2,468 exon spanning events, 336 intron retentions, and 2,761 other alternative splicing.