For the 454 and Illumina methods, NGS raw data is comprised of luminescent or fluorescent images, respectively, recorded after each iterative sequencing step. Each platform has its own algorithm for image processing, image compilation, signal-to-noise ratio assessment, and conversion into linear sequence. Nucleotide bases are assigned individual “quality” scores which share a conceptual analogy to Phred scores introduced in 1998 for estimating accuracy and quality of the Sanger sequence data [11, 12]. A Phred score provides quality values, q, using a mathematical scale converting the estimated probability of an incorrect call, e, to a log scale:

Miscall probabilities of 0.1 (10 %), 0.01 (1 %), and 0.001 (0.1 %) yield increasing confidence Phred scores of 10, 20, and 30, respectively. Factors that contribute to lower quality bases in NGS include nonoptimal signal-to-background noise levels, cross talk between signals from adjacent beads or clusters, and a phenomenon termed dephasing due to unequal extension of individual reads in a clonal cluster or on a bead surface [13]. The signal output from Ion Torrent is a change in pH, and obtaining optimal signal-to-background pH ratios requires careful control of reagent pH. As with the 454 and Illumina methods, Ion Torrent signals are converted to linear sequence with quality scores assigned to individual bases. While each platform has a unique signal-to-sequence processing algorithm, instrument-specific sequence files are convertible to a common text-based file format termed FASTQ that contains sequence reads and individual base quality scores.

FASTQ-formatted sequence reads are used for two major computational purposes: assembly and alignment. Assemblies are performed when no reference genome exists for the sequenced DNA, an example being a previously uncharacterized pathogen. Assembly algorithms search for overlapping sections of sequence reads and join them together to generate longer length “contigs” that serve as a scaffold for genome assembly and subsequent alignments. To achieve more complete assemblies, longer sequence read lengths and sequencing from both ends of library fragments (termed paired-end read sequencing) are desirable. For most clinical applications, however, reference genomes exist (e.g., the human genome reference or a previously sequenced pathogen), and the primary computation is alignment of reads to the reference sequence. The unique features of NGS, notably shorter read lengths and massive read numbers, have spurred the development of many new alignment, variant calling, and annotation algorithms. While the reader is referred to a recent review on the topic of NGS bioinformatics [14] for greater details, major process steps are described next (shown in Fig. 59.3).

Alignment is the process of determining the best match between the sequencing reads and the reference sequence. Due to the large number of sequence reads that must be processed, NGS alignment algorithms typically employ one of two approaches to decrease computational time, and these include a data compression method referred to as a Burrows-Wheeler Transform (BWT) and a method based on the concept of a hash table. The BWT approach is integral to the popular open-source Burrows-Wheeler Aligner (BWA) algorithm that has become a standard for sequence alignment [15–17]. For hash table-based alignment, a first step is conversion of the reference sequence or the sequence reads into a population of shorter length sequences (termed an index), in which each sequence is given a read identifier that can be computationally tracked. By using indexes of shorter length (referred to as “seeds”), alignments proceed more rapidly [14]. An example of a popular aligner that utilizes a hash table function to create an index from the reference genome is NovoAlign (Novacraft, Selangor, Malaysia) [18]. For both BWT and hash table-based aligners, key criteria for initial mapping and alignment are established, including the number of nucleotide mismatches permitted across a given read or seed length and whether gaps in alignment are allowed to accommodate insertions and deletions (indels). With initial criteria applied, a first mapped and aligned read data set is generated. This data set is known to contain inaccuracies, and therefore additional algorithms are used to yield a more refined and accurate final set of alignments as described below.

Most alignment algorithms can provide results in a sequence alignment/map (SAM) file format, which contains information about read position and orientation in relationship to the reference sequence and confidence in the alignment. A binary version of SAM is the BAM format, and the SAM/BAM formats are the most commonly used file formats for storing and secondary manipulation of read alignments. After initial alignment, SAM/BAM files are used as inputs into secondary algorithms that refine and increase alignment accuracy prior to identifying differences between the sequencing reads and the reference sequence, a process referred to as variant calling. Popular open-source software for refining alignments and calling variants are the Genome Analysis Toolkit (GATK) [19] and SAMtools [20]. Algorithms within the GATK program are the most widely used for refining Illumina reads. Three major refinement steps are (1) marking and removal of reads with the same start and end points, referred to as PCR duplicates; (2) local realignment to improve accuracy in identifying indels; and (3) recalibration of base quality scores. PCR duplicates arise from sequencing identical fragments generated by PCR during library preparation. PCR nucleotide errors can be introduced and propagated through unequal amplification of fragment library templates, and these errors can manifest as false-positive variants. During removal of PCR duplicates, only a single read with overall highest base qualities is retained for analysis. The purpose of recalibrating base quality scores is to compensate for the fact that the Phred-like quality scores generated by the Illumina platform have been demonstrated to deviate from the true error rate. Figure 59.4 shows an example of how local realignment improved the accuracy of detecting a 3 bp deletion in the FOXP3 gene. The mapped and aligned reads are visualized in the open-source Integrative Genomics Viewer (IGV) [21, 22]. Subsequent to initial and refined alignment steps, variants are called and tabulated in the Variant Call Format (VCF) File that contains variant chromosomal position, reference base, and the alternate base(s) (e.g., single-nucleotide variant (SNV), indel).

By virtue of NGS libraries being comprised of randomly overlapping fragments, multiple reads align to the reference in a staggered fashion. This multiplicity provides quantitative information reflected in the number of reads aligned to the reference sequence at a given location and is designated as read coverage depth. Another feature of randomly fragmented libraries is that fragments are sequenced in both directions to yield forward and reverse strand sequences, and ideally their percentages should be approximately equal. If a variant is present, then the percentage of reads that contain the variant can be expressed as an “allelic read percentage,” and the number of forward and reverse sequence strands containing the variant can be enumerated. Figure 59.5 shows a representative example of a heterozygous SNV and illustrates the concepts of coverage and allelic read percentage. Several points and caveats are important when considering coverage and variant calling. First, greater variant calling accuracy is achieved when there is a consensus among aligned reads that agrees with predictions for the variant. For example, a sample with a heterozygous SNV would ideally display a 50/50 ratio of reads containing variant and reference nucleotides, and forward and reverse strand reads would be equivalently represented in both variant- and reference-containing reads. By extension, a homozygous variant would be expected to be present in 100 % of aligned sequences with equal representation of forward and reverse strands. In practice, a wider range of allelic read percentages for true heterozygous and homozygous variants is observed. Further, forward and reverse strand read percentages divergent from expected are commonly seen (referred to as “read strand bias”). These read strand biases arise from technical and bioinformatic sources including differential PCR amplification of library fragments, sequencing errors in difficult-to-sequence regions, and misalignment of related sequences (e.g., pseudogenes or highly homologous genes).

Second, the minimum number of reads required to call a variant has not been standardized and varies depending on application needs and sequencing platform accuracy (i.e., lower sequencing error rates require fewer reads to accurately establish a variant call). In the 1000 Genomes Project [23], coverage as low as four- to six-fold has been used in a large-scale screening mode, whereas in other research studies 20- to 30-fold coverage has been empirically found to offer a balance between sensitivity, specificity, and sequencing costs for detection of constitutional variants [24–26]. By extension, clinical laboratories are typically requiring a minimum of 20- to 30-fold coverage for the detection of constitutional variants and are often designing diagnostic assays so that much higher coverages are routinely obtained in an effort to increase variant detection and improve variant call confidence.

Third, coverage depth across sequenced regions is variable, and factors that contribute to this are differential ligation of adapters to fragments during library preparation and differential amplification during clonal expansion [27]. NGS assays can be designed to achieve sufficient sequencing to meet clinically necessary minimum coverage for areas prone to lower coverage.

Fourth, whereas the above noted 50/50 and 100 % ratios are pertinent to identifying constitutional heterozygous and homozygous variants, respectively, they do not apply when identifying somatic variants in cancer samples comprised of a mixture of tumor and normal cell populations (where allelic read percentages for true somatic variants of interest can be quite low), as well as for tumor heterogeneity in which not all tumor cells have acquired the same variants. To identify low allelic read percentage somatic variants, increased read coverage in the several hundred- up to thousand-fold range are being applied empirically for clinical testing [28, 29]. In sum, multiple variables ultimately need to be considered in variant calling, and more sophisticated algorithms are beginning to emerge for the detection of variants in both constitutional and somatic testing applications. For example, the open-source GATK program factors in coverage depth, read bias, recalibrated base quality scores, and variant prior and posterior probabilities based on comparison to control samples. In combination, these factors are used by GATK to generate statistical measures that predict variant accuracy and generate a “refined” VCF file.

A final step is to input the VCF file information into an “annotation” program that ascribes additional features to variants. Open-source software that contain annotation functions are ANNOVAR [30, 31], GATK, and snpEff [32], among others. Annotation outputs can include many features, but those that are common are chromosomal location of base change from reference, whether the variant is in a gene and its respective location (e.g., exon, intron, splice site), the consequence of the change to a codon (e.g., synonymous vs nonsynonymous, missense vs frameshift), and zygosity (e.g., homozygous or heterozygous). Often incorporated into annotation software programs are algorithms that predict the functional impact of variants on proteins such as Sorting Intolerant from Tolerant (SIFT) [33, 34], PolyPhen-2 [35], and MutationTaster [36, 37].

Investigations into the genetic basis of a growing number of inherited disorders have revealed that a clinical phenotype can be due to multiple causative genes with a broad mutational spectrum. Examples include X-linked mental retardation, mitochondrial disorders (secondary to mutations in both the mitochondrial genome and nuclear genes), congenital hearing loss, cardiomyopathies, and primary immune deficiencies [8, 10, 38–41]. A comprehensive diagnostic approach can require the analysis of a few dozen to over one hundred genes. Similarly, molecular analyses of hematopoietic and solid malignancies continue to identify a growing number of genes relevant to diagnostic stratification, prognosis, and therapeutic response [42–44]. The technical complexity of multiple gene analysis is very challenging with the Sanger sequencing, whereas it can be addressed by targeted enrichment of multiple genes followed by NGS.

In general, enrichment strategies can be categorized as either amplification-based or oligonucleotide-array-capture-based [45, 46]. For amplification-based enrichment, single locus or multiplex PCR remains a mainstay method. Resequencing of exons by targeted PCR enrichment has been demonstrated for a number of genes wherein amplicons are generated from a sample, pooled in equimolar ratios, then ligated to adapter oligonucleotides to generate an NGS library (without overlapping fragments). This approach integrates well with the longer read length technologies (e.g., Roche 454). Another strategy for target enrichment is to use overlapping long-range PCR with 5–10 kilobase amplicons, which are then pooled, fragmented, and converted into a library. Multiple samples can be sequenced together if index or barcode sequences are incorporated into adapter sequences. After sequencing, indexed reads are assigned to their respective sample of origin by bioinformatic algorithms, which is sometimes referred to as “deconvolution” [47–49]. To automate PCR amplification of many targets, a highly parallel PCR microdroplet technology has been commercialized by RainDance Technologies (Lexington, MA) [50, 51]. In this technology, individual primer pairs for targets are designed (for amplicon lengths of 200–600 bp), synthesized, and sequestered in individual, stable, emulsion microdroplets. Up to a few thousand primer pairs can constitute a microdroplet primer population. For targeted amplification, genomic DNA is first fragmented into a size range of 2,000–4,000 bp and then randomly distributed into a separate microdroplet population. On an automated microfluidic platform, individual primer pair and genomic DNA microdroplets are merged by exposure to a voltage potential so that one microdroplet containing fragmented genomic DNA is associated with one microdroplet containing an individual primer pair. The fused microdroplets are collected into a single microfuge tube and thermocycled to achieve highly parallel PCR amplification. After amplification, the droplets are disrupted and the pooled amplification products are converted into a library for NGS. An additional PCR-based enrichment approach, Fluidigm (South San Francisco, CA), commercialized their Access Array platform based on a microfluidic chip that contains nanoliter-scale reaction chambers separated by valves. Forty-eight samples can be loaded onto the chip, and each sample can be distributed into 48 chambers with unique or multiplex primer pairs for subsequent PCR. After amplification, valves are reversed and samples are returned to their original wells and processed for NGS.

Oligonucleotide array capture methods constitute the second major enrichment strategy [52–59]. With array methods, genomic DNA, or genomic DNA converted into an NGS library, is hybridized to oligonucleotides complementary to target regions of interest. After hybridization, the enriched material is eluted from the array and processed for NGS. Originally, capture oligonucleotides were formatted on solid surface arrays, but this approach has been supplanted by in-solution formats offered by Roche NimbleGen (Madison, WI), Agilent (Santa Clara, CA), Illumina, and most recently Integrated DNA Technologies (IDT, Coralville, IA). In-solution oligonucleotide array capture procedures require several hours to days of hybridization depending on methodology, and an important array design consideration is the specificity of capture probes and the potential for co-capture and enrichment of nontarget sequences, notably from closely related genes and pseudogene analogs which can lead to false-positive variants in alignments. The inherent coverage variability in NGS libraries can be compounded by uneven and inadequate capture using in-solution probes, sometimes necessitating oligonucleotide array redesign.

A novel hybrid between capture and amplification approaches termed HaloPlex (Halo Genomics, Uppsala, Sweden and acquired by Agilent Technologies) is being adopted for clinical use. In this method, genomic DNA is digested with different combinations of restriction endonucleases and then hybridized to biotinylated oligonucleotide probes whose right and left sides are complementary to the fragment ends. The probes also contain a sequencing motif and index. Probe hybridization yields a genomic DNA circular structure that is isolated by incubation with magnetic streptavidin beads and then closed by ligation. Ligated circles are amplified by PCR yielding a library ready for sequencing.

Each enrichment method has advantages and disadvantages. For clinical applications, the choice of enrichment technology will depend on target size, content, and laboratory workflow preferences. Array capture is better suited for enrichment of larger target regions up to the scale of the human exome but is prone to co-capture of highly homologous sequences. Although amplification-based strategies are not as scalable, their increased target enrichment specificity may be preferable or required depending on target content. As an example of in-solution array capture data, Fig. 59.6 shows two genes, FOXP3 and HLA–DOB, and their respective capture results. The left panel shows capture probes and coverage of seven exons from the FOXP3 gene located on the X chromosome. The right panel shows capture probes and coverage of two of four exons from the HLA–DOB gene located on chromosome 6. HLA–DOB belongs to a large gene family with highly homologous sequences, making the design of unique capture probes difficult for all exons.

Currently, the majority of clinical NGS testing is comprised of multigene panels for inherited disorders and oncology applications (discussed further in the Considerations for Clinical Implementation of NGS). The adoption of multigene panels has been greatly facilitated by lower cost platforms, notably the Illumina MiSeq and the Ion Torrent PGM. In parallel, exome sequencing (and to a lesser degree genome sequencing) is being increasingly applied to identify causal and candidate genes in patients with disorders suggestive of a genetic etiology. We next discuss exome and genome sequencing, and the Illumina technology is highlighted as it has been used in the majority of clinical exome and genome publications.

In the scenario of a suspected inherited disorder, the overarching goal of exome or genome sequencing is to generate a list of variants from which a presumptively rare and pathogenic variant (or variants depending on mode of inheritance) can be identified. The complexity of this goal is reduced if a family with more than one affected member or even a trio of both parents and the affected child can be tested, as this allows segregation analyses of putative pathologic variants to be performed. While complex, the reported success of exome and genome sequencing for identifying candidate and causal genes in patients and families whose disorders strongly suggest a genetic etiology is 25–40 % [60–63]. Technical and bioinformatic considerations are presented before describing how to “search” for candidate genes.

As shown in Fig. 59.7, initial preparatory steps are the same for exome and genome libraries and include genomic DNA fragmentation and conversion of fragments into an oligonucleotide adapter-tagged library. For exome sequencing, in-solution capture with probes complementary to coding regions is performed as described above to yield an exome-enriched library. Several commercial vendors offer in-solution exome capture probe reagents including Agilent which uses RNA probes and Roche NimbleGen and Illumina which use single-stranded DNA oligonucleotide probes. While there is considerable overlap in the performance characteristics of these three approaches, each differs in targeted capture areas and capture probe sequence composition [64–66].

Prior to Illumina sequencing, the genome or exome library concentration needs to be determined to guide the preparation of an appropriate dilution for loading onto the flow cell. Quantitative PCR (qPCR) is performed on the library using amplification primers complementary to the adapter oligonucleotide sequences that bind to the flow cell surface capture oligonucleotides. A qPCR standard curve is generated using artificial templates of known concentrations that also contain the same primer target sequences complementary to the flow cell surface capture oligonucleotides. The library concentration is determined by comparing the PCR crossing threshold of the library to those of the standard curve templates. The use of PCR primers complementary to adapter sequences that bind to flow cell surface oligonucleotides specifically assesses the number of library molecules that should form clusters. Based on the library concentration determined by qPCR, a dilution is prepared with the goal of achieving optimal cluster density formation.

The bioinformatic steps for generating annotated variants from Illumina exome or genome sequencing data are similar to those described above. After sequencing, FASTQ read files are aligned to the human genome reference sequence using BWA followed by refined alignment and variant calling with GATK. Annotation is performed with either GATK or ANNOVAR. The process of alignment and variant calling generates, for a genome, approximately 3–3.5 million positions that differ from the human reference sequence. In comparison, approximately 15,000–30,000 differences from the reference will be observed in the coding regions in a genome or exome data set depending on bioinformatic filtering parameters. Within these ranges, the greatest number of variants from the human genome reference is observed in individuals of African and African-American descent, an intermediate number are observed in individuals of Asian descent, and the fewest variants are observed in Caucasians of European descent. At this time, genome sequencing costs approximately three times more than exome sequencing. Operationally, analyses of either genome or exome data typically focus on coding region variants and variants in close proximity to coding regions because this currently is the most “interpretable” portion of the genome. In considering the use of exome vs genome sequencing, exome sequencing does not provide as comprehensive an analysis of coding regions as genome sequencing. This reflects the difficulty of designing specific capture probes for all coding regions due to the variability of capture efficiency and the presence of cross homologous sequences. Figure 59.8a shows a comparative example of these differences for the RET gene in exome and genome data sets. The genome data shows comprehensive coverage across the gene region, whereas the exome data has sequencing reads only in locations enriched by capture probes. Figure 59.8b illustrates in exome sequencing of the ABCF1 gene that some exon regions do not have probes and therefore are not covered (Fig. 59.8b), while the genome sequencing data is complete across this region. Another technical consideration of exome sequencing that impacts data analysis is that coding regions with high or low GC content are captured less efficiently under single temperature hybridization conditions. This variability of capture was described by Clark et al. [65] and is inherent to all three commercially available exome capture reagents. The effect of GC content is particularly pronounced for the first exon of human genes, which on average have higher GC content. This can lead to reduced coverage of first exons in exome data sets as illustrated for the MAZ gene in Fig. 59.9 where read coverage for exon 1 is greatly reduced compared to other exons although a capture probe is included for this exon. While these limitations exist, exome sequencing is routinely performed at higher coverage depths compared to genome sequencing (due to cost considerations). The increased coverage depth can allow for identification of coding region variants missed in lower coverage genome sequencing data [26, 65].