Patent Application
20130173177
1. Technical Field
This document relates to materials and methods involved in nucleic acid sequence analysis. For example, this document relates to methods and materials for distinguishing sequence analysis errors (e.g., false sequence calls and/or missed true sequence calls) from true polymorphic sequence variations (e.g., single-nucleotide polymorphisms, sequence insertions, sequence deletions, or combinations thereof).
2. Background Information
Knowledge of DNA sequences has become indispensable for basic biological research, other research branches utilizing DNA sequencing, and in numerous applied fields such as diagnostic, biotechnology, forensic biology, and biological systematics. The advent of DNA sequencing has significantly accelerated biological research and discovery. For example, the discovery of disease related regions can aid in diagnosing and treating such diseases.
This document relates to materials and methods involved in nucleic acid sequence analysis. For example, this document relates to methods and materials for distinguishing sequence analysis errors (e.g., false sequence calls and/or missed true sequence calls) from true polymorphic sequence variations (e.g., single-nucleotide polymorphisms, sequence insertions, sequence deletions, or combinations thereof), present in a population. Such methods and materials can be used to provide highly accurate sequence information from large data sets that can provide insight into human evolution, aid in the discovery of disease related regions, and provide knowledge of currently unexplored areas of a genome.
In general, one aspect of this document features a method for assessing nucleic acid sequence information. The method comprises, or consists essentially of, (a) obtaining a collection of at least five sequence output data sets, wherein each of the sequence output data sets comprises a determined sequence that is assembled from a collection of sequence reads of a nucleic acid region and that is aligned to a reference sequence to identify a sequence difference between the determined sequence and the reference sequence, wherein at least one assembly or alignment parameter used to assemble or align the determined sequence is different for each of the sequence output data sets, and (b) determining whether the sequence difference is (i) a processing artifact or (ii) a true sequence difference present in the nucleic acid region as compared to the reference sequence based on a rule set established for the collection of at least five sequence output data sets. The nucleic acid region can be a region of a human chromosome. The collection of sequence reads can be a collection obtained using a second generation sequencing technique. The collection of sequence reads can comprise sequence reads ranging from about 25 to 250 nucleotides in length. The determined sequence for each of the sequence output data sets can be different. The collection of at least five sequence output data sets can be a collection of nine or more sequence output data sets. The at least one assembly or alignment parameter can be selected from the group consisting of a mutation percentage parameter, a coverage parameter, an alignment method parameter, and a matching base parameter. The determined sequence of at least one of the sequence output data sets can be assembled or aligned using a matching base parameter of between 40 and 60 percent. The determined sequence of at least one of the sequence output data sets can be assembled or aligned using a matching base parameter of greater than 90 percent. The determined sequence of at least one of the sequence output data sets can be assembled from a collection of forward paired end sequence reads. The determined sequence of at least one of the sequence output data sets can be assembled from a collection of forward paired end sequence reads and not reverse paired end sequence reads. The determined sequence of at least one of the sequence output data sets can be assembled from a collection of forward paired end sequence reads and reverse paired end sequence reads. The sequence difference can be a single nucleotide difference. The sequence difference can be a single nucleotide deletion. The sequence difference can be a multiple nucleotide deletion or insertion. The sequence difference can be a complex deletion.
In another aspect, this document features a method for assessing a mammal for homozygosity or heterozygosity. The method comprises, or consists essentially of, (a) obtaining a collection of at least five sequence output data sets, wherein each of the sequence output data sets comprises a determined sequence that is assembled from a collection of sequence reads of a nucleic acid region, wherein at least one assembly parameter used to assemble the determined sequence is different for each of the sequence output data sets, and (b) determining whether the mammal is homozygous or heterozygous for a sequence within the nucleic acid region based on a rule set established for the collection of at least five sequence output data sets.
In another aspect, this document features a method for assessing a mammal for homozygosity or heterozygosity. The method comprises, or consists essentially of, (a) obtaining a collection of at least five sequence output data sets, wherein each of the sequence output data sets comprises a determined sequence that is assembled from a collection of sequence reads of a nucleic acid region and that is aligned to a reference sequence of the nucleic acid region, wherein at least one assembly or alignment parameter used to assemble or align the determined sequence is different for each of the sequence output data sets, and (b) determining whether the mammal is homozygous or heterozygous for a sequence within the nucleic acid region based on a rule set established for the collection of at least five sequence output data sets.
Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. Although methods and materials similar or equivalent to those described herein can be used to practice the invention, suitable methods and materials are described below. All publications, patent applications, patents, and other references mentioned herein are incorporated by reference in their entirety. In case of conflict, the present specification, including definitions, will control. In addition, the materials, methods, and examples are illustrative only and not intended to be limiting.
The details of one or more embodiments of the invention are set forth in the accompanying drawings and the description below. Other features, objects, and advantages of the invention will be apparent from the description and drawings, and from the claims.
This document provides materials and methods involved in nucleic acid sequence analysis. Any appropriate sample can be used to obtain a nucleic acid for nucleic acid analysis. For example, nucleic acids can be obtained from blood samples or tissue samples. In some cases, a blood sample, a cheek swab sample, or a hair sample can be used to obtain nucleic acid. Any type of nucleic acid can be used including, without limitation, genomic DNA, cDNA, or plasmid DNA. In some cases, genomic DNA obtained from a human can be used.
Once the nucleic acid sample is obtained, the nucleic acid can be amplified. For example, a portion of a chromosome, a portion of a gene of interest, or a non-coding region within a genome can be amplified. In some instances, introns, exons, 3′ untranslated regions, 5′ untranslated regions, and/or promoter regions can be amplified. Any appropriate method can be used to amplify a region of nucleic acid. For example, long-range PCR or short-range PCR can be used to amplify a region of nucleic acid. In some cases, nucleic acid can be sequenced without performing a nucleic acid amplification process.
Once the nucleic acid is obtained and/or amplified, it can be fragmented into smaller segments. Any appropriate method can be used to fragment nucleic acid. For example, adaptive focused acoustics (e.g., sonication), nebulization, and/or enzymatic digestion with, for example, DNAse I can be used to generate nucleic acid segments. In some cases, restriction enzymes (e.g., BglII, EcoRI, EcoRV, HindIll, etc.) can be used to fragment nucleic acid. In some cases, more than one reaction enzyme (e.g., a combination of two, three, four, five, or more restriction enzymes) can be used. The resulting fragmented nucleic acid can range in length from about 20 to about 1500 base pairs (e.g., about 50 to about 1200 base pairs, about 100 to about 1000 base pairs, about 150 to about 800 base pairs, about 150 to about 500 base pairs, or about 150 to about 300 base pairs). In some cases, the fragmented nucleic acid can be separated based on size. For example, fragments between about 100 and about 300 base pairs (e.g., about 200 base pairs) in length can be separated from larger and smaller fragments using standard fractionation techniques. In some cases, nucleic acid can be sequenced without performing a nucleic acid fragmentation process.
Once the nucleic acid is obtained, amplified, and/or fragmented, it can be sequenced using any appropriate sequencing techniques. For example, adaptors can be added to the nucleic acid which is then subjected to, for example, Illumina®-based sequencing techniques. Such adaptors can provide each fragment to which they are added with a known sequence designed to provide a binding site for a primer that is used during the sequencing process. Other examples of sequencing techniques that can be used include, without limitation, Sanger sequencing, Next Generation Sequencing (or second generation sequencing), high-throughput sequencing, ultrahigh-throughput sequencing, ultra-deep sequencing, massively parallel sequencing, 454-based sequencing (Roche), Genome Analyzer-based sequencing (Illumina/Solexa), and ABI-SOLiD-based sequencing (Applied Biosystems). In some cases, Illumina®-based sequencing techniques are used to sequence a large number of nucleic acid fragments that were generated from long range PCRs. In some cases, nucleic acid from different individuals (e.g., two, three, four, five, six, seven, eight, nine, ten, eleven, twelve, or more different humans) can be sequenced at the same time. In such cases, unique adaptors can be used for each individual such that each sequenced fragment can be assigned to the particular individual from which the fragment originated.
Once the nucleic acids are sequenced, the resulting sequence reads can be assembled and aligned to a reference sequence. Any appropriate sequence can be used as a reference. In some cases, a reference sequence can be obtained from the National Center for Biotechnology Information (e.g., GenBank). Any appropriate software program can be used to assemble and/or align sequences, including, for example, NextGENe® software. In some cases, alignment methods such as BLAT and/or BLAST can be used.
As described herein, the alignment and/or assembly can be performed with stringency and other settings or parameters, such that multiple outputs (e.g., four, five, six, seven, eight, nine, ten, eleven, twelve, 13, 14, 15, 20, 25, or more outputs) are generated. Each output can include a determined sequence that is based on a different set of alignment and/or assembly parameters. For example, a collection of five or more (e.g., six or more, seven or more, eight or more, nine or more, ten or more, eleven or more, twelve or more) output data sets can be obtained with each determined sequence being based on either a highly stringent, a moderately stringent, and a less than moderately stringent set of assembly/alignment parameters.
Each output can include one paired end (e.g., a forward paired end read or a reverse paired end read) in the absence of the other paired end or both paired ends assembled together. For example, a collection of seven output data sets can include (1) a first output data set of forward paired end sequence reads that were aligned and assembled using a first set of parameters, (2) a second output data set of the reverse paired end sequence reads that were aligned and assembled using the same first set of parameters, (3) a third output data set of forward paired end sequence reads that were aligned and assembled using a second set of parameters, (4) a fourth output data set of the reverse paired end sequence reads that were aligned and assembled using the same second set of parameters, (5) a fifth output data set of forward paired end sequence reads that were aligned and assembled using a third set of parameters, (6) a sixth output data set of the reverse paired end sequence reads that were aligned and assembled using the same third set of parameters, and (7) a seventh output data set of both forward and reverse paired end sequence reads that were aligned and assembled using a fourth set of parameters. Any combination of parameters can be used to generate additional output data sets, whether using forward sequence reads, reverse sequence reads, and both forward and reverse sequence reads. In some cases, the order for each output of an output data set can be maintained for each analyzed sample.
Once the output data sets are generated, comparison of each determined sequence to a reference sequence can be performed to identify any sequence differences. These sequence differences can be assessed across each output data set to determine whether the sequence difference is a true difference with respect to the reference sequence, or whether the sequence difference is a false difference (e.g., a false sequence call and/or a missed true sequence call). For each collection of output data sets, a rule set can be established using a known nucleic acid sample having various known sequence differences, e.g., SNPs and indels, as compared to a reference sequence. This established rule set can be used to assess additional sequences to distinguish true sequence difference (e.g., a SNP) from sequence analysis errors (e.g., false sequence calls and/or missed true sequence calls).
In some cases, patterns can be identified that correspond to true sequence differences (e.g., SNPs or indels) as opposed to sequence analysis errors (e.g., false sequence calls and/or missed true sequence calls). For example, when using a collection of nine output data sets, as disclosed herein, the presence of a single nucleotide difference in all output data sets can indicate that that single nucleotide difference is a true SNP. Likewise, the presence of a single nucleotide difference in the first eight outputs (paired ends run separately) and not the ninth output (paired ends run together) can indicate that that single nucleotide difference is a true SNP. In addition, the presence of a single nucleotide difference in only the ninth output, and not in the first eight outputs, can indicate that that single nucleotide difference is a true SNP. Other patterns can be found in Table 1. It is understood that other rule sets can be used for other different collections of output data sets.
As shown in Table 1, pattern D is when experiment 1, paired end 1 called a variant, and experiment 5 called a variant. No other experimental settings called a variant. This pattern was seen for true variant sites nine times in our sample set.
In some cases, when analyzing deletions, the deletions can be grouped into (a) simple or single (non-multi), (b) simple (multi), and (c) complex deletions. A simple (multi) deletion is a deletion of more than one by of the same nucleotide. A single (non-multi) deletion is a deletion of one bp. A complex deletion is a deletion of more than one by of different nucleotides. For simple or single deletions, the collection of output data sets can be analyzed by percent nucleotide content. For complex deletions, the collection of output data sets can be analyzed by the entire sequenced unit.
Table 2 contains an example of a complex deletion. The deletion is 9 bp. The person is heterozygote for this TGAGCCGAG deletion. The percentages in Table 2 range from 41% for the last G to 60% for the “CCG”. Because of these frequency differences, complex deletions are analyzed as an entire unit. The TGAGCCGAG unit travels together in the experimental settings.
Any SNP, insertion, or deletion that does not meet a rule established for that collection of output data sets can be removed as a sequence or PCR artifact, thereby establishing the final sequence for the analyzed nucleic acid.
In some cases, the methods and materials provided herein can not only be used to differentiate true polymorphisms from artifacts, but also can be used to determine homozygosity and heterozygosity at any particular nucleotide position or positions.
In some cases, a reference sequence presented in GenBank® may represent only a haploid consensus sequence. If the reference shows a “G” at a chromosomal position, a mammal (e.g., a human) could carry the homozygote form and inherited G/G from mother/father, but that mammal may be heterozygote and carry G/A from mother/father. The father having the “A”. In this case, the millions of fragmented DNA reads would consist of both father's and mother's alleles. They would consist of G's and A's. If 50 reads aligned to the reference at that position, 25 would be G's and 25 would be A's. This does not happen very often due to the randomness of this process. The percentages can be different. For example, if 50 reads aligned, only 10 may be A's, and 40 may be G's. In such cases, these results still indicated that the mammal is a heterozygote at this position.
DNA samples from 96 Caucasian-Americans were obtained from the Coriell Cell Repository (Camden, N.J.), Human Variation Panel—Caucasian Panel of 100 (Internet site: “www” dot “coriell” dot “org” slash). In addition, ten tumor samples and four anonymized clinical samples were used. Written and informed consent was obtained from all subjects on their use.
The human reference genome was obtained from the National Center for Biotechnology Information, (Build 36 v3; NT—007592.14; subsequence 26,398,617-26,558,272 and NT—016354.19; subsequence 89,146,844-89,218,953). HapMap data for the Centre d'Etude du Polymorphisme Humain (Utah residents with ancestry from northern and western Europe) was downloaded from internet site “http” colon “hapmap” dot “org.” 1000 Genomes Project data was obtained from internet site (“http” colon slash slash “browser” dot “1000genomes” dot “org” slash) and the Single Nucleotide Polymorphism Database (dbSNP) Build 130.
Eleven amplicons, totaling 9.6 kb which targeted 1000 bp of the 5′FR, all exons, and 152 bp of the 3′UTR were produced. Each of the 11 reactions was performed in 20 μL containing 10˜15 ng genomic DNA, five pmol each of forward and reverse primers (Table 3) and FastStart Taq DNA polymerase (Roche). PCR cycling parameters included 95° C. for five minutes, 30 cycles at 95° C. for 30 seconds, 55˜59° C. for 30 seconds, 72° C. for 30˜120 seconds, and a final extension at 72° C. for seven minutes. PCR products were subsequently purified with ExoSAP-IT (USB corp). Amplicons were sequenced on both strands with an ABI 3730 DNA sequencer using ABI BigDye Terminator sequencing chemistry. Four additional regions, totaling 1.1 kb were amplified. All chromatograms were analyzed using Mutation Surveyor v 2.2 (SoftGenetics, LLC, State College, Pa.). The forward and reverse reads were manually inspected.
Long-range PCR (LR-PCR) was performed producing a total of 21 amplicons for each of the 96 Caucasian Coriell samples. The amplicons ranged in size from 3000 bp to 14,581 bp. The 21 LR-PCR reactions used 20-100 ng genomic DNA, 0.4 μM each forward and reverse primers (Table 4), in a total reaction volume of 20-50 μL. The 21 amplicons produced were quantified by the PicoGreen dye binding assay, combined in equimolar amounts and used to create libraries for Illumina GA. The genomic region on 6p was divided into two sections. The first section consisted of nine amplicons and the last section consisted of twelve. An overlap of 19,003 bp resulted from two of the amplicons. The overlap reactions were designated Rxn 20 and Rxn 21.
Paired-end indexed libraries were prepared following the manufacturer's protocol (Illumina). Briefly, 2-5 μg of genomic DNA in 100 μL TE buffer was fragmented using the Covaris E210 sonicator. Double-stranded DNA fragments with blunt or sticky ends were generated with a fragment size mode between 400-500 bp. The overhangs were converted to blunt ends using Klenow and T4 DNA polymerases, after which an “A” base was added to the 3′ ends of double-stranded DNA using Klenow exo− (3′ to 5′ exo minus). Paired-end index DNA adaptors (Illumina) with a single “T” base overhang at the 3′ end were then ligated and the resulting constructs were separated on a 2% agarose gel. DNA fragments of approximately 500 bp were excised from the gel and purified (Qiagen Gel Extraction Kits). The adaptor-modified DNA fragments were enriched by PCR. Indexes were added by 18 cycles of PCR using the Multiplexing Sample Prep Oligo kit (Illumina). The concentration and size distribution of the libraries was determined on an Agilent Bioanalyzer. Four indexed libraries per lane were mixed at equimolar concentrations. Clusters were generated at a concentration of 4.5 μM using the Illumina cluster station and Paired-end cluster kit version 2, following Illumina's protocol. This resulted in cluster densities of 130,000-160,000/tile. The flow cells were sequenced as 51×2 paired-end indexed reads on Illumina's GA and GAIIx using SBS sequencing kit version 3 and SCS version 2.0.1 data collection software. Base-calling was performed using Illumina's Pipeline version 1.0. Reads were converted to FASTA, aligned to the reference and analyzed using NextGENe software v1.04 and v1.10 (SoftGenetics, LLC, State College, Pa.).
An exact test was used to test Hardy-Weinberg equilibrium. Linkage disequilibrium was calculated as the D′ and r2 measures. Tajima's D measures and π (average difference between nucleotide pairs) were estimated as described elsewhere (Tajima, Genetics, 123:585-595 (1989)). Agreement of next-generation sequencing and other genotyping techniques was calculated as the number of sites in agreement between the platforms over total number of sites considered. A confidence interval for this agreement measure was constructed using a sandwich estimator assuming compound symmetric covariance, clusters were individual samples.
JAVA and Perl languages were used in a PolyX program. Excel (2007) VBA and VLOOKUP can also be used for merging the output spreadsheets. This is not in replacement of the automated program. A Perl program parsed the nine NextGENe reports produced by the five experiments for each sample, merged them, and applied “column-based” rules to filter out non-true polymorphic sites. A summary report of the polymorphisms that met the thresholds was produced for each sample. A Java program then collected all of the sample summary reports and applied “population-based” rules to further determine the true polymorphic sites across the population. Input into the rule-set for determining deletions included the so-called “poly-X” program; a Java program that interrogated the reference sequence identifying the length of homopolymers. A structured flowchart (Nassi-Schneiderman diagram) of the overall algorithm, and column and population rules are is set forth in
Five bioinformatic experiments were designed to manipulate two basic settings: reads chosen for alignment, and alignment strategies for chosen reads. All experiments used a median quality score threshold greater and equal to 20, and any reads containing more than three uncalled bases were removed. For the first four experiments, the paired ends were separated and run individually through one cycle of “consolidation.” Consolidation corrects errors in the original reads and elongates them. It also reduces the number of reads by eliminating redundancy. Consequently, the read count and coverage was lower. Although there was less coverage, when using consolidation, it was found to be more uniform across the entire region (
The experimental settings are provided in Table 5.
Because Homopolymers (HPs) are associated with microdeletions and microinsertions, it was attempted to determine how many were within the 160 kb location on chromosome 6p (Denver et al., Abundance, distribution, and mutation rates of homopolymeric nucleotide runs in the genome of Caenorhabditis elegans, Department of Biology, Indiana University, Bloomington, Ind., USA). A HP was defined as being a single nucleotide repeat greater or equal to five by (Ball et al., Human Mutation, 26(3):205-13 (2005)). There were a total of 1,403 HPs within this region, and the lengths ranged from 5-37 bp and decreased in number with increasing length, with only one 37 bp single nucleotide run found. Since the majority of HPs fell within the 5-11 bp range, and PCR and sequencing of these can be difficult, this method was designed to only detect a homopolymer indel, if it fell within a nucleotide run less than or equal to 11 bp. As seen in
A “poly-X program” was written to locate the homopolymers within the genomic region used as a reference and to record their length. This information was integrated into the detection of deletions. Deletions were separated into three categories and paths. Simple (multi) deletions were defined as a greater or equal to two by deletion of the same nucleotide. If it was within a homopolymer region greater than 11 bp, it was ignored. If it was not within the region, the percentages of the nucleotides had to be within 1% of each other, since if they are both deleted, they would be appearing as a unit within the reads most of the time. “Column” and “population” rules were then applied, and the prospective indel put off for manual inspection. Single (non-multi) deletions were defined as a one by deletion. Again, if it was within a homopolymer region greater than 11 bp, it was ignored. If it was not within the region, it was subjected to column and population rules and put off for manual inspection. Complex (multi) deletions were defined as unique, non-repetitive, nucleotide sequences of any size, which consistently appeared as a unit in each experiment. If the frequencies of the nucleotides within this unit were within two percent of each other, it was considered highly reliable. If the frequencies were not within two percent of each other, it was still considered worthy of inspection, as beginnings and ends of reads vary within the alignment, especially if the unit is large (Table 2). Genotypes were determined, the units subjected to column and population rules, and then put off for manual inspection. Insertions had their genotypes determined based on percentages. They were subjected to column and population rules and manually inspected. The actual nucleotide(s) inserted was manually determined.
Because it is preferable to analyze a group of individuals all at one time, it was decided to determine the inherent variability in a real vs. simulated dataset as well as the variability of 96 samples versus one sample. The hypothesis was that the five experiments would establish some consistency and that a set of patterns would emerge among true polymorphic sites. If all settings detected a SNP, a pattern of “9 columns” would result. This would indicate adequate coverage and unambiguous alignments. On the other hand, if only some settings detected a SNP, this would indicate difficulty in alignment or a lack of quality reads, and low coverage. To verify this hypothesis, prior Sanger data was used on 6% of the genomic region under study for all 96 samples. Over the 519 verified markers, patterns emerged. The most frequent pattern for a true polymorphic site was pattern A (
After the column rules were applied to each individual merged datasets, all datasets across the population were combined for each putative polymorphic locus. In a subset verified by Sanger, the total percentage of failed experiments tolerated to maintain reasonable genotype accuracy across the entire population fell within the 0-31% range for SNPs and 0-50% for indels. Higher percentages of failed experiments showed to be inaccurate across 96 subjects, indicative of systematic alignment difficulties within a region which would compromise correct zygosity determinations per sample.
If a position is called to be a polymorphic site, but when looking at the experimental results for all 96 people, most experiments showed no calls, then it can be a difficult area. If it was a good area, all the experimental settings should have mostly picked up a variant at that position. An example of this is position 44049 discussed below. 44049 is a true SNP, but it is in a GC rich area and many experimental settings were not able to detect a variant at that position.
The following population rules were developed:
If a SNP is seen only once in the population, and the percent of failed experiments is greater than 0.25, then remove it.
If a SNP is seen twice in the population, and the percent of failed experiments is greater than 0.25, then remove it.
If a SNP is seen three times in the population, and the percent of failed experiments is greater than 0.30, then remove it.
If a SNP is seen four times or more in the population, and the percent of failed experiments is greater than 0.31, then remove it.
In each case, a failed experiment is one where the parameters selected did not detect a variant at that chromosomal location. For instance, experiment 1-paired end 1 may detect a variant. Experiment 1-paired end 2 may not detect a variant. In this cases, experiment 1-paired end 2 is a “failed experiment.”
Population Rules for Indels (All Indels are Subject to Manual Inspection)
For Simple (multi) and Single (non-multi) deletions, if the percent of failed experiments is greater than 0.50, then remove it.
For Complex deletions, if all members of the unit have a percent of failed experiments greater than 0.50, then remove it. If some members of the unit have a percent of failed experiments greater than 0.50 and other members of the unit have a percent of failed experiments less than 0.50, then do not remove it.
For Insertions, if the percent of failed experiments is greater than 0.50, then remove it.
The following is an example of the implementation of the population rules. The
RefNum is the location of the variant within the subsequence of a contig. It is equivalent to a chromosomal location (Table 6). In this example, a variant was called at position 44049. Hyphens indicate that variant was not called by the experimental setting. For instance, sample CA03 shows that Exp.1 PE2, Exp2 PE1 and PE2, Exp.3 PE 1 and PE2 and Exp.4 PE1 and PE2 did not detect a variant at position 44049.
In this case, the experiments were performed in order. Position 44049 is a true variant site. The rs10947564 is the ID given by dbSNP on NCBI. For the experimental results for CA01 (Caucasian sample 1), the exp. 1-PE1 settings detected a SNP. Exp. 1-PE2 did not. There is a hyphen there to indicate it did not. Exp.2-PE1 also did not detect a SNP. Exp2-PE2 also did not detect a SNP at that position. There are 7 hyphens, which indicate the parameters that did not detect a SNP. Only the settings for Experiment 1-paired end 1 and Experiment 5 were able to detect a SNP at that location.
Twenty-three failed calls out of 36 (four samples X nine possible)=64 percent total failed. This percentage is greater than 0.31, so the putative marker was removed from the final set.
Parameters for genotype calls were developed using NextGENe software and comparison to prior Sanger data. The parameters were as follows:
Parameters for Deletion Variant Calls
Simple (multi): Example CGTTTTACTG (SEQ ID NO: 1) (two by deletion of the same nucleotide).
Homozygous Variant:
A homozygous variant is assigned if any of the five experiments are showing the same nucleotide consecutively less than or equal to ten times, AND that nucleotide equals Ref, and the Ref(s) is within a homopolymer less than or equal to 11 bp OR not within a homopolymer, AND the consecutive Ref nucleotides are within one percent of each other AND Del is greater than or equal to 0.80.
Heterozygote:
A heterozygote is assigned if any of the five experiments are showing the same nucleotide consecutively less than or equal to ten times, AND that nucleotide equals Ref, AND the Ref(s) is within a homopolymer less than or equal to 11 bp OR not within a homopolymer, AND the consecutive Ref nucleotides are within one percent of each other AND Del is less than 0.80.
Single (non-multi): Examples ATCGTCAAT (one by deletion) or ATCGGGGGGTACGC (SEQ ID NO: 2) (one by deletion within a homopolymer less than or equal to 11 bp).
Homozygous Variant:
A homozygous variant is assigned if Ref is within a homopolymer less than or equal to 11 bp OR not within a homopolymer AND Del is greater than or equal to 0.80 AND Ref equals the highest percentage (A, C, G, T).
Heterozygote:
A heterozygote is assigned if Ref is within a homopolymer less than or equal to 11 bp OR not within a homopolymer AND Del is greater than 0.80 AND Ref equals the highest percentage (A, C, G, T).
Homozygous Variant:
A homozygous variant is assigned if any of the five experiments are showing the same consecutive unit (series) of nucleotides AND Del (deletion) percent is greater than or equal to Ref (reference; A, C, G, T) plus 0.40 OR Del percent is greater than or equal to (highest percentage of A, C, G, T, which must equal Ref) plus 0.40.
A homozygous variant is also assigned if some of the nucleotides within the unit show Del percent less than Ref and some show Del percent greater than Ref, then find the member of the unit which has the highest coverage. If the corresponding member of the unit has Del percent greater than the Ref nucleotide, then the entire unit is a homozygote.
Heterozygote:
A heterozygote is assigned if any of the five experiments show the same consecutive unit (series) of nucleotides AND Del percent is less than Ref (A, C, G, T) plus 0.40 OR Del percent less than (highest percentage of A, C, G, T, which must equal Ref) plus 0.40.
A heterozygote is also assigned if some of the nucleotides within the unit show Del percent less than Ref and some show Del percent greater than Ref, then find the member of the unit which has the highest coverage. If the corresponding member of the unit has Del percent is less than Ref nucleotide, then the entire unit is a heterozygote.
Parameters for Insertions
Homozygous Variant:
A homozygous variant is assigned if the Ins percent is greater than or equal to 0.80 AND Ref equals the highest percentage (A,C,G,T).
Heterozygote:
A heterozygote is assigned if the Ins percent is greater than 0.80 AND Ref equals highest percentage (A,C,G,T).
Parameters for SNP Variant Calls
Homozygous Variant:
A homozygous variant is assigned if Alt is greater than or equal to 0.98 and Ref equals 100 minus Alt.
A homozygous variant is also assigned if there are multiple percentages and neither of the two highest percentages equals Ref, then default to the highest percentage variant nucleotide as being homozygous.
A homozygous variant is also assigned if there are multiple percentages and one of the highest percentages equals Ref, and the other highest percentage is greater than or equal to 0.98.
Heterozygote Variant:
A heterozygote variant is assigned if Alt is greater than 0.98, and Ref equals 100 minus Alt.
A heterozygote variant is also assigned if there are multiple percentages and one of the highest two percentages equals Ref, and the other highest percentage is less than 0.98.
Once the genotypes are determined for all five experiments (nine columns), the consensus genotype across all experiments is chosen as the correct one. With this, there is consistency across nine putative duplicate genotypes as a built-in quality control. Replicates can be important. If there is not a clear majority, and the ratio is 50:50, the genotype with the highest coverage is designated as true. In some instances, the reference homozygous genotype is not calculated, and therefore it is not considered in the majority rule to determine the genotype. The reference homozygous genotype is a default genotype to be added at the end of the method.
If a variant is within a region less than the first nucleotide of the forward primer and greater than the last nucleotide of the reverse primer, remove it.
A discordant genotype can be defined if the Next Generation Sequencing (NGS) genotype does not equal the Applied Biosystems, Inc. Sanger genotype.
A discordant genotype can be defined if the NGS genotype does not equal the Illumina genotype.
A discordant genotype can be defined if the NGS genotype does not equal the Affymetrix genotype.
A false variant site is defined as within the boundaries of the PCR forward and reverse primers used for Sanger sequencing if NGS detects either a heterozygote or homozygote variant and Sanger has a homozygous reference. The zygosity, whether true or false, is not considered in this definition. There can be a genotype (zygosity) that is discrepant between the platforms for one or more individuals, but the SNP/Indel marker was still found by NGS since one or more individuals did have the variant.
A missed variant site is defined as within the boundaries of the PCR forward and reverse primers used for Sanger sequencing if NGS did not detect either a heterozygote or homozygote variant among all the individuals and Sanger did detect a heterozygote or homozygous variant. In some instances, the SNP genotype array cannot detect true false variant sites or missed variant sites. It can only determine discordance or concordance. The array can have pre-selected SNPs which are of tested quality and frequency and do not allow for detection of de novo variants (Harismendy et al., Genome Biol., 10(3):R32 (2009)).
A common polymorphism is defined as a DNA variant that is greater than 1% in a population (Roden and Altman, Ann. Intern. Med., 145:749-757 (2006)).
When comparing the five experiments, experiment four, with the different alignment method produced the largest number of called variants with an average (over both paired ends) of 1,113.5 calls. Experiment one resulted in 158.9 calls. Experiment five resulted in 142.5 calls. Experiment two resulted in 128.4 calls. Experiment three, which had the most stringent parameters, resulted in 96.7 calls (
Overall, 613 SNPs and 57 indels were detected (Table 7). Of the 57 indels, 16 were insertions, 41 were deletions, 21 were singletons and 35 had frequencies over 1%. Thirty-four of the indels were within genomic regions of repetitive elements, and 22 were within or immediately next to a homopolymer. The largest complex microdeletion was nine by in length, and the largest structural variant was 3.3 kb in size. Both of these were verified with Sanger sequencing methods. Of the 613 SNPs, 313 were singletons, and 300 were common polymorphisms.
When the 160 kb region was amplified, it was divided into two sections: chr6:35805383-35749634 and 35768636-35648407. This created an overlap of 19,002 bp. Because the 19 kb area was situated between two of the LR-PCR amplicons, there were two different amplifications of the same area. Consequently, the duplicates were compared and used as a built-in quality control. Overall, there were 66 polymorphisms within this region, and most were duplicates. The duplicates mostly showed identical genotypes, and the non-duplicates could be explained by inconsistent coverage in either of the two amplicons. However, there were some exceptions. For example, it was discovered that eight samples had primer SNPs, and preferential amplification of an allele was occurring (Mutter and Boynton, Nucleic Acids Research, 23(8):1411-8 (1995); Walsh et al., PCR Methods Applications, 1(4):241-50 (1992); Quinlanand and Marth, Nature Methods, 4(3):192 (2007)). The LR-PCR conditions were also different between these two amplifications, causing some differences, especially when there was an indel in the vicinity, thus causing additional PCR bias. Where this occurred, all 18 columns plus coverage were taken into account, and the majority genotype with the highest coverage was designated as correct. Because there was not coverage over the entire 160 kb with more than one amplicon, the missed heterozygote rate (MHR) was minimized by building in leniency to the threshold frequencies for alternate and reference alleles in the variant parameters. This was done by excluding a cut-off for the reference.
Several measures were taken to validate results. On a broader scale, heterozygosity (r), was calculated for the entire 160 kb region, as well as regions within the gene structure (Table 8). These values were striking, showing a 40-fold difference (π=0.00002−0.0008) between flanking regions, introns and Untranslated Regions (UTR), intimating unique genetic histories at these loci. Higher heterozygosity in GC-rich areas agreed with previous reports of similar findings (Sachidanandam et al., Nature, 409(6822):928-33 (2001)).
For the entire 160 kb region, π=4×10−4 was obtained (Table 8). This value is much lower than expected, and the negative Tajima's D value of −1.44 conflicts with previous reports of this region on chromosome six as being under balancing selection. Upon inspection, the dissimilar reports were based on small data sets which disregarded low frequency polymorphisms. The complete NGS data shows a dramatic increase in low frequency polymorphisms, thus changing the landscape of evolutionary conclusions.
The Sanger method of sequencing, long considered the “gold standard” for accuracy, was performed. Three thousand three hundred and sixty genotypes were interrogated. No comparison could be made for 53 genotypes because the Sanger results failed. All of these were in intronic areas. Five genotypes were discordant between NGS and Sanger, and these too were in introns. This resulted in a 99.8% concordance between the two methods. With this, it also was possible to determine the number of false variant sites and missed variant sites over the areas amplified. The results showed no false variant sites and no missed variant sites with the exception of “gap 4” (
There was genotyping data on these same samples from the Illumina 550 Kv3 and 510S SNP chips, as well as Affymetrix 6.0. Of the 5,088 genotypes, 17 were failed calls in one or the other technology, and 81 were discordant between the two. This resulted in a 98.4% concordance. Two of the SNPs, one from Illumina (rs7749607) and one from Affymetrix (rs9470065) had been previously genotyped in 96 Coriell Caucasian samples, and they were not found with NGS. To validate this further, Sanger sequencing was used and found the NGS results in agreement for rs9470065 but not for rs7749607. This was not surprising since the single sample in which rs7749607 was found had a reliability index of 3/31, indicating numerous gaps and consequent alignment ambiguities (Table 9). Overall, when combining the two, this method revealed a 98.97% concordance (95% CI: 98.6-99.2). Only one of the Sanger SNPs (rs2143404) overlapped with the Illumina genotyping set. All three methods (NGS, ABI Sanger, and Illumina) revealed 100% concordance across all 96 genotypes.
From the first and last parts of the gene amplified, the gaps greater than or equal to 100 bp were designated “major,” and any gaps less than 100 bp were designated “minor.” Major or minor values were assigned to each of the 96 Caucasians. For instance, NA17259 had values of 16/20-3/2, showing that this person, in experiment five, had 16 coverage gaps of size greater than or equal to 100 bp and 20 smaller gaps less than 100 bp for the first part of the gene. This same person had three major and two minor gaps for the last part of the gene. These values are just indicators of possible trouble and do not represent precise locations.
As a fourth form of validation, dbSNP130 was examined. Two hundred fifty-eight of the SNPs/Indels found were also in dbSNP, although the genotypes across all 96 individuals, utilizing this database, were not available to compare. In several cases, the dbSNP variant, although at the same chromosomal location, did not agree. For instance, at rs35311317 dbSNP has a C/T SNP while NGS found a C insertion. This was validated and the Sanger results agreed with NGS (
The fifth means of validation was assessing whether these genotypes conformed to Hardy-Weinberg equilibrium (HWE) expectations. Deviations from HWE can be due to inbreeding or population stratification, but also can be due to problems with genotyping (Weinberg and Morris, American Journal of Epidemiology, 158(5):401-3 (2003)). Using P>0.001 (Wigginton et al., A Note on Exact Tests of Hardy-Weinberg Equilibrium, Center for Statistical Genetics, Department of Biostatistics, University of Michigan, Ann Arbor Mich.; and Institute of Genetic Medicine, Johns Hopkins University School of Medicine, Baltimore Md.), 25 loci were found to be out of HWE, and none of them were in linkage disequilibrium with each other, indicating the reason they deviated from HWE was most likely due to genotyping errors among one or more samples. Of the loci, 13 were discovered to be within areas of poor coverage and adjacent to large gaps in sequencing. The remaining 12 loci were indels, indicating the zygosity threshold determinations for indels may not be optimal.
The sixth form of validation was testing the method on two additional sample sets. The first set consisted of ten anonymized tumor samples over the same 160 kb region on chromosome 6. One hundred ninety-two kb were verified with Sanger sequencing. All polymorphic sites were detected except where there was no coverage in one of the CpG islands. All genotypes were concordant with Sanger sequencing with the exception of two where the Sanger results failed and therefore a comparison could not be made. The second set consisted of four anonymized and pooled DNA samples over a 5.5 kb region on chromosome 4. All variant sites were detected with no missed sites.
The seventh and final means of validation was the production of a so-called population reliability index based on experiment five (Table 9). The reliability index was to ascertain the number of gaps and therefore explain the remaining discordant calls. Experiment five, unlike the other experiments maintained the original read counts and therefore assured the gaps were not caused by lower coverage because of consolidation of the reads. From the first and last parts of the gene amplified, the gaps greater than or equal to 100 bp were designated “major,” and any gaps less than 100 bp were designated “minor.” Major or minor values were assigned to each of the 96 Caucasians. For instance, NA17259 had values of 16/20-3/2, showing that this subject, in experiment five, had 16 coverage gaps of size greater than or equal to 100 bp and 20 smaller gaps less than 100 bp for the first part of the gene. This same sample had three major or two minor gaps for the last part of the gene. Since these values are just indicators of possible trouble and do not represent precise chromosomal locations, a visual representation was made of the reliability index (
The gap map of
The eight gaps and their genomic contents; short interspersed elements (SINEs), long interspersed elements (LINEs), GC-rich areas, and simple tandem repeats (STRs) are shown in Table 10 below. These repetitive elements were previously reported as causing difficulty with this sequencing technology. Region four has a 77 percent GC content, and region five showed SNPs that were out of HWE.
This accounted for many of the discordant results. However, some results appeared as outliers, not close to any of these areas. Three results, rs6926133, rs9348979 and rs4713904 gave consistent genotype results in all five experiments for three individuals. To verify further, these three individuals were sequenced using Sanger. For two of them, the Sanger results agreed with NGS, but for rs4713904, the Sanger results did not agree. This particular sample had a reliability index of 11/102, thus explaining the discordant genotype for that individual.
In summary, three reasons for discordant genotypes were found: low coverage unique to an individual or common to all samples across a genomic region, other platforms were in error, and preferential amplification.
For the first four in silico experiments, the initial read length of 49 bp increased on average to 66 bp after consolidation, and the percent of alignable reads decreased from 94% to 84%. The correlation between read count and percent of alignable reads were not as expected. For example, NA17222 with a lower read count had 95% alignable reads before consolidation and 91% after. NA17290, with a higher read count had 95% alignable reads before consolidation and 74% after, thus intimating that although original read count is important and a certain minimum threshold is necessary, the quality of those reads, as well as the insert size (Harismendy and Frazer, Biotechniques, 46:229-231 (2009)), is of equivalent importance. The percent alignable reads diminished on average from 68% to 44% after elongation for experiment 5. When comparing the five in silico experiments for numbers of called variants, parameter 4 produced the largest number (1113 calls), and parameter 3 produced the lowest (97 calls). The sensitivity and specificity of individual parameters was also assessed before application of the pattern recognition methodology. While parameter 1 displayed the highest sensitivity (90%)—meaning 90% of the called variants were correctly classified as true, it also showed 65% specificity—an indicator of too many FP. Parameter 3 displayed the highest specificity (86%). Parameter 5 introduced reads with sequencing errors into the alignment, resulting in numerous tri-allelic calls and consequent FP. After application of the pattern recognition methodology, a significant improvement was observed with both specificity and sensitivity increased to 98%.
MAQ (http://maq.sourceforge.net/) is an open source and easy-to-use software that has been used extensively for variation discovery (Clement et al., Bioinformatics, 26:38-45 (2010); Bansal et al., Genome Res., 20:537-545 (2010); Ahn et al., Genome Res., 19:1622-1629 (2009); and The 1000 Genomes Project Consortium, Nature, 467:1061-1073 (2010)). It maps short reads and calls genotypes. MAQ, version 0.7.1 was used to assess 20 of the 96 samples over the 120 kb region on chr6: 35,768,636-35,648,407. Using the default parameters, the SNP filter and loading both paired ends, the SNP and indel calls from MAQ were compared to the results obtained using the pattern recognition methodology. Overall MAQ detected a total of 435 SNPs and 13953 indels in the 20 samples. The pattern recognition methods provided herein identified a total of 292 SNPs and 24 indels. A variant was considered validated if it was seen in Sanger traces, Illumina/Affymetrix data, or dbSNP. From a set of 887 validated sites, the numbers of FP and FN between the two methods were compared. The methods provided herein exhibited 0% FP for both SNPs and indels. MAQ showed 9% FP for SNPs, with only 1.1% of the indels verified as true. As for false negatives, the methods provided herein showed 0.75% and 0.13% for SNPs and indels, respectively. MAQ showed 11% FN for SNPs and 0.26% for indels. To further evaluate the methods provided herein, the SNP and indel calls made using the methods provided herein were compared to those made on the same 20 samples over the same 120 kb region using the SAMtools, version 0.1.16 (Li et al., Bioinformatics, 25:2078-2079 (2009)) and GATK, version 1.1-10 (McKenna et al., Genome Res., 20:1297-1303 (2010)). Using BWA, version 0.5.9 as the aligner and the “mpileup”, “varfilter” and “Unified Genotyper” tools, FP and FN were obtained. The results, using SAMtools, showed 7% FP and 55% FN for SNPs. GATK showed 18% FP and 7% FN for indels. The high FN rate was likely due to this software's very stringent default parameters for calling a SNP or indel.
SNP In Hormone Response Element in LD with Silent (Synonymous) SNP
The method identified a SNP (rs73746499:T>C) at a critical position within a HRE (Hubler et al., Cell Stress Chaperones, 9:243-252 (2004) and Paakinaho et al., Mol. Endocrinol., 24:511-525 (2010)). rs73746499 was found to be at relatively high frequency in the study, with 3.1% of the 96 Caucasian subjects carrying the variant. Further inspection showed 22 additional SNPs and one 5 bp deletion in LD (r2=1) with rs73746499 (
Since the Exon10 and 3′UTR variants were part of the mRNA and both synonymous SNPs and 3′UTR variants have been shown to have functional consequences such as inducing structural changes which could affect protein binding (Nackley et al., Science, 314:1930-1933 (2006); Duan et al., Hum. Mol. Genet., 12:205-216 (2003); Hunt et al., Methods Mol. Biol., 578:23-39 (2009); and Sauna et al., Cancer Res., 67:9609-9612 (2007)), drug interactions or alter mRNA stability, Mfold 3.1 (Zuker, Nucleic Acids Research, 31:3406-3415 (2003)) was used to predict the secondary structures for the full-length wild-type, Exon 10, 3′UTR, and (Exon 10−3′UTR) haplotype mRNA transcripts. The Exon 10 synonymous SNP showed a change in calculated free energy and secondary structure, whereas the wild-type, 3′UTR and (Exon 10−3′UTR) haplotype SNPs showed no changes (
Since RNAs generally adopt multiple conformations, SNPfold (Halvorsen et al., PLoS Genet., 6:e1001074 (2010)) was used to determine whether the SNPs had a large effect on the RNAs structural ensemble. SNPfold computes all the possible suboptimal conformations of the RNA strand and determines the probability of base-pairing for each nucleotide. By evaluating all possible mRNA structures, it was predicted if the SNPs had an affect on the probability of base-pairing (accessibility) of critical interaction sites on the mRNA when compared to the wild-type. According to SNPfold, the Exon 10, 3′UTR, and haplotype (Exon 10−3′UTR) variants significantly disrupted the RNA structural ensemble in specific regions of the mRNA (
Because RNA-binding proteins (RBPs) and ribonucleoprotein complexes (RNPs) partly control gene expression by regulating RNA transcript translation and stability, data obtained by the PAR-CLIP (Photoactivable-Ribonucleoside-Enhanced Crosslinking and Immunoprecipitation; Hafner et al., Cell, 141:129-141 (2010)) method were used to explore whether the FKBP5 mRNA was bound by RBPs and RNPs. Data showed Argonaute (AGO) and trinucleotide repeat-containing (TNRC6) proteins, both part of the miRNA induced silencing complexes (Chen et al., Nat. Struct. Mol. Biol., 16:1160-1166 (2009)), binding to segments of RNA within the 3′UTR of FKBP5. AGO and miR-124, one of the most conserved and abundantly expressed miRNAs in the adult brain (Lagos-Quintana et al., Curr. Biol., 12:735-739 (2002)), were bound to the same site in Exon 9. Insulin-like growth factor 2 mRNA-binding proteins (IGF2BFs) was the most abundant RBP, binding to sites predominantly in the 3′UTR. The methodology provided herein uncovered genetic variants within seven of these binding sites; 5 of which appear to be novel (Tables 11 and 12).
The methods provided herein detected 267 novel rare variants (<1%) within the chromosomal region encompassing FKBP5. The negative Tajima's D value of −1.44 conflicted with previous reports of this region on chromosome 6 as being under balancing selection and upon inspection, the dissimilar reports were based on small datasets which disregarded low frequency variants (Kreitman and Di Rienzo, TRENDS in Genetics, 20:300-304 (2004) and Zan et al., J. Hum. Genet., 51:451-454 (2006)). The complete next generation sequencing data showed a dramatic increase in low frequency polymorphisms, thus changing the landscape of evolutionary conclusions.
Comparison with HapMap and 1000 Genomes Project
Realizing the genetic variation in the CEPH samples may not be identical to that found in these samples, and that the sample sizes are different, it was decided to see if the common polymorphisms detected by this method for this genomic region on chromosome 6 were also present in the HapMap and the 1000 Genomes (1KG) Project data deposited in dbSNP130. All the HapMap CEU common polymorphic sites were in agreement with these findings with the exception of rs3734257, which in the CEU population had a 1.7% frequency in 120 alleles and was monomorphic in our 192 alleles. One hundred sixty-eight common polymorphisms, of which 36% were supported by other platforms such as dbSNP, Sanger, Illumina and Affymetrix, were detected by this method. These polymorphisms were not found in the low or deep coverage 1KG pilots as noted in dbSNP130. Eighty-three of these markers had frequencies greater than 3%. Furthermore, two large gap areas, one of which there was prior Sanger data on, contained high frequency SNPs. These correspond to gaps four and five on the reliability gap map. Gap four is a GC-rich area, and this method was able to detect three out of the three high frequency SNPs within this region: rs9462103, rs13215797 and rs10947564, although because of very low coverage across the entire population and therefore unreliable genotypes, all three were excluded from the final data set. 1KG project detected rs13215797 alone. Gap five contained an Alu, and although a 3.3 kb deletion was detected in some of the samples using Sanger, 1KG also did not detect anything in this area.
These results demonstrate that running more than one experiment reduces the chance of false variant calls in low coverage areas because if one setting does not detect a SNP, another setting may pick it up. This assures that a putative SNP is not disregarded just because it is in a low coverage area, as it would be if only one set of parameters was used.
It is to be understood that while the invention has been described in conjunction with the detailed description thereof, the foregoing description is intended to illustrate and not limit the scope of the invention, which is defined by the scope of the appended claims. Other aspects, advantages, and modifications are within the scope of the following claims.
This application claims priority to U.S. Provisional Application Ser. No. 61/376,641, filed on Aug. 24, 2010. The disclosure of the prior application is considered part of (and is incorporated by reference in) the disclosure of this application.
This invention was made with government support under GM061388 awarded by the National Institute of Health. The government has certain rights in the invention.
| Filing Document | Filing Date | Country | Kind | 371c Date |
|---|---|---|---|---|
| PCT/US2011/048925 | 8/24/2011 | WO | 00 | 3/6/2013 |
| Number | Date | Country | |
|---|---|---|---|
| 61376641 | Aug 2010 | US |
