Patent Application
20120191364
The present invention relates to a method of DNA sequencing and in particular but not exclusively to methods and arrays for nucleotide base calling.
Every year there is an exponential growth in the amount of DNA sequence information generated and deposited into Genbank. Many of the current sequencing technologies use a form of sequencing by synthesis (SBS), wherein specially designed nucleotides and DNA polymerases are used to read the sequence of chip-bound, single-stranded DNA templates in a controlled manner. To attain high throughput, many millions of such template spots are arrayed across a sequencing chip and their sequence is independently read out and recorded. Devices, equations, and computer systems for making and using arrays of material on a substrate for DNA sequencing are known. However, there is a continued need for methods and compositions for increasing the fidelity and accuracy of sequencing nucleic acid sequences.
Sequencing of viral genomes in particular has historically been performed using standard dye termination technologies. In recent years, many researchers have migrated away from traditional capillary sequencing instruments and towards high-throughput DNA sequencing technologies that provide higher accuracy at a lower cost. However, these technologies are still too slow, costly and labour-intensive to obtain genomic sequences of viruses that mutate ever so frequently and for large-scale epidemiologic or evolutionary investigations in viral outbreaks. For example, the currently available sequencing technology is not suitable for sequencing the genomic sequences of H1NA influenza A virus and in particular the 2009 influenza A (H1N1) virus from the ever-increasing pool of infected individuals.
In April 2009, a novel swine-origin H1N1 influenza A virus erupted in Mexico and spread swiftly across the world at unprecedented speed, forcing the World Health Organization (WHO) to raise its pandemic alert to phase 5. As of September 13th, WHO had reported over 2,96,471 laboratory-confirmed cases of pandemic (H1N1) 2009 in 135 countries. However, these figures are likely to be an underestimate as surveillance has been focused on severe cases. Fortunately, despite the high transmissibility of this outbreak, there has been a low number of fatalities (3,486 reported deaths). This suggests that the virulence of the 2009 influenza A (H1N1) virus may be relatively low.
The influenza pandemics of 1918, 1957, and 1968 that killed millions of people remind us that the most recent 2009 influenza A (H1N1) virus outbreak should not be taken lightly. This virus will continue to evolve through mutations and/or recombination that may increase its virulence and/or drug resistance of the virus. As drug companies rush to supply the world with antiviral drugs for this pandemic outbreak, isolated cases of drug-resistant H1N1 flu strains have already emerged. These drug-resistant strains usually have mutations near drug-binding sites that reduce the binding affinities and effectiveness of certain drugs. Thus, it is absolutely vital that the evolution of the 2009 influenza A(H1N1) viruses be closely and continually monitored for any genetic variations.
Oligonucleotide resequencing microarrays that are capable of identifying nucleotide sequence variants may offer an alternative solution to the standard dye termination technologies and in recent years, have been used for detecting and subtyping influenza viruses. By analysing sequences generated from tiling probes across targeted regions of various strains of the influenza virus (e.g. partial fragments of the haemagglutinin (HA) and neuraminidase (NA) genes), important information such as viral subtypes, lineages and sequence variants can be determined. Analysis of the sequences is usually done using platform accompanying software that employs probabilistic base-calling algorithms such as ABACUS and Nimblescan PBC. Although statistically sound, these methods are susceptible to hybridization noise caused by factors such as poor probe quality, poor amplification or mutations. This results in numerous ambiguous and false positive base calls that may affect the accuracy of downstream evolutionary analysis. Efforts have been made to improve the call rates and accuracies of existing probabilistic base-calling algorithms but the methods mostly result in the base call rates suffering.
Also, ideally during sequencing, a perfect match (PM) probe used in the sequencing, would be expected to gain a hybridization intensity multi-fold that of its corresponding mismatch (MM) probes, making base calling a straight-forward task. However, two types of errors are prevalent in practice:
A myriad of factors, such as weak PCR products, suboptimal annealing temperatures, CG biases, poor probe quality, and non-specific binding of MM probes have been attributed to be the causes of these two types of errors. With the use of better primers, optimization of annealing temperatures and the use of variable length probes, certain factors such as weak PCR products and CG biases can be overcome. However, some factors are unavoidable. This implies that even under optimal experimental conditions, there may still exists MM probes that do not exhibit a significant reduction in hybridization intensity relative to the PM probe, causing a type I error. The tiling requirement of a resequencing array also greatly inhibits the exclusion of poor quality probes from the array. For example, the inclusion of probes that are of low complexity or containing consecutive runs of the same nucleotide (homopolymers) are likely to cause type II errors since they have a higher tendency to exhibit non-specific cross-hybridization.
These factors affect the hybridization intensities of the PM/MM probes has proved useful in designing probes for microarray experiments however, the accuracy of sequence calling has yet to be improved.
The present invention is defined in the appended independent claim. Some optional features of the present invention are defined in the appended dependent claims.
In general terms, the invention sequencing a first polynucleotide strand (e.g. a strand of a virus which is believed to have mutated) using the known polynucleotide structure of a second polynucleotide strand (e.g. the virus before mutation). For each of a number of fragments of the second polynucleotide strand, and for each position along each fragment, we obtain (i) “first probe data” describing the hybridization activity of the first polynucleotide strand with a “first probe” designed to bind with a portion of the second polynucleotide strand centred at that position, and (ii) “second probe data” describing the hybridization of the first polynucleotide strand with “second probes” which differ from the first probe only at that position. In positions where the hybridization with the first probe is much greater than with the second probe, it is likely that the first and second polynucleotides are the same. In other positions, there is a higher chance of a mutation.
In one specific expression, the present invention relates to a method of sequencing a first polynucleotide strand comprising a first polynucleotide sequence, the first polynucleotide strand resembling a second polynucleotide strand having a known second polynucleotide sequence, the method employing a data set which, for one or more fragments of the second polynucleotide sequence, contains:
The method of the present invention may enable large-scale identification of variations in polynucleotide sequences. In particular, it may enable large-scale identification of variations in viruses. This may be advantageous especially with H1N1 (2009) viruses which mutate easily and frequently and may vary in multiple patient samples. The method of the present invention may provide a means for rapidly whole-genome sequencing the H1N1 samples.
The term “fragment” is used here to refer to a part (i.e. a sub-set) of the second polynucleotide strand, with no implication that the fragment has been separated from the rest of the second polynucleotide strand. Preferably the set of fragments collectively span the entire second polynucleotide strand (in the sense that every base in the second polynucleotide strand is included within at least one of the fragments), so that if the first polynucleotide strand differs from the second polynucleotide strand only by mutations, the method may be used to sequence substantially the whole of the first polynucleotide strand (also, in some instances, as discussed below, at certain isolated positions, the method may determine that no identification of the base is possible). Alternatively, the fragments may be selected such that they do not span the entire second polynucleotide strand (e.g. to omit portions of the polynucleotide strand which are not believed to be of clinical importance).
The first probe is “designed to bind to a portion of the second polynucleotide strand” in the sense of having a sequence complementary to that portion of the second polynucleotide strand.
The one of the first and second probes which is complementary to the first nucleotide strand at the central position (i.e. the probe with the highest hybridization, activity) is called the “perfect match probe”, and the other probes are called “mismatch probes”. In the case that the corresponding portion of the first polynucleotide strand does not contain a mutation, the “first probe” is the “perfect match probe”, and the second probes are the mismatch probes. Conversely, if there is a mutation at the central position, then the corresponding one of the second probes is the “perfect match probe”, and the first probe and the other second probes are the mismatch probes.
In one embodiment, the method further comprises at each said position,
The said at least one second numerical parameter for each said position may include a parameter comparing the mean and the standard deviation of the corresponding first probe data and second probe data. If either of said determinations is negative, a verification algorithm may be performed using data (“perfect match data”) describing the hybridization intensity of the perfect match probe of neighbouring positions.
The verification algorithm may comprise a first determination of whether the perfect match data for the neighbouring positions is indicative of a divergence between the first and second nucleotide sequences at said position. The first determination may be positive if the average of the perfect match data for one or more nearest neighbouring positions is lower than the perfect match data for neighbouring positions further from said position than said nearest neighboring positions.
Alternatively or additionally, the verification algorithm may comprise a second determination of whether there is a likelihood of a substitution bias at said position. One of said second numerical parameters may be obtained from the hybridization intensity-based order of the PM probe and mismatch probes for the site. Suppose that, for a given position, we say that a given probe encodes base b if b is located at the centre of the region. We denote the base encoded by the PM probe as b1 and the mismatch probes encode b2, b3 and b4 where {b1, b2, b3, b4}={A, C, G, T}. Without loss of generality, we will assume that hybridization intensity reduction order is b1b2b3, b4. The second numerical parameter may then be obtained as a ratio fobs/frand, where fobs is a probability of observing the hybridization intensity reduction order b1b2b3b4 given that the perfect match probe encodes b1, and frand, is the probability of observing the hybridization intensity reduction order b1b2b3b4 by chance.
The values fobs and frand may be obtained by calculating:
wherein, for any order of the bases denoted by wxyz, the function #(wxyz) denotes the number of times, in a number t of other positions, that the hybridization intensity reduction order was wxyz. Preferably the t positions are those in which the first numerical parameter indicated that the first and second nucleotide strands were both b1, and #(wx) denotes the number of times, in the t positions that the hybridization order began wx. For example, #(b1b2)=#(b1b2 b3b4)+#(b1b2 b4 b3).
Upon said first determination being positive and said second determination being negative, it may be determined that the nucleic acid of the first polynucleotide sequence differs from the nucleic acid of the second polynucleotide sequence at said position.
In another specific expression, the present invention relates to a method of sequencing a pair of first polynucleotide strands, which are complementary strands having complementary first polynucleotide sequences. In particular, in, the pair of strands, one strand has the first polynucleotide sequence and the other strand has a polynucleotide sequence complementary to the first polynucleotide sequence. The method comprises performing a method according to any aspect of the present invention for each first polynucleotide strand using a respective second polynucleotide strand, the second polynucleotide strand having complementary respective second polynucleotide sequence, for each corresponding position in the second polynucleotide sequence, said verification algorithm may be performed upon a determination that said first numerical parameters are indicative of the two first polynucleotide sequences not being complementary in that position.
As mentioned above, the set of fragments of the second polynucleotide sequence may collectively span the entire polynucleotide strand. Preferably, the fragments overlap to some degree, so that the dataset contains multiple sets of perfect match data and mismatch data for locations in the overlap regions. This data may be averaged before calculating the first numerical parameter in respect of such positions. Preferably, the overlap regions are selected to include regions considers to be critical in the sense given below, so that more accurate sequencing of the critical regions is possible.
In one expression, the present invention relates to a method of producing an array for sequencing a first polynucleotide strand having a first polynucleotide sequence, the method employing data encoding a second polynucleotide sequence of a polynucleotide strand resembling the first polynucleotide strand, the method comprising:
Step (a) of defining the one or more fragments may include:
The method above may be implemented by a computer (e.g. any general purpose computer, such as a PC) having a processor and a data storage device containing program instructions operable by the processor to carry out the method. Furthermore, a computer program product (e.g. a software download, or a tangible data storage device, such as a CD-ROM) may be provided containing such program instructions.
In another expression, the present invention relates to an array for sequencing a first polynucleotide strand having a first polynucleotide sequence and resembling a second polynucleotide strand having a second, known polynucleotide sequence, the array comprising, for each of one or more fragments of the second polynucleotide sequence:
These arrays may be used as a practical, large-scale re-sequencing tool. Also, the sequences obtained from the arrays may also be highly reproducible.
The dataset may be derived using an array which may be produced by a method according to any aspect of the present invention and/or an array according to any aspect of the present invention.
The second polynucleotide strand may be a RNA or DNA of a virus. In particular, the virus may be influenza A virus. More in particular, the virus may be H1N1 influenza A virus.
In another expression, the present invention relates to a kit comprising:
Preferably, the computer readable medium may be fully-automated and may provide a comprehensive graphical report that shows the first polynucleotide sequence quality and the location of all mutations with their associated confidence and proximity to the important regions in the first polynucleotide strand. The short turnaround time from sample to sequence and analysis results may also be short. For example, it may take approximately 30 hours for 24 samples, making this kit an efficient large-scale evolutionary surveillance tool.
The array may be a 12-plex array. The kit may be used for sequencing H1N1 influenza A virus. In particular, the H1N1 influenza A virus may be 2009 influenza A(H1N1) virus. More in particular, the computer readable medium may be used for automatic base-calling and variant analysis, capable of interrogating all eight segments of the 2009 influenza. A(H1N1) virus genome and its variants. The array according to any aspect of the present invention may be able to detect all sequence variations with respect to a second polynucleotide strand with a second polynucleotide sequence. In particular, the second polynucleotide sequence may be a consensus 2009 influenza A(H1N1) virus sequences with added focus on important regions such as drug-binding sites, structural components and previously reported mutations.
The consensus 2009 influenza A (H1N1) may comprise at least one sequence selected from the group consisting of SEQ ID NO:1 to SEQ ID NO:8, fragment(s), derivative(s), mutation(s), and complementary sequence(s) thereof. In particular, the consensus 2009 influenza A (H1N1) may consists of nucleotide sequences SEQ ID NO:1 to SEQ ID NO:8.
In another expression, the present invention relates to isolated oligonucleotide comprising at least one nucleotide sequence selected from the group consisting of: SEQ ID NO:1 to SEQ ID NO:8, fragment(s), derivative(s), mutation(s), and complementary sequence(s) thereof. The sequences may be derived from H1N1 influenza A.
As will be apparent from the following description, preferred embodiments of the present invention allow an optimal use of the method of the present invention to take advantage of the accuracy, speed and reproducibility. This and other related advantages will be apparent to skilled persons from the description below.
Preferred embodiments of a method of DNA sequencing will now be described by way of example with reference to the accompanying figures in which:
The term, “resembling” is used herein to refer to a measure of similarity. In particular, it refers to the measure of similarity between the first polynucleotide strand and the second polynucleotide strand: For example, the polynucleotide sequence of the first strand may vary from the polynucleotide sequence of the second strand by 1-20 nucleotides. In particular, the polynucleotide sequence of the first strand may vary from that of the second strand by 1, 2, 3, 4, 5, 10 or 15 nucleotides. The polynucleotide sequence of the first strand may be 95-99% similar to the polynucleotide sequence of the second strand.
The term “fragment” is used herein to refer to a portion of the second polynucleotide strand. In particular, the fragment may refer to a sequence of the polynucleotide that is at least 5 nucleotides long. More in particular, the fragment may refer to a sequence of the second polynucleotide strand that is 5, 8, 10, 15, 20, 25, or 25 nucleotides long. It may also refer to a longer fragment, such as an entire segment of the virus, and thus be up to several hundred or thousand nucleotides long.
The term “second polynucleotide strand” is used herein to refer to a reference sequence or part thereof. The second polynucleotide strand may be a consensus sequence and/or a known sequence used as a reference to determine the polynucleotide sequence of the first nucleotide strand.
The term “nucleic acid” is used herein to includes, but is not limited to, a monomer that includes a base linked to a sugar, such as a pyrimidine, purine or synthetic analogs thereof, or a base linked to an amino acid, as in a peptide nucleic acid (PNA). A nucleotide is one monomer in a polynucleotide. A nucleotide sequence refers to the sequence of bases in a polynucleotide.
The term “polynucleotide” is used herein to refer to a nucleic acid sequence (such as a linear sequence) of any length. Therefore, a polynucleotide includes oligonucleotides, and also gene sequences found in chromosomes. The term “polynucleotide” also encompassed RNA or DNA, as well as mRNA and cDNA corresponding to or complementary to the RNA or DNA. A fragment of a polynucleotide is a shortened length of the polynucleotide.
The term “mutation” of a position in the first polynucleotide sequence, refers at least one nucleic acid that varies from at least one reference (second) sequence via substitution, deletion or addition of at least one nucleic acid. In particular, the mutants may be naturally occurring or may be recombinantly or synthetically produced.
This method of sequencing is a platform-independent automated method for sequence calling that analyzes data from results of any array. The method adopts a gain-of-signal approach which assumes that the signal intensity of the perfect match (PM) probe (which matches exactly to the polynucleotide sequence in a sample) will be significantly higher than that of the corresponding mismatch (MM) probes. Hence, base calls are made by quantifying the gain in hybridization intensities of a PM probe over its corresponding MM probes. Using this method, an indication of the type of error in a suspicious base call is determined and the true PM probe may be discerned from the noisy MM probes.
The flowchart of the two-step process for base-calling is shown in
The remaining bases (i.e. Base queries with hybridization intensity abnormalities) are then passed to step 2 of
If not, the method passes to a sequence correction step.
The terms “base query” and “query base” are interchangeably used and are herein used to refer to a nucleic acid in a sequence that is not known and/or shows signs of hybridization intensity abnormalities. The base query refers to a position in the first polynucleotide strand that is to be determined using the method according to any aspect of the present invention.
All base queries with type I or II errors are assumed to have the following characteristics:
1. The base derived from the PM probe in the forward strand is not the same as the base derived from the PM probe in the reverse strand,
2. In either or both of the forward or reverse strands, the putative PM probe (the probe with the highest hybridization intensity) does not have hybridization intensity significantly higher than that of its MM probes,
3. One or more of its eight querying probes at any one position have unusually low signal-to-noise ratio. For a probe, its signal-to-noise ratio is defined as the ratio of the mean to the standard deviation of the intensities of the 9 pixels on the array encoding the probe.
Under optimized experimental conditions, the average percentage of high confidence calls made per sample is approximately 93%. Thus the number of non-high confidence calls (7%) can still seriously undermine the reliability of sequences generated by an array. Thus, it is imperative that these problematic queries be identified and subjected to further analysis.
The second step specifically comprises mutation confirmation and recovery of unreliable base queries through: neighbourhood hybridization intensity profile (NHIP) analysis and nucleotide substitution bias analysis.
In step 2, to extract any information out of noisy base calls, and unreliable base calls and to obtain more assurances of putative mutation calls, hybridization intensity patterns are used. Since a high-confidence mutation call may be a result of coincidental non-specific hybridization of the same MM probe in both strands, it is important to validate the mutation.
Many factors that cause noise in resequencing arrays do not only affect a single isolated query base. For example, if a region of the sample sequence is not amplified efficiently by PCR, the query bases in the region will be erroneous. As another example, when a single nucleotide mutation occurs at a particular query base, it may affect the hybridization intensities of probes belonging to neighbouring query bases as well.
The nature of a suspicious query base is determined by analyzing the hybridization intensities of its PM and MM probes together with its neighbouring (±6 bases from query base) PM and MM probes. Collectively, the hybridization intensities of these probes form a NHIP of the query base. Each query base is analysed to be classified as an isolated error, part of a poor quality region or real sequence variation based on its NHIP.
NHIP analysis results in a more informative decision on base-calling. Five distinct types of NHIP belonging to true non-mutations (wild-type), true mutations, isolated errors/‘N’s, long consecutive errors/‘N’s, and unknown errors/‘N’s, respectively are present and shown in
On the other hand, query bases with NHIP shown in
Query bases with NHIP shown in
Finally, to confirm the mutation and/or to identify the nucleic acid at the base query, nucleotide substitution bias analysis is carried out on these query bases.
Viral RNA from diagnostic swabs or RNA extracted from MDCK cell cultures was extracted using the DNA minikit (Qiagen, Inc, Valencia, Calif., USA) according to manufacturer's instructions. RNA was reverse-transcribed to cDNA using customized random primers designed using LOMA (Lee, 2008) and then amplified by PCR using proprietary H1N1 (2009) specific primers. The presence of 2009 influenza A (H1N1) virus in the samples was confirmed using a separate real-time PCR assay based on the published primer sequences from the Centre for Disease Control and Prevention (CDC), USA.
36 mutation hotspots were found in the alignments where mutations occurred near one another (within 20 bp). A perfect match (PM) probe residing in a mutation hotspot may contain mismatches that will have a detrimental effect on its hybridization intensity. To avoid this problem, additional mismatch probes were designed that contain all possible combinations of mutations found in each mutation hotspot. Thus, if two mutations are found within 20 bp of each other in the alignments, then in total four (22) additional mismatch probes were needed to encode them. In general, 2x additional mismatch probes are needed to completely encode a cluster of x mutations that occur within 20 bp of one another in the alignments.
The 2009 Influenza A (H1N1) virus resequencing array was designed based on eight consensus sequences (one for each segment; SEQ ID NO:1-8) derived from 1715 complete and partial sequences of 2009 Influenza A (H1N1) virus isolates deposited in NLM/NCBI H1N1 flu resources database (http://www.ncbi.nlm.nih.gov/genomes/FLU/SwineFlu.html) as of Jun. 11, 2009. Each consensus sequence of a segment was generated by aligning all available sequences of the segment using MAFFT (Koh, 2008) with high accuracy option. At the time of production (June 2009), no deletions, insertions or significant evidence of recombination in the alignments of the eight segments were found. There has also been no reports of any deletions, insertions or recombination in 2009 Influenza A (H1N1) virus sequences deposited in NCBI up to September 2009. This suggests that, at the present stage, mutation is the only evolutionary mechanism driving changes to the 2009 Influenza A (H1N1) virus.
Probes encoding all possible combinations of such mutations (as mentioned in the Design of probes in mutation hotspots section, subject to the maximum probe limit of the array) were included. Lastly, to enhance the usability of the array not only as an evolutionary surveillance tool but also as an evolutionary alarm, genomic sequences of the drug-binding pocket targeted by neuraminidase inhibitors (Maurer-Stroh S, 2009) such as oseltamivir (Tamiflu®) and zanamivir (Relenza®) were included onto the array. In this way, any nucleotide mutations that might cause a change in the amino acids in the drug-binding pocket and consequently render current neuraminidase inhibitors ineffective, will be accurately detected and reported by the array.
The complete list of consensus sequences, mutational hotspots, structural important sites and drug-binding sites of the 2009 Influenza A (H1N1) virus used for the design of the array of the preferred embodiment is given in Table 1. The sequence of the 8 segments of the consensus sequence is in Table, 2. There are 54 sequences of total length 16,861 bases. In order to interrogate both strands of the 54 sequences for all possible single nucleotide substitutions, the array consists of 8×16,861 probes (of variable length 29-39 nucleotides with optimized annealing temperature). There are 4 probes (‘A’, ‘C’, ‘G’ and ‘T’ probes) to interrogate each base of the 54 sequences on each strand. Among these 4 probes, the one that matches exactly to the given sample sequence is known as the perfect match (PM) probe, while the rest are mismatch (MM) probes. The correct base is deduced by analyzing the differences in hybridization signal intensities between sequences that bind strongly to the PM probe and those that bind weakly to the corresponding MM probes. As such, probes are designed such that the location of the interrogated target base is in the centre-most position of the probe, and thus provides the best discrimination for hybridization specificity. The array design ensures that bases that reside in the important regions of the virus are queried at least 4 and up to 8 times each and at least 2 times otherwise, and provides 99.9 percent coverage of the 2009 Influenza A (H1N1) virus (dated June 2009).
Due to the small amount of virus present in samples relative to human or cell-line total RNA, it was necessary to amplify the viral RNA through PCR. A combination of sequence-specific and random PCR approaches using LOMA-optimized primers (Lee, 2008) were used. The addition of random primers ensured complete genome amplification, even if mutations were present at the specific-primer binding sites. PCR conditions were optimized by conducting five duplicate hybridizations of the same virus sample cultured from a patient sample under different PCR conditions. The optimized method was then tested on RNA isolated directly from nasal swabs obtained from the same patient and from virus grown in cell culture. Microarray sequences generated from these replicate experiments were compared with capillary sequencing to estimate sequencing accuracy. Results not shown.
Identification of Base Queries with Suspicion of Type I or II Errors (Step 1)
The array specifies that eight probes (four for the forward strand and four for the reverse strand) were used to query each base. For each probe, the hybridization intensity is given by the mean and standard deviation of the fluorescence intensities of 9 individually scanned pixels associated with the probe on the microarray.
The signal-to-noise ratio (SNR) of a probe is defined as the ratio of the mean to the standard deviation of the intensities of the nine pixels associated with the probe. >95% of all probes had SNR less than TSNR (TSNR=μ
Base queries with one or more probes with ≧TSNR are analysed further in step 2. All base queries whose PM probe in the forward strand and PM probe in the reverse strand are non-complementary, or have weak PM/MM hybridization intensity differentiation (<1.4-fold) are also passed to step 2.
All putative mutation calls are also passed to step 2 for confirmation. In particular, all high confidence calls resulting in a mutation (different from the corresponding base in the reference sequences used to design the array) were also considered to as a putative type II error. Since mutations may have far-reaching implications in epidemiology studies and drug development against the 2009 Influenza A (H1N1) virus, they were subject to further hybridization intensity analysis in step 2 to confirm the mutation.
Based on empirical observations, 1.4 was set as the minimum fold-change threshold for PM/MM hybridization intensity since ≧99% of the bases called using this threshold are consistent with capillary and 454 generated sequences from the same sample (
This step is used to extract any information out of noisy base calls and to determine the validity of a mutation call.
Determination of Neighbourhood Hybridization Intensity Profile (NHIP) Types
Due to the use of tiling probes in re-sequencing arrays, a single nucleotide mutation at a particular query base could cause a dramatic reduction in the hybridization intensities of neighbouring PM probes up to six bases away. This effect can be measured by studying the NHIP of each query base. The NHIP of each query base is defined as the observed pattern of hybridization intensities of its PM and MM probes and neighbouring (±6 bases from query base) PM and MM probes.
To study the effects of sequence variation (mutation) and noise on the NHIP of a query base, RNA from H1N1 (2009) patient 380 was sequenced by capillary sequencing and on duplicate microarrays. The sequence calls were compared with those generated using Nimblescan or capillary sequencing and a list of true (correct) calls, error calls and ‘N’ (unknown) calls was compiled.
In total, of the expected 13,588 bases of the H1N1 virus (based on genome described at http://www.ncbi.nlm.nih.gov/genomes/taxg.cgi?tax=211044) the microarray according to a preferred embodiment of the present invention called 13,449 bases while capillary sequence was only able to call 12,832 bases. The microarray according to a preferred embodiment of the present invention is thus more reliable, accurate and efficient.
Unlike the NHIPs of wildtype and true-mutation calls, the NHIPs of most errors and ‘N’ calls appear haphazard (
Long chains of consecutive error and ‘N’ calls (especially at the 50- and 30-end of the sample sequences) often have NHIPs where the PM probe of the query base together with neighbouring PM probes, have poor hybridization differentiation with their MM probes (
NHIP analysis showed that all true mutation calls had a characteristic profile (
Nucleotide Substitution Bias Analysis
Re-sequencing arrays rely on the difference in hybridization intensity between a specific hybridization of a PM probe and non-specific hybridization from its MM probes to make a base-call. However, there is evidence that non-specific binding by MM probes depends upon the individual nucleotide substitutions they incorporate. This nucleotide substitution bias implies that a general order in terms of hybridization intensity reduction may exist among the MM probes of each PM probe such that it is possible to compute the likelihood that an observed PM probe is indeed the true PM probe of the sample sequence given the hybridization intensity-based ordering of its MM probes. The key idea is to build a likelihood model of the substitution bias among the probes of non-ambiguous calls on the array; then use this to call bases with ambiguous signals.
The effects of nucleotide substitutions was determined using PM and MM probes (both strands) from high confidence base calls without suspicion of having type I or II errors. There was clear evidence of nucleotide substitution biases shown. The findings from an experiment (305M_A06) is shown in Table 3.
Regardless of strand,
From Table 3, there is strong indication that there exist general orders in terms of hybridization intensity reduction for each PM probe encoding. For example, it is expected that the most frequent hybridization intensity reduction order for PM probes encoding an ‘A’ is TGC since 58% of their MM probes with the substitution ‘T’ suffered the least reduction in hybridization intensity, 50% of their MM probes with the substitution ‘G’ suffered intermediate reduction in hybridization intensity and 65% of their MM probes with the substitution ‘C’ suffered the most reduction in hybridization intensity. There are hybridization intensity reduction orders that are observed primarily for certain PM probes encoding. Thus, if characteristic hybridization intensity reduction orders are identified for each PM probe encoding, then it can be used to ascertain the correctness of a PM probe encoding with some statistical confidence.
Using the same experimental dataset as Table 3, Table 4 shows the enumeration of all possible hybridization intensity reduction orders for each PM probe encoding and their respective frequencies. For each hybridization intensity reduction order, the fraction, fobs, that a hybridization intensity reduction order is observed in the PM probe encoding it belongs to and the random fraction, frand, that the particular hybridization intensity reduction order is seen in other PM probe encodings was computed. Formally, given a PM probe encoding b1 and a hybridization intensity reduction order b2b3b4 where b2, b3, b4≠b1 and b2 has the least reduction in hybridization while b4 has the most reduction in hybridization, then
where t is the total number of hybridization intensity reduction orders excluding b1b2b3b4 obtained from high confidence base calls. Finally, the likelihood that an observed PM probe is indeed the true PM probe of the sample sequence given the hybridization intensity-based ordering of its MM probes is estimated by fobs/frand. Hybridization intensity reduction orders with likelihood scores>2 are statistically significant and are used to discern the PM probe encoding.
For each of the query bases with NHIP of type described in
For query bases that results in a mutation call but have NHIP of type described in
The remaining query bases that have NHIP of type described in
Since nucleotide substitution biases may vary depending on the experimental conditions, experimental reagents or input samples, for each experiment, a set of high-confidence base-calls are obtained and used to infer the hybridization intensity reduction orders for each PM probe encoding. This is then used to compute likelihood “l” scores for base-calling non-high-confidence query bases and mutation confirmation.
The substitution bias on this platform was determined by comparing the PM and MM probes (of both strands) of 25,028 true calls made by PBC from two replicate microarray experiments of patient sample 380. For each true call, a hybridization intensity reduction order was generated by ranking the PM and MM probes of a particular strand in decreasing order of hybridization intensity and recording their respective frequencies (Table 5). Table 5 shows that for each PM probe encoding, certain hybridization intensity reduction orders occur much more frequently than others. For example, if the PM probe encoding is ‘A’ (regardless of strand), then it is most likely that the hybridization intensity reduction order is ‘TGC’ or ‘GTC’. Thus, by matching the hybridization intensity reduction orders of its PM/MM probes with that in Table 5, the likelihood that the putative base call for a query base was determined. In this way, base calls of ambiguous query bases exceeding a reasonably high likelihood threshold and achieve better accuracy and call rate than PBC was recovered.
Another heat map based on the percentage identity of the call sequence to the reference sequence measured at 50 bp windows generated from EvoISTAR is shown in
The map template consists of all eight segments of the 2009 influenza A(H1N1) virus and the locations of known drug binding sites (marked with grey lines) on the NA gene. Locations of all mutation calls are denoted by dark grey triangles beneath the heat map bar. Sequences that are of low coverage (<90%) are automatically flagged, and the overall PM/MM discrimination ratio for each segment is displayed. The heat map bar allows the technician to rapidly assess the quality of the sequence data obtained from the microarray and identify regions where PCR did not work well, or presence of potential recombination/reassortment events. Other details such as coverage, number of base calls successfully made, number of mutations and number of ‘N’ calls for each sequence call are also shown on the visualization map.
Six pairs of replicate experiments consisting of one pair of nasal swab (305 A01, 305_A02) and five pairs of cell culture isolates (305_A03, 305_A04; 305_A05, 305_A06; 305A07, 305_A08; 305_A09, 305_A10; 305_A11, 305_A12), belonging to the same patient sample (305) were employed, to determine the robustness of EvoISTAR sequence calls. Of the experiments, two pairs of replicates (305_nasal and 305_cell_cond1) were amplified under the same optimal experimental conditions while each of the other pairs (305_cell_cond2, 305_cell_cond3, 305_cell_cond4, 305_cell_cond5) were amplified under different sub-optimal experimental conditions (simulating experimental volatility). The results were compared with that of the propriety Probabilistic Base Caller (PBC) algorithm used by Nimblegen. This results are shown in Table 6.
On average, EvoISTAR was successful in calling 99.6% of the 13,449 sites of the 2009 Influenza A(H1N1) virus in the six pairs of replicates. Among the sites EvoISTAR called in each pairs of replicates, >99.9% of sites are called identically. In total, there are 10 mutations (compared to the reference sequences) in the genomic sequences of the 2009 Influenza A (H1N1) virus in patient sample 305 and all of them were correctly called by EvoISTAR in each experiment. The error rate was 6.22e-06 (i.e. 1 error in 1,60,750 bases called) since only one base was wrongly called by EvoISTAR in all 12 replicate experiments. By comparison, PBC was successful in calling only 94.3% of the total possible sites. Although PBC managed to correctly call all 10 mutations present in sample 305, it has a relatively high error rate of 0.006 (i.e. 1 error in 165 bases called). In particular, PBC performed badly on nasal swab replicates 305_A01 and 305A02, achieving only up to 86% coverage and >1.5% error rate. There may have been two likely causes: (1) nasal swab samples have much less concentration of virus RNA than cell cultures, and (2) abundance of human DNA in the nasal swab samples. In comparison, EvoISTAR suffered only a slight drop in performance 98.9% coverage) when analyzing these nasal swab samples.
The comparison was repeated and it was shown that compared with the available capillary sequences for sample 305, EvoISTAR had an average error rate of 0.0012% and 28 ambiguous calls per sample (338 in total). On the other hand, Nimblescan PBC obtained a relatively higher average error rate of 0.169% and 237 ambiguous calls per sample (2855 in total). EvoISTAR is thusrobust and performs well when samples are prepared under sub-optimal conditions. Even for nasal swab samples that tend to have much less concentration of virus RNA than cell cultures, EvoISTAR suffered only a slight drop in performance compared to Nimblescan PBC.
To further validate the software, 14 patient samples were hybridized in duplicate onto the microarray. The microarrays were analysed in parallel using Nimblescan (PBC algorithm) and EvoISTAR, and the sequences obtained were compared to Sanger capillary sequencing. The number of true-non-mutation calls, true-mutation calls, error calls and ambiguous (‘N’) calls were counted for both methods. The substitution bias was also confirmed in all 14 duplicate hybridization experiments (Table 7) to be consistent with that found in Table 5. Compared with the available capillary sequences for the 14 samples, EvoISTAR had an average error rate of 0.0029% and 12 ambiguous calls per sample (346 in total). This is far superior to Nimblescan PBC, where had an average error rate of 0.083% and 158 ambiguous calls per sample (4,434 in total). EvoISTAR also called all true mutations correctly. The genome coverage attained by EvoISTAR (99.02±0.82%) was also much higher than that of Nimblegen PBC (94.3±6.06%).
More than 70% of the 65 error calls (false mutation calls) made by PBC did not have the characteristic NHIP of a true-mutation shown in
To investigate the effects of a re-assortment event on the array, independently amplified segments 1, 2, 3, 5, 6 and 7 of the 2009 influenza A(H1N1) virus and segment 4 of a H3N2 influenza A virus, were hybridized onto an array according to the preferred embodiment of the present invention. The visualization map of this experiment is shown in
The sequence call for segment 4 [based on PM/MM probes from the segment 4 consensus of the 2009 influenza A(H1N1) virus] is poor in quality and coverage. Good base calls from region 1150-1547 was obtained. This region turns out to be the only significantly similar (70% matched) region between the segment 4 (SEQ ID NO:4) consensus of the 2009 influenza A(H1N1) virus and segment 4 of a H3N2 virus (CY039087). This shows that identifying regions of high similarity between the 2009 influenza A(H1N1) virus with other influenza viruses and checking if these regions have good sequence calls may be a plausible way of detecting re-assortments.
| Number | Date | Country | Kind |
|---|---|---|---|
| 200906588-9 | Sep 2009 | SG | national |
| Filing Document | Filing Date | Country | Kind | 371c Date |
|---|---|---|---|---|
| PCT/SG2010/000371 | 9/29/2010 | WO | 00 | 3/29/2012 |
