A total of 128 blood samples of known thalassemia genotypes were collected from Fujian Medical University Union Hospital (Fujian, China) and the First Affiliated Hospital of Sun Yat-sen University (Guangzhou, China). The Samples were collected from October 2016 to January 2017. There were 79 female and 49 male patients, with an age range of 4 months to 86 years (mean age, 26). Peripheral blood samples (~5 ml) were collected into tubes that contained ethylenediaminetetraacetic acid. For each sample, genomic DNA was extracted from 100 µl whole blood using the DNeasy Blood and Tissue kit (Qiagen, Inc., Germantown, MD, USA) according to the manufacturer's protocol. Briefly, the blood samples were hydrated with 200 µl Buffer AL and 20 µl proteinase K followed by incubation for 10 min at 56°C, the contents were transferred to a DNeasy Mini Spin Column placed in a 2-ml collection tube following the addition of 200 µl ethanol (96–100%). The samples were centrifuged at 6,000 × g for 1 min at room temperature, at 6,000 × g for 1 min at room temperature following the addition of Buffer AW1, and again at 20,000 × g for 3 min at room temperature following the addition of Buffer AW2. Finally, the samples were eluted with 200 µl Buffer AE, quantified on a Qubit^® fluorometer (Life Technologies; Thermo Fisher Scientific, Inc.), and stored at −20°C prior to use.

The primer sequences were designed using reference sequences of the HBA2 and HBB gene loci [accession nos. NC_000016.9 (222846.223709) and NC_000011.9 (5246696.5248301)] from the NCBI database (http://www.ncbi.nlm.nih.gov/nuccore/NC_000016.9 and http://www.ncbi.nlm.nih.gov/nuccore/NC_000011.9) with Ion AmpliSeq™ Designer (https://www.ampliseq.com/). The Ion AmpliSeq™ Thalassemia Panel, which consists of two primer pools made up of 82 pairs of primers (72 pairs of primers for HBA2 and 10 pairs of primers for HBB), was designed by Life Technologies; Thermo Fisher Scientific, Inc.

Each sample was used to construct the library using the Ion AmpliSeq™ Library Kit 2.0 (Life Technologies; Thermo Fisher Scientific, Inc.). In brief, 10 ng genomic DNA, 4 µl 5X Ion AmpliSeq™ HiFi mix, 10 µl 2X Ion AmpliSeq™ primer pool, and 4 µl nuclease-free water were mixed to amplify the target regions. Subsequently, 2 µl FuPa reagent was added to each amplified sample to partially digest the primer sequences, and each library was ligated into a unique barcode and a universal adapter provided in the Ion Xpress™ barcode adapters (Life Technologies; Thermo Fisher Scientific, Inc.). Each library was purified using AMPure XP beads (Beckman Coulter, Inc., Brea, CA, USA). The purified libraries were quantified on a Qubit^® 3.0 fluorometer. The size distributions of the libraries were verified using the Agilent High Sensitivity DNA kit on a 2100 Bioanalyzer (Agilent Technologies, Inc., Palo Alto, CA, USA).

Each library was diluted to 100 pM according to its quantified concentration as determined on the Qubit^® 3.0 fluorometer. Subsequently, one test making up 14 or 15 libraries of 100 pM was emulsion PCR-amplified with Ion PGM™ Hi-Q™ ion sphere particles (ISPs) using the Ion OneTouch™ 2 Instrument (Life Technologies; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. The template-positive ISPs were enriched using the Ion OneTouch™ ES instrument (Life Technologies; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol.

The enriched templates were loaded onto one Ion 318™ chip V2 and sequenced on the Ion Torrent Personal Genome Machine (PGM; Life Technologies; Thermo Fisher Scientific, Inc.), a semiconductor sequencing platform.

Sequencing data was mapped to the human reference sequence hg19 (Genome Reference Consortium GRCh37). The variants were called (Torrent Suite v.4.4.3; Life Technologies; Thermo Fisher Scientific, Inc.) using variant calling software with optimized parameters for the thalassemia panel. The variants were annotated using Annovar (16) and the system's software. The detected variants were subjected to a rigorous manual curation process, which included querying variant databases, including the SNP database (www.ncbi.nlm.nih.gov/snp/), Exome Aggregation Consortium (exac.broadinstitute.org/), 1000 Genomes database (www.internationalgenome.org/1000-genomes-browsers) and Clinvar database (www.ncbi.nlm.nih.gov/clinvar/) and a literature review.

The CNV was calculated by counting the reads in each amplicon with MAQ>10. The reads were first uniquely mapped to the hg19 sequence from the RAW.bam file. The CIGAR index in each read was then trimmed. The read counts with >50% uniquely mapped in one amplicon were calculated. For certain amplicons, the reads were calculated with MAQ>10, as their low mapping quality would lead to multiple hits, which included amplicons in the HBA1 and HBA2 genes. The same protocol was performed again using MAQ>0 as the control group.

A novel algorithm was developed to identify the CNV types of α-thalassemia. In these cases, the target amplicons were related to different types of α-thalassemia regions, as described above. The algorithm consisted predominantly of four tests: A ratio, which revealed the α-thalassemia-^SEA deletion type; B ratio, which revealed the α-thalassemia-α^4.2 deletion type; C ratio, which revealed the α-thalassemia-α^3.7 deletion type; and D ratio, which represented the compound heterozygous or homozygous deletion.

Initially, several basic parameters were defined, including the sequence read numbers of the ith reference amplicon (ref-reads-i) and the ith target amplicon (AMPL reads-i). The ref-reads was defined as the average number of reads of the five reference amplicons. The control reads ratio (AMPL-i) was defined as AMPL-i=(AMPL reads-i/ref-reads), and 28 control reads ratio values were obtained from 34 normal samples as a baseline. Other parameters were defined as follows: Median (median value of a cluster of numbers): Reads ref=Σ (ref-reads-i)/4 (i=3, 4, 8, 9 and 10); test reads ratio=(AMPL reads-i/ref-reads); A ratio=median (test reads ratio/control reads ratio) (i=4, 7, 8, 9, 10, 12, 13, 15, 44, 45, 48, 49, 50, 51, 55, 57, 58, 60, 65, 66 and 67); B ratio=median (test reads ratio/control reads ratio) (i=20, 21 and 22); C ratio=median (test reads ratio/control reads ratio) (i=32, 33 and 35); D ratio=median (test reads ratio/control reads ratio) (i=27). GraphPad Prism (v. 5.0; GraphPad Software, Inc., La Jolla, CA, USA) software was used for all statistical analyses. Data are expressed as the mean ± standard error of the mean and were analyzed using an unpaired Student's t-test (two tailed).

All samples were amplified using the α-Thalassemia Genetic Diagnostic kit (gap-PCR method; DaAN Gene Co., Ltd., Sun Yat-sen University, Guangzhou, China). The target products were detected by agarose gel electrophoresis.

The HBA2 and HBB target genes were amplified using specific primers, and the target products were sequenced by Sanger sequencing. New primers were designed using Primer Premier 5.0 software. The primer sequences for HBA2 were: Forward 5′-CCCCACATCCCCTCACCTACATTC-3′ and reverse 5′-CGGGCAGGAGGAACGGCTAC-3′; the primer sequences for HBB were: Forward 5′-CAGAAGAGCCAAGGACAGGTACGGCT-3′ and reverse 5′-AAGGGCCTAGCTTGGACTCAGAATAATCC-3′.

The present study aimed to establish a method of simultaneously detecting CNVs and SNVs of thalassemia that can be applied to other diseases, including autism spectrum disorder (ASD), spinal muscular atrophy (SMA), and Duchenne muscular dystrophy (DMD). The samples for the study were selected with the aim of including as many different types of thalassemia as possible. Samples with sequencing reads ranging between 100 and 500 M were selected for analysis. The average number of total raw bases was 45,236,536 (range: 22,632,244–100,007,680). The average read length was 155 bp. The mean percentage of sequencing reads mapped to the reference hg19 genome was 98%. Following filtering of the low-quality reads, polyclonal reads and primer dimer reads, sequenced bases with ³Q20 values ranged between 21,042,640 and 94,594,388 (Fig. 1).

In the present study, SNVs were identified in the target region. Approximately 11 SNVs were identified according to the reference human genome hg19, including eight SNVs with clear definition of pathogenic alleles recorded in the Clinvar database, of which five SNVs were located in HBB exons, one in HBA2 exons, and two in introns or upstream of the HBB gene. The most frequently mutated gene locus was NM_000517.4:c.427T>C (HBA2), which resulted in a termination codon mutation in the HBA2 gene to glutamic acid, making it difficult to continue synthesis of the polypeptide chain until the next stop password. Two samples carried a nonsynonymous variant causing p.Val114Glu, and three samples carried frameshift variants. The results were consistent with the results of Sanger sequencing (Figs. 2A and B and 3).

There were three other nonsynonymous variants identified in exons which were not recorded in the Clinvar database. The SIFT score (17–19), Polyphen2_HDIV_score, Polyphen2_HVAR_score (17,18), and PROVEAN score (20), which were used to predict whether an amino acid substitution or indel affected the biological function of a protein, were calculated to evaluate the possible adverse effects (i.e., deleterious or possibly damaging nature) of the nonsynonymous variants on protein function. However, the SNV carriers exhibited symptoms of thalassemia, supporting the prediction results (Tables I and II).

An additional 66 SNVs were identified that did not result in an amino acid change and were located in an intron, an intergenic region, or upstream or downstream of genes. The minor allele frequency (MAF) value of 20 SNVs in the 1000 Genomes database was <0.01, indicating that these SNVs occur less frequently in the normal population. However, their potential adverse effects require further evaluation. The MAF value of 32 additional SNVs in the 1000 Genomes database was >0.01, indicating a probable polymorphism, and 1–94 of the 128 samples in the present study carried these SNVs. These SNVs may form their own polymorphism in Chinese individuals, providing evidence for gene haplotype and crowd site distribution. Of these SNVs, 14 had no information in the 1000 Genomes database or other databases (Table III).

The human HBA gene cluster, located on chromosome 16, spans ~30 kb and includes seven loci: 5′-zeta-pseudo-zeta-mu-pseudo-alpha-1-alpha-2-alpha-1-theta-3′. The α-2 (HBA2) and α-1 (HBA1) coding sequences are identical. The similarity of these gene sequences is almost 97%. They differ only marginally in their 5′untranslated region and introns and differ significantly in their 3′untranslated region. The target CNVs of HBA2 depend on an accurate alignment algorithm to avoid ambiguity between HBA2 and HBA1.

The present study introduced the concept of mapping quality, a measure of the confidence that a read actually comes from the position it is aligned to by the mapping algorithm. Align MAQ can build assemblies by mapping shotgun short reads to a reference genome using quality scores to derive genotype calls of the consensus sequence of a diploid genome (21). In the present study, six -^SEA/-α^3.7 samples and three -^SEA/-α^4.2 samples were analyzed using Align MAQ=10. The aligned reads of AMPL-25, AMPL-29, AMPL-30, AMPL-37, and AMPL-38 in the -^SEA/-α^3.7 samples using Align MAQ>10 were close to 0 compared with those using Align MAQ=0 (Fig. 4A and B). Similar results were found in the -^SEA/-α^4.2 samples.

Applying reference amplicons is key to constructing an algorithm to detect CNVs. For an algorithm to be accurate, the reference gene region should be a stable diploid with minimal variation in the amplicon sequencing depth of different samples. According to thalassemia disease-associated genes, regions of the HBB gene (ref-03-chr11: 5246753–5246986, ref-04-chr11: 5246976–5247184, ref-08-chr: 5248047–5248296, ref-09-chr11: 5248286–5248485, and ref-10-chr11: 5248475–5248641, hg19), which encodes β-globin, were selected as reference amplicons. The HBB gene was used as the endogenous reference gene as β-thalassemia is predominantly caused by SNVs in the HBB gene, rather than a CNV. For thalassemia of the HBB CNV types, other genes require selection as the reference gene. The sequencing depth at each base pair position in these five regions was counted in all 128 samples divided by seven different groups. Normalized ref-reads were generated and are shown for each sample. The values varied between 1,136 and 13,282, with no significant differences in the amplicon sequencing depth among the samples in the seven groups, with the exception of -^SEA/αα (Fig. 5A). The abnormal value in the -^SEA/αα group may have been caused by deletions in the HBA2 gene region, although without influence on the final results.

The following step was to investigate the consistency of the samples. A cluster of reference reads ratios of 28 amplicons were built as a baseline across 33 normal samples (Fig. 5B). The reads ratio was defined as the ratio of the target region reads to the reference region reads of each sample. Examination of the coefficient of variation (CV) of the reference samples revealed that 24 of the 28 amplicons had CVs with values <41.1% (Fig. 5C).

To identify an indicator for CNV detection, a novel algorithm was developed based on the ratio of the median reads ratio of the target sample to that of the reference. The median ratio value, but not the mean ratio value, was used to evaluate the CNV type as the middle value is less vulnerable to a deviation as a result of a sequencing error. The A ratio, B ratio, C ratio and D ratio revealed the copy numbers of the region related to the Southeast Asia deletion, the -α^4.2 deletion, the -α^3.7 deletion, and the compound deletion type of α-thalassemia, respectively. The A ratio ranged between 0.741 and 1.298 in the normal group, between 0.263 and 0.899 in the -^SEA/αα group, between 0.246 and 0.898 in the -^SEA/-α^3.7 group, and between 0.232 and 0.707 in the -^SEA/-α^4.2 group. The discrepancy in the A ratio was significant (P<0.0001) between the normal (αα/αα) samples and heterozygous Southeast Asia deletion type (−^SEA/αα) samples according to Student's t-test (Fig. 6A). Consistent with the heterozygous Southeast Asia deletion type (−^SEA/αα, -^SEA/-α^4.2, and -^SEA/-α^3.7), the fluctuations in the B ratio and C ratio associated with the -α^4.2 deletion and -α^3.7 deletion were similar to that of the A ratio. The discrepancies were also significant (Fig. 6B and C). The D ratio was defined as the ratio of the AMPL-27 reads ratio of the target sample to the reference median reads ratio. The AMPL-27 ranged between chr16: 223333 and chr16: 223548 in the HBA2 gene (HBA1 and HBA2 genes encode ~97% of the total Hb). A homozygous deletion in this region indicates a severe type of thalassemia. AMPL-27 is a common deletion region in these three types. Therefore, the D ratios in the -α^3.7/-α^3.7, -^SEA/-α^3.7, and -^SEA/-α^4.2 groups were close to zero (Fig. 6D).

The Southeast Asia, -α^4.2, and -α^3.7 deletions were identified using the following criteria: A ratio <0.8, B ratio <0.4, and C ratio <0.8, respectively. Subsequently, gap-PCR was used to evaluate the sensitivity and specificity of the approach. A total of 61 heterozygous Southeast Asia deletion (−^SEA/αα, -^SEA/-α^4.2, and -^SEA/-α^3.7) samples were detected with 96.72% (59/61) sensitivity and 93.94% (31/33) specificity, 12 heterozygous -α^4.2 deletion (−α^4.2/αα and -^SEA/-α^4.2) samples were detected with 83.33% (10/12) sensitivity and 100% (33/33) specificity, and 38 -α^3.7 deletion (−α^3.7/αα, -α^3.7/-α^3.7, and -^SEA/-α^3.7) samples were detected with 97.37% sensitivity (37/38) and 93.94% (31/33) specificity. Compound homozygous thalassemia was identified using the following criterion: D ratio <0.002. In total, 20 homozygous deletions of AMPL-27 were detected with 95% (19/20) sensitivity and 100% specificity (33/33).

As NGS technology is able to simultaneously detect target gene CNVs and SNVs, their correlation was investigated in the present study. When a gene exhibits a loss of heterozygosity, only a haploid gene exists, not a diploid. Once this gene acquires SNVs, 100% frequency can be detected; this abnormal sample is defined as a compound heterozygous CNV and SNV. In the present study, certain samples had compound heterozygous CNVs and SNVs (e.g., -^SEA/αα and CD122).