Sequencing Quality Control

Go Back

06/07/2017

Figure 1: Illumina sequencing process. DNA fragments bind the flow cell through the adaptors. Thereafter, bridge amplification with unlabeled nucleotide creates the clusters. Nucleotides labeled with different dyes are incorporated and detected with a CCD camera that gives the base call. Figure modified from ThermoFisher scientific.

In the course of the last decade, genomic research has involved the use of massively parallel sequencing technologies in order to achieve deep sequencing coverage, fast turnaround times, high efficiency and resolution, all at relatively low cost per base. Unlike the Sanger sequencing technique, Next generation sequencing (NGS) platforms can perform the sequencing of millions of DNA molecules simultaneously generating massive amounts of sequence data1. Among these platforms, the Illumina platforms are the most widely used for genome sequencing with the least error rate per base2.

The chemistry behind the Illumina platforms adopts the strategy of “reversible terminator sequencing”. The DNA template strand is immobilized onto the flow cell containing oligos which are complementary to the adaptors ligated to the DNA strand. Each fragment is then copied by bridge amplification3,4 forming clusters with thousands of identical DNA sequences. The clustering process can occur either directly on the platform or separately on the cBot, depending upon the chosen Illumina system. Each cluster is sequenced using fluorescently-labelled, reversible terminator nucleotides. With the addition of dNTPs, fluorescence is emitted and imaged. Before the next chemistry cycle proceeds, the blocked 3’ end and the fluorophore from each incorporated base is removed, to allow incorporation of the next base. Since all four reversible terminator-bound dNTPs are present during each sequencing cycle, natural competition minimises incorporation bias (Figure 1). The final result is obtained through a base-by-base sequencing.

 

 

Quality control (QC) procedures

As part of the routine at the OGC, quality checks of all the sequencing runs are carried out daily in order to guarantee the best performance for each processed project. The Quality Controls (QC) are namely accomplished through three main steps:

  1. Monitoring the run and initial evaluation.
  2. Preliminary assessment of the pool quality.
  3. Main QC summary statistics.

 

Monitoring the run and initial evaluation

During the run, metrics are shown on the Sequencing Analysis Viewer (SAV) interface (Figure 2, below). The SAV software provides an easy interface that allows monitoring of the main run parameters during and after the sequencing.

For the initial cycles the intensity values of the four incorporated bases can be assessed to determine if a clustering failure has occurred. No intensity (or very low intensities) would indicate a failure which can be due to issues with the library itself or reagent and technical problems (e.g. fluidic blockages, camera issues etc.).

The next metric to assess is the Cluster Density K/mm2 (represents the clusters for each tile in thousands per mm2). This value will be available on the SAV after cycle 5, 7 or 20 for HiSeq 2500, MiSeq (v3 kit) and MiSeq (v2 kit) platform, respectively. Successful results depend upon a correct prediction of how many clusters will be generated. An accurate prediction is dependent on the library type and method for quantifying. Overloading of the library may cause merging of clusters, and hence an increased density and lowering of the % PF with a reduction of the quality score (Q30). Conversely, a low loading input of library may result in low cluster density and insufficient data yield (Figure 3.1 and 3.2, below). Please note that this metric is not meaningful for the HiSeq 4000 platform (an ordered flow cell is utilised on this system so the Cluster Density K/mm2 is the same for all lanes).

After cycle 25 for all systems, some important further metrics for Read 1 will be available for quality control on the “summary” section of the SAV:

  • % PF (percentage of reads passing filter which reflect the chastity of the intensity signal). The purity of the signal from each cluster is examined over the first 25 cycles and calculated as “Chastity” of a base call. This value is reported as the ratio of the brightest intensity divided by the sum of the brightest and second brightest intensities for each cycle. Optimal values for this parameter depends upon the platform. For MiSeq and HiSeq 2500, % PF varies from 80-95% depending on such factors as base complexity, cluster density etc.; for HiSeq 4000, % PF fluctuates from 60-80%.
  • % phasing and pre-phasing (percentages of bases that fell behind or jumped ahead the current cycle within a read, respectively). The proportion of sequences in each cluster that are affected by phasing and pre-phasing increases with cycle number, hampering correct base identification for long reads.5,6 Optimal values for these parameters are usually below 0.5 or 0.2% depending on platform. Issues may arise when read length is above Illumina specifications, paired end reagents are improperly mixed, cross-contamination between SBS reagents, fluidics issue, etc.

  • % Q-score >= Q30 (percentage of bases that have a Q-score above or equal to 30; Q30 is a probability of incorrect base calling of 1 in 1000).
  • % align (percentage of the reads that aligned to the PhiX control genome). The presence of the PhiX is important in order to estimate the error rate (see below). Depending upon the platform, the % align may provide an indication on the success of the clustering process, e.g. higher values may signify underclustering or imprecise cluster density prediction. However, this is not always true for all platforms.
  • Error rate (calculated for each read based on the aligned % of PhiX control sample).

It is important to monitor the % Q-score >= Q30 and Error rate throughout the run to determine that the sequencing is proceeding as expected.

 

Preliminary assessment of the pool quality

Most runs performed use a library which contains a number of different samples pooled together. This approach is very common due to the existence of specific combinations of oligos (or tags) which can be used to uniquely identify samples. The samples are generally pooled equimolarly for each tag within a pool. During the run (after the index reads are completed), a report provides information about the counts of the index tags following “de-multiplexing”. This is performed by our bioinformatics pipeline for a single tile. Loss of data yield is expected to vary from 3-10% depending on library type and platform. However, loss of data can be higher if the indexing strategy is suboptimal, contamination occurred with other libraries, etc. The report also provides a bar chart in order to visually assess the success (balance) of pooling.

 

Main QC summary

A pre-installed software – the Real-time analysis (RTA) software – performs the primary analysis on all the Illumina platforms during the sequencing run. The RTA collects the temporary image files at every cycle from the machine and extracts the intensities by detecting the signal of the generated clusters. The intensity files are subsequently converted to base calls files (.bcf). During the secondary analysis, the base calls files are uploaded on the servers located at the WTCHG, and de-multiplexed into “fastq” files. The fastq files contain information about the sequence of the read and Q-scores for each base in all reads. At last, the fastq files will be converted to alignment files (BAM and bai) during mapping to the reference genome and generate the main QC summary statistics. As standard, either Stampy in combination with BWA7, or BWA mem, is used to map reads. The main QC file is a compendium of explanatory graphs and tables containing information for all the mapped or unmapped reads that are fundamental to evaluate the overall quality of a run. All the information is reported in graphs and tables containing values related to each pool belonging to a specific lane of the sequenced flow cell.

A table of the most important QC metrics and their meaning can be found below:

Figure 2: Overview of the Sequencing Analysis Software (SAV) for a HiSeq 4000 example run. In the “Flow Cell Chart”, on the left, is shown a schematic flow cell with the lanes, tiles, and swaths. It is possible to observe the distribution of different parameters (intensity, Q30, Error rate, Density PF) choosing from the dropdown list. This schematic view of the flow cell is very useful to check if any specific area of the flow cell is affected. The “Data By Cycle”, in the middle, refers to a graph in which different parameters chosen from the dropdown list (e.g. intensities, percentage of bases and, the Q30) can be inspected lane by lane per each cycle. The “Data by Lane”, below, shows the density of reads that passed the filter against the raw data. The “QScore Distribution” graph shows the percentage of bases with a Q-Score above 30 which is logarithmically correlated to the base calling error probability. Therefore, a Q-score of 30 (Q30) indicate the probability of inserting one incorrect base every 1,000 times (base call accuracy of 99.9%).

Figure 3: (3.1) Flow cell charts evaluation for overclustering. A) & B) Significant overclustering on HiSeq 2500 in Rapid run mode, C) on HiSeq 2500 High Output mode and, D) on MiSeq. (3.2) Thumbnail images showing cluster densities gradient from underclustered to overclustered flow cell. The above images can be investigated in the “Imaging” tab of the SAV software. The images will reflect numerical values shown in the “imaging metrics table” that can be found in the same section. (Image modified from Illumina).

 

 

 

Most important QC metrics and their meanings
G+C histogram Describes the GC distribution of the reads (the graph will depend on the genome and library type being sequencing). The dotted line refers to all reads that have passed chastity filters whereas solid lines refers to only mapped reads.
Insert size histogram The insert size distribution is summarized by the median and median absolute deviation. Note that this is for mapped reads and therefore may give questionable results if the % mapping is low.
Genomic coverage by G+C The G+C fraction is computed from the reference genome, over the approximate fragment regions with coverage in the top 0.1 percentile being excluded. The dotted line shows the G+C histogram of the reference. This graph is useful to compare these in order to inspect for GC skewing of the library.
Fraction N/lowQ, read 1 (read 2) Shows the fraction of N bases against the cycle number. Multiply the peak height by 100 to obtain the percentage of Ns at any given cycle. Please note that the scale of the X-axis varies substantially on this graph. Spikes of N bases below 2% are common and should be ignored and occasional spikes of 25-50% are expected.
G+C by cycle (PF), read 1 (read 2) This graph can be used to see if there are any unexpected complexity issues or large scale GC skews at specific cycles.
Mean Q by cycle, read 1 (read 2) The solid line is the numerical mean Q score and is a measure of the average information content per read. This is calculated using mapped reads only, whereas the dashed line in this plot uses all reads. This graph is useful for assessing the quality of sequencing throughout the read.
Q score histogram, read 1 (read 2) Describes the dotted line information in the Mean Q by cycle graph as a histogram.
Yield (Mrd) This describes the Yield in Mrd (million reads) for all the libraries in the pool.
Mean G+C Shows the the overall GC% for each of the libraries. The colour blue indicates the value for read 1 and blue indicates read 2.

Authors: Maria Lopopolo and Lorne Lonie

 


References

  1. Metzker ML: Sequencing technologies – the next generation. Nat Rev Genet 2010, 11(1):31–46.
  2. Xi Yang1, Di Liu1, Fei Liu1, Jun Wu1, Jing Zou1, Xue Xiao1, Fangqing Zhao2 and Baoli Zhu1 HTQC: a fast quality control toolkit for Illumina. BMC Bioinformatics 2013, 14:33 http://www.biomedcentral.com/1471-2105/14/33.
  3. Bennett S: Solexa Ltd. Pharmacogenomics 2004; 5: 433-8.
  4. Bentley DR, Balasubramanian S, Swerdlow HP, Smith GP, Milton J, Brown CG, Hall KP, Evers DJ, Barnes CL, Bignell HR, Boutell JM, Bryant J, Carter RJ, Keira Cheetham R, Cox AJ, Ellis DJ, Flatbush MR, Gormley NA, Humphray SJ, Irving LJ, Karbelashvili MS, Kirk SM, Li H, Liu X, Maisinger KS, Murray LJ, Obradovic B, Ost T, Parkinson ML, Pratt MR, Rasolonjatovo IM, Reed MT, Rigatti R, Rodighiero C, Ross MT, Sabot A, Sankar SV, Scally A, Schroth GP, Smith ME, Smith VP, Spiridou A, Torrance PE, Tzonev SS, Vermaas EH, Walter K, Wu X, Zhang L, Alam MD, Anastasi C, Aniebo IC, Bailey DM, Bancarz IR, Banerjee S, Barbour SG, Baybayan PA, Benoit VA, Benson KF, Bevis C, Black PJ, Boodhun A, Brennan JS, Bridgham JA, Brown RC, Brown AA, Buermann DH, Bundu AA, Burrows JC, Carter NP, Castillo N, Chiara E Catenazzi M, Chang S, Neil Cooley R, Crake NR, Dada OO, Diakoumakos KD, Dominguez-Fernandez B, Earnshaw DJ, Egbujor UC, Elmore DW, Etchin SS, Ewan MR, Fedurco M, Fraser LJ, Fuentes Fajardo KV, Scott Furey W, George D, Gietzen KJ, Goddard CP, Golda GS, Granieri PA, Green DE, Gustafson DL, Hansen NF, Harnish K, Haudenschild CD, Heyer NI, Hims MM, Ho JT, Horgan AM, Hoschler K, Hurwitz S, Ivanov DV, Johnson MQ, James T, Huw Jones TA, Kang GD, Kerelska TH, Kersey AD, Khrebtukova I, Kindwall AP, Kingsbury Z, Kokko-Gonzales PI, Kumar A, Laurent MA, Lawley CT, Lee SE, Lee X, Liao AK, Loch JA, Lok M, Luo S, Mammen RM, Martin JW, McCauley PG, McNitt P, Mehta P, Moon KW, Mullens JW, Newington T, Ning Z, Ling Ng B, Novo SM, O’Neill MJ, Osborne MA, Osnowski A, Ostadan O, Paraschos LL, Pickering L, Pike AC, Pike AC, Chris Pinkard D, Pliskin DP, Podhasky J, Quijano VJ, Raczy C, Rae VH, Rawlings SR, Chiva Rodriguez A, Roe PM, Rogers J, Rogert Bacigalupo MC, Romanov N, Romieu A, Roth RK, Rourke NJ, Ruediger ST, Rusman E, Sanches-Kuiper RM, Schenker MR, Seoane JM, Shaw RJ, Shiver MK, Short SW, Sizto NL, Sluis JP, Smith MA, Ernest Sohna Sohna J, Spence EJ, Stevens K, Sutton N, Szajkowski L, Tregidgo CL, Turcatti G, Vandevondele S, Verhovsky Y, Virk SM, Wakelin S, Walcott GC, Wang J, Worsley GJ, Yan J, Yau L, Zuerlein M, Rogers J, Mullikin JC, Hurles ME, McCooke NJ, West JS, Oaks FL, Lundberg PL, Klenerman D, Durbin R, Smith AJ: Accurate whole human genome sequencing using reversible terminator chemistry. Nature 2008; 456: 53-9.
  5. Erlich Y, Mitra PP, delaBastide M, McCombie WR, Hannon GJ: Alta-Cyclic: a self-optimizing base caller for next-generation sequencing. Nat Methods. 2008, 5: 679-682. 10.1038/nmeth.1230.
  6. Kircher M1, Stenzel U, Kelso J. Improved base calling for the Illumina Genome Analyzer using machine learning strategies. Genome Biol. 2009;10(8):R83. doi: 10.1186/gb-2009-10-8-r83. Epub 2009 Aug 14.
  7. Lunter G, Goodson M. Stampy: a statistical algorithm for sensitive and fast mapping of Illumina sequence reads. Genome Res. 2011 Jun;21(6):936-9