Method of genotyping by determination of allele copy number

[0001] This application claims priority under 35 U.S.C. 119(e) to U.S. Provisional Patent Application No. 60/280,727, filed Mar. 31, 2001, and U.S. Provisional Patent Application No. 60/313,578, filed Aug. 17, 2001, the entirety of both of which are incorporated by reference herein.

[0003] The invention relates to a method that is used in genotyping as well as to methods of processing data, typically photometric data, which can be used in genotyping and related applications, and more particularly to a method for determining allele copy number using one or more such methods.

[0004] Genetic analysis requires tools for accurately and precisely defining the genetic composition of individuals. Many different genetic technologies have been developed for this purpose, but no single technique is universally applicable to all types of genetic analyses. For example, DNA sequencing technologies recover the precise nucleotide sequence of a targeted (and often very small) region of the entire genome. Numerous studies have shown that this approach readily detects variation among individuals at or above the generic level. However, as typically practiced, nucleotide sequencing generally does not detect variation below the generic level (e.g. populations of a species) due to increased relatedness among individuals. A new class of genetic markers, collectively known as DNA fingerprinting technologies, overcomes many of the current limitations of polymorphism detection inherent to small-scale nucleotide sequencing.

[0005] In its most general sense, DNA fingerprinting, (or DNA profiling), refers to any technology that characterizes an individual's DNA. Most often, this characterization is embodied as some manipulation of the host's chromosomal material followed by a detection procedure, which enables visualization of the test results. Any differences observed between individuals indicate a genetic difference (polymorphism) at a level of resolution dictated by the original characterization procedure, most often individual nucleotides. Generally, the test results of most characterization and detection procedures consist of DNA fragments separated and visualized in an electrophoretic gel. In one best-case scenario, observed DNA fragments will form a pattern unique to each individual sampled. However, depending on the type of genetic characterization employed, not all DNA fingerprints may uniquely identify every individual.

[0006] The ability of any genetic technology to discriminate among individuals depends on many important criteria including: (1) the accuracy and precision of the genetic fingerprinting technology employed, (2) the underlying basis of marker polymorphism, (3) the number of markers generated per reaction (multiplex content) and the number of alleles per marker (information content), and especially (4) the ability to discriminate between homozygote and heterozygote genotypes. Despite the large number of DNA fingerprinting technologies currently available, few, if any, exhibit even a majority of the aforementioned criteria of a genetic marker system. This fact is further exacerbated when considering costs of time and financial investment associated with a particular genetic marker technology.

[0007] For example, in RAPD (Randomly Amplified Polymorphic DNA) markers, random segments of DNA are generated with PCR amplification using short primers of arbitrary nucleotide sequence. Amplification products may be present or absent between individuals reflecting changes in the presence or absence of presumably homologous priming sites. The short lengths of the primers demand that the PCR amplification occur under conditions of low selectivity (stringency). However, low stringency greatly increases the chance of nonspecific priming (primer mismatches) thus generating artifactual DNA fragments ultimately leading to erroneous estimates of polymorphism. Moreover, a great deal of experimental evidence suggests that nonspecific priming events are extremely dependent on individual PCR reaction conditions making reproducibility very difficult. As with most PCR-based fingerprinting technologies, RAPD markers cannot discriminate between homozygote and heterozygote genotypes (cf. dominant markers). This inability significantly reduces the information content of RAPD markers although in certain cases, detailed pedigrees may allow identification of heterozygous individuals.

[0008] Simple Sequence Repeat (SSR) markers are another PCR-based technique where differences in the number of tandemly repeated nucleotide motifs are generated by PCR amplification and detected by gel electrophoresis. Polymorphism in the number of repeats is visualized as length variation across the region containing the repeated nucleotide motifs. Although SSR markers typically exhibit high variation in repeat number thus providing a genetic marker with high information content, SSR markers have a multiplex ratio of one, i.e. only a single SSR marker can be generated per reaction. SSR markers can also discriminate among homozygote and heterozygote genotypes (cf. codominant markers). Despite these advantages, development of SSR markers requires considerable amounts of time. SSR markers developed in one taxon can rarely be used in others. Therefore, SSR markers usually need to be developed de novo for each additional taxa. More importantly, the molecular basis of SSR polymorphism is not completely understood creating under-appreciated difficulties with the analytical and/or statistical treatment of SSR-based data.

[0009] Restriction Fragment-Length Polymorphism (RFLP) markers also detect genetic variation among individuals. RFLP is a DNA hybridization-based technique that reveals polymorphisms as differences in fragment lengths after treatment with a restriction enzyme and electrophoretic separation of the resulting fragments. Two sources of polymorphism contribute to RFLPs: (1) presence and/or absence of restriction sites that determine the number of RFLP fragments generated and (2) length variation caused by insertions or deletions between restriction sites. At first, RFLPs may appear as an attractive fingerprinting methodology especially as RFLP can discriminate among homozygote and heterozygote genotypes. However, RFLP analysis is a very time consuming process, requires relatively large quantities of target and probe DNA, and, as traditionally performed, uses hazardous radioactive labeling of probe DNA. These shortcomings, further magnified with a multiplex content of one and generally low information content, make RFLP analysis a poor DNA fingerprinting choice.

[0010] Alloenzymes, variants of a protein encoded by different alleles at the same locus, have also been widely used to study genetic variation in populations. Although not technically a DNA fingerprinting methodology, protein polymorphisms result from nucleotide base substitutions and thus provide a crude estimate of polymorphism at the DNA level. Nucleotide substitutions in protein-coding loci often lead to the incorporation of different amino acids into the protein. Because of electrical charge variation among individual amino acids, alloenzymes can be differentiated by their relative migration distance during gel electrophoresis. Although alloenzymes are time and cost efficient for research, many protein-coding loci appear invariant within populations (or even at higher taxonomic levels), and most polymorphic enzymes have only a few variants. Recent evidence suggests that this phenomenon may be, in fact, due to experimental artifacts produced by pH gradients affecting the electrical characteristics of the enzyme. In any case, these problems greatly limit the power of alloenzyme analysis to resolve genetic differences among individuals.

[0011] As detailed above, each of these fingerprinting techniques exhibit both strengths and weaknesses. In practice, the choice of a particular technique is often a compromise that depends on the research question pursued and the degree of genetic resolution required, as well as on financial constraints and the technical expertise available.

[0012] A relatively new technique that shows promise as being an ideal DNA fingerprinting technology is Amplified Fragment-Length Polymorphism (AFLP, Vos et al., 1995, Nucleic Acids Research, 23(21) pp. 4407-4414). Briefly, the AFLP technique involves (1) the restriction digestion of total genomic DNA, (2) ligation of adapters, with known sequence, to the digested fragments, and (3) PCR amplification of a subset of these fragments. Products (amplicons) resulting from the AFLP procedure are typically visualized on a denaturing polyacrylamide gel. The resulting amplicons represent a complex DNA fingerprint derived from cleavage sites distributed throughout the entire genome. These amplicons can be isolated and sequenced, but are more often used for establishing unique DNA fingerprints of individuals in a population, without the need for prior sequence information.

[0013] The AFLP procedure exhibits almost all of the desirable qualities in a genetic marker system including: (1) a random distribution of sample sites throughout the genome, (2) a high degree of abundance and polymorphism, (3) a mode of evolution that is primarily selectively neutral, (4) a high multiplex and high information content (due to the large number of loci generated), (4) a well understood basis of marker polymorphism, (5) a capacity for high-throughput amenable to automation, and (6) cost effectiveness and safety.

[0014] However, a serious limitation to AFLP is its inability to directly distinguish whether an observed amplicon was originally derived from one or two doses of template DNA. For example, in a diploid organism, when DNA is processed with the AFLP procedure and an amplicon is generated by PCR amplification, it may be comprised of DNA segments originally amplified from two identical templates, one present on each homologous chromosome pair. Alternatively, the amplicon may be comprised of DNA segments originally amplified from a single template present on only one homologue of the chromosome pair (on which homologue the template resides on is not distinguished). Thus, any detectable amplicon generated by the AFLP process may have been originally derived from either one or two doses of the (assumed) identical DNA template. A third situation also exists in that when neither template is present, no amplicon is produced and is therefore undetectable (null). The above situations are intended to provide examples for the most commonly observed results of the AFLP process and do not reflect all possible template configurations.

[0015] Genetically, the template configurations described above are equivalent to defining the DNA templates as alleles and the position they occupy on a chromosome as a locus. Defined in this manner, the aforementioned configurations simply reflect diploid expectations of the three possible genotypes at a single locus containing two alleles. Therefore, respectively, the three genotypic classes may be identified as having 2, 1 or 0 copies (doses) of the allele (DNA template) distributed at the identical locus (position) on both, one, or neither homologous chromosome pair. The three genotypes of amplicons are herein referred to as homozygote present (having been derived from 2 copies of the allele), heterozygote (having been derived from 1 copy of the allele), and homozygote absent (having neither copy of the allele, null).

[0016] While it is common for multiple alleles to be segregating at a single genetic locus, the AFLP technique cannot generally be used to distinguish among such alleles. This is a direct result of the AFLP technique relying only on template-length variation to discriminate among amplicons. For example, consider six identically-sized amplicons derived from the identical locus in six different individuals. It is possible that one or more of these six amplicons may be comprised of a nucleotide sequence unique to that amplicon. Genetically, each amplicon with a different nucleotide sequence would be considered a unique allele provided it was derived from the same locus. However, since there is no difference in amplicon length, and AFLP-generated amplicons are not routinely sequenced, all six amplicons would be considered derived from identical alleles. In other words, the definition of allele, as is used herein, refers to the original DNA template from which the PCR-generated amplicon was derived. Additionally, allele identification must be restricted to size-homology rather than identity of nucleotide sequence. It is further recognized that in most cases, the exact chromosomal position (i.e. locus) of a given DNA template that results in the production of an amplicon is unknown.

[0017] Therefore, all AFLP-generated amplicons present at a given locus, provided they have identical lengths (as indicated by their molecular weight or molecular mobility) are assumed to be identical-in-state (i.e. allelic). Furthermore, it is assumed that the absence of an amplicon, where others are known to exist, indicates a null allele. As a result, AFLP, as traditionally practiced, can only provide an indication of whether or not a particular amplicon is present or absent (cf. phenotype). Thus, even if the amplicon is present, AFLP, as traditionally practiced, still cannot provide any indication of genotype. For most genetic analyses, the ability to distinguish among each genotype is clearly superior to situations where this is not possible: it allows direct estimation of allele frequencies and for any given level of statistical power, requires smaller sample sizes.

[0018] A recent improvement in the AFLP technique has been the implementation of fluorescent labels and automated DNA sequencers or analyzers to aid in the visualization of AFLP-generated fingerprints. Fluorescent AFLP (fAFLP) is similar to AFLP but implements a modification to the original AFLP procedure where selective amplification is performed with a fluorophore-labeled primer. fAFLP methodologies have greatly increased the throughput for processing samples. However, the same inability to discriminate among homozygote and heterozygote genotypes in AFLP also exists with fAFLP. That is, current techniques available to analyze fAFLP data generated with fluorophore-labeled primers and visualized on automated DNA sequencers still cannot readily determine whether an amplicon was originally derived from one or two copies of an allele. Thus, despite this enhancement to the original AFLP technique, fAFLP still only provides a phenotypic indication of the presense or absence of any amplicon and still cannot provide any indication of genotype.

[0019] While fAFLP remains a powerful genetic technology, it would benefit greatly from a method and system of discriminating between homozygote and heterozygote genotypes, thereby realizing its full potential as a genetic marker system. This capability would dramatically increase the power of the AFLP technology with broad applications in pathotyping, population genetics, quantitative trade loci (QTL) mapping, forensic DNA analysis and many other uses. Such a method is described herein.

[0020] The present invention is directed toward methods of processing data used in genotyping. A preferred implementation of a method of the invention utilizes one or more procedures that enable genotyping to be performed using an emission-based marker system.

[0021] One or more procedures of the method can be implemented to practice some or all of a method of genotyping and the procedures can be used in other biotechnology related data processing applications. Some procedures may not be necessary depending on various factors that include, for example, the type of marker or label system, the type of detection system, the need, convenience, or desire to display specific types of data, as well as other factors. Indeed, it is contemplated that one or more of the procedures of the method can be used or adapted for use in other genotyping methods.

[0022] In a preferred implementation, photometric data generated an apparatus that detects emissions from fluorophore-labeled DNA fragments or other types of labels attached thereto is processed. Preprocessing of recorded data can be performed, such as where photometric data must be reformatted or separated and extracted from other data.

[0023] The photometric data can be baselined to shift all spectral data to a common baseline of zero spectral amplitude. In one preferred baselining procedure, the data is baselined using a regression technique that preferably is multipoint (piecewise) linear regression. For example, the data can be arranged into subsets of datapoints. The datapoint exhibiting the smallest value in each subset is located and the slope between minimum values in adjacent subsets is computed. A baseline offset correction is interpolated for each datapoint using the computed slope. The data is then baselined by subtracting this value from each respective datapoint.

[0024] Where desired, spectral overlap among different fluorophores can also be removed. For example, in a multicomponent system (a system containing more than one detectable fluorophore), the components of fluorescence intensity contributed by fluorophore emissions at wavelengths outside of its principal wavelength or range of wavelengths are removed using a n-fluorophore multicomponenting procedure.

[0025] Noise reduction of the photometric data preferably is performed. In a preferred noise reduction procedure, noise is attenuated and artifacts, if present, removed, reduced or compensated. Noise reduction that preferably is spectral-based noise reduction is performed to reduce or eliminate noise spikes that can be of a random nature. Noise reduction can also be performed to reduce low amplitude background noise, such as what remains after baselining adjustments are performed

[0026] One preferred procedure for reducing spectral-based noise transforms the photometric data from the time domain into a frequency domain using a Fourier Transform that preferably is a Discrete Hartley Transform. To remove high-amplitude spectral noise (e.g. spectral spikes), one or more high-frequency components of the transformed photometric data are truncated at a truncation point that can be empirically determined. To remove low-amplitude spectral noise (e.g. spectral background), partial truncation of the amplitude of one or more low-frequency components can be performed.

[0027] After truncation is completed, the photometric data is transformed back into the time domain using an inverse Fourier Transform that preferably is an inverse Discrete Hartley Transform.

[0028] To compensate for, or attenuate artifacts that can be created by manipulating the frequency components, an apodization function can be applied to the frequency components after truncation but before transformation back into the time domain. In a preferred implementation, a Gaussian apodization function is multiplied to the frequency components of the photometric data that remain after truncation.

[0029] Where it is desired to display the photometric data, the data can be rescaled into a different format that preferably is a more standardized format than produced by commercially available detection systems. It may also be desirable to rescale the data for other purposes. In one preferred implementation, the photometric data is rescaled into an 8-bit hexadecimal data format that can be read by commercial graphic and image-processing programs. If needed, photometric data can be truncated or compressed in order to facilitate resealing and/or storage.

[0030] Where it is desired to display the photometric data where a number of different colors or shades of gray are needed, the data can be mapped accordingly, such as into a false color image. For example, where four different fluorophores are used, the intensities of the four fluorophores are mapped into a three-color space.

[0031] Where lane tracking is needed to define individual samples (e.g. in slab gel electrophoresis), it can be performed. Using graphically displayed photometric data, a tracking spline can be constructed that identifies each sample lane. Such a procedure can be manually or algorithmically performed in software. Where the photometric data is generated from a capillary electrophoresis procedure or the like, lane tracking may not be needed. Where present, lane-tracking information can be extracted from the original set of photometric data, such as from raw data derived from the detector. There may be no need to rescale or display data where lane tracking is not performed manually.

[0032] In a procedure designed to reduce error in determining the relative position of an amplicon within a sample, a function that preferably is an idealized response function, such as a Gaussian function, is fitted to each amplicon in the photometric data. The maximum spectral intensity (peak) of each amplicon is located and fitted with an idealized mathematical function that, if desired, preferably replaces the original data. In locating the spectral peak of an amplicon, a leading edge and/or trailing edge of the amplicon is detected and thereafter the peak apex is located. The apex position is designated when a location, typically a location designated by a scan number, is found where the spectral amplitude of a single amplicon is at a maximum intensity.

[0033] Three parameters are used in fitting the peak with an idealized peak-response function. Such parameters preferably include the position of the peak apex, the amplitude at the peak apex, and the variance (e.g., width) of the amplicon. Preferably, amplicon variance is defined as the width of the amplicon where the amplitude of the peak apex is one-half its maximum amplitude.

[0034] In a preferred implementation of a procedure for fitting the peak with an idealized peak function, estimates of the three peak parameters are inputted into a regression procedure. A preferred regression procedure uses orthogonal distance regression to produce an idealized peak function based on a Gaussian function. Upon convergence, the procedure produces more accurate estimates of peak apex location, peak apex amplitude, and amplicon variance. If desired, the fitted Gaussian function then replaces the data to which it was fit. Preferably, this procedure is repeated for each amplicon located in the photometric data and can be implemented in a manner that fits the Gaussian function simultaneously to all of the amplicons. In order to estimate the relative size of size-unknown sample amplicons, a size-standard is preferably electrophoresed concurrently with the size-unknown amplicons, preferably within each sample lane.

[0035] After fitting both sample amplicons and DNA size standards with the Gaussian function, a sizing curve or sizing function is generated using the fitted parameters for each DNA standard in the set of DNA size standards. The sizing curve or function is thereafter used to estimate the size (e.g. in units of molecular weight or molecular mobility) of each sample amplicon that is of unknown size.

[0036] In a preferred procedure for generating such a sizing curve or sizing function, a locally quadratic, weighted regression approach is used. In a preferred implementation, the estimated peak parameters for a number of standardized DNA fragments less than the total number of standards (a neighborhood) is used in determining sizing functions for each individual standard fragment. The number of DNA fragments used in determining a sizing function for a particular DNA standard is selected so as to minimize the residual error of a least-squares minimization fit to differently-sized subsets of standard fragments. Once the sizing curve is generated, the sizes of size-unknown sample amplicons may be predicted by using the estimated peak apex position of each size-unknown sample amplicon.

[0037] To transform the continuously-valued sizes into discrete representations of identically-sized amplicons, a binning operation preferably is performed using the estimated sizes of all sample amplicons. Any two amplicons are considered homologous (identical-in-state) if the absolute difference in their estimated size is less than or equal to the largest residual of size in all of the size estimates for a given sizing curve. In a preferred implementation, a matrix is created that has an element corresponding to each amplicon and holds a value (e.g., 1 or 0) depending upon whether, at a particular amplicon size-class (cf. locus), an amplicon is either present (1) or absent (0) in each sample analyzed. Each such locus holding uniform values of 1 across all samples is designated as being putatively monomorphic. Any locus holding at least one value of 0 is designated as being polymorphic.

[0038] Using the parameters obtained from the photometric data that correspond to each amplicon, the spectral intensity of the amplicon is estimated using a model that describes the expectations of ssDNA diffusion and displacement in a denaturing electrophoretic gel matrix. The spectral intensities of all amplicons at each putatively designated monomorphic locus are estimated to determine a locus that can be designated as a reference monomorphic locus. This reference locus is used to estimate normalizing coefficients that, when applied to the spectral intensities of other amplicons, will enable direct comparison of relative spectral intensities of amplicons at other loci. Once normalized, the spectral intensities of all other amplicons are able to be genotyped.

[0039] In a preferred implementation of a procedure for estimating the spectral intensity of an amplicon, the estimated parameters of peak apex amplitude and amplicon variance are used. One preferred implementation of a two-dimensional procedure produces an analytical result equal to the square root of four multiplied by multiplied by the estimated peak variance associated with the amplicon multiplied by its peak apex amplitude. The estimated spectral intensity preferably is computed for each amplicon in the sample.

[0040] In another preferred implementation, the spectral intensity of each amplicon is estimated by fitting a thin-plate spline to the spectral data and then determining the area under the surface by numerical integration.

[0041] In selecting a reference monomorphic locus, the variability in spectral intensity of each amplicon of each putatively monomorphic locus is determined. The locus having the least amount of variability in spectral intensity, as determined by a coefficient of variation, amongst all amplicons at that locus, is designated as a reference monomorphic locus. It is assumed that the reference monomorphic locus is truly monomorphic.

[0042] In a preferred procedure for calculating normalizing coefficients for estimated spectral intensities, the spectral intensities of the amplicons at the reference monomorphic locus are used to determine a normalizing coefficient for all amplicons in each individual sample. The computed normalizing coefficients of each sample are applied to the estimated spectral intensity of each amplicon within that sample.

[0043] In one preferred implementation, the average spectral intensity of all of the amplicons at the reference monomorphic locus is determined. The average spectral intensity is then used to compute the normalizing coefficient for each sample. Preferably, the estimated spectral intensity of each amplicon at the reference monomorphic locus is divided by the average spectral intensity to obtain a normalizing coefficient for each individual sample. Thereafter, all of the spectral intensities of amplicons at other loci are normalized by multiplying the spectral intensity of each amplicon by the normalizing coefficient associated with that sample.

[0044] In a preferred implementation of a genotyping procedure, the normalized values of spectral intensity of each amplicon are scaled to unity. Thereafter, the scale values are used to assign a genotype to each amplicon at each locus. Preferably, this procedure clearly separates the scaled values of normalized spectral intensity such that they can be placed into three genotypic classes by visual inspection: homozygote present, heterozygote, and homozygote absent (null). Preferably, a genotype of homozygote absent is assigned to each amplicon exhibiting a normalized spectral intensity of approximately zero. Preferably, a genotype of homozygote present is assigned to each amplicon exhibiting a normalized spectral intensity of about one. A heterozygote genotype is assigned to each amplicon exhibiting an intermediate normalized spectral intensity that preferably is about one-half of the normalized spectral intensity of amplicons having a homozygote present genotype.

[0045] In one preferred implementation of a method of genotyping, each procedure is used with one or more of the other procedures implemented in software or firmware using a computer equipped with a processor, memory and other computer components, such as a keyboard, a mouse, a display, etc., as needed. Such a computer can be a stand-alone device or can be integrated with some other apparatus, such as a detector or some other apparatus. Where lane tracking information is not needed or can be obtained in another manner, it may not be necessary to rescale and/or graphically display the data. If desired, one or more of these procedures can be used in conjunction with other genotyping methods.

[0046] In a preferred implementation of a method of genotyping, one or more of these procedures are used to process photometric data obtained from fragments of DNA prepared with fAFLP. While such photometric data will contain spectral intensity information from detection of detectable wavelengths emitted by markers attached to the fragments, the method and procedures of the method are also well suited for use where the photometric data contains information or data obtained using other detection processes. For example, the method and procedures of the method can be used to process photometric data obtained from detection of light that need not be visible to the naked eye.

[0047] It is an object of the invention to provide an improved genotyping method that can be used to determine allele copy number.

[0048] It is another object of the inveniton to provide an improved genotyping method that can be implemented in software.

[0049] It is still another object of the invention to provide stand-alone procedures that can be implemented in software that can be used alone or in combination with one or more procedures that make up the above-described genotyping method or another genotyping method.

[0050] It is a further object of the invention to provide an improved genotyping method compared to traditionally-practiced AFLP and other methods of polymorphism screening and detection which do not directly allow discrimination between homozygote and heterozygote genotypes.

[0051] It is a still further object of the invention to provide a method of genotyping to that includes one or more procedures that are implemented in software. Depending on the specific polymorphism screening and detection procedures, one or more of the implemented procedures may not be necessary to fully realize the improved genotyping process.

[0052] It is an object of the invention to provide a method of genotyping that increases the certainty of genotyping when compared to other methods of polymorphism screening and detection, which do not directly allow discrimination between homozygote and heterozygote genotypes.

[0053] It is an object of the invention to provide a method of genotyping that baselines photometric data, removes spectral overlap, attenuates spectral noise, identifies spectral peaks in the photometric data as amplicons to be genotyped, estimates a plurality of peak parameters for each amplicon, provides amplicon sizing information, bins amplicons, assigns polymorphic and putatively monomorphic loci, estimates the spectral intensity of the amplicons, assigns a reference monomorphic locus, normalizes amplicon spectral intensity, and assigns genotypes to amplicons based on the relative intensity of the normalized spectral intensities.

[0054] It is another object of the invention to provide a procedure that can be used to more accurately estimate the size (e.g. in units molecular weight or molecular mobility) of a DNA fragment that can be a PCR-generated amplicon.

[0055] It is still another object of the invention to reduce marker assay error, loading volume error, electrophoresis error, and/or detector error using one or more of the aforementioned procedures.

[0056] It is a further object of the invention to estimate the spectral intensity of an amplicon to infer allele dosage (2, 1, or 0 doses), which in turn is a direct indicator of genotype (homozygote present, heterozygote, or homozygote absent null).

[0057] It is an advantage of the invention that processing a plurality of amplicons in a plurality of samples provides more accurate genotyping results.

[0058] It is another advantage of the invention that allele copy number determination using the method as well as one or more select procedures that make up the method or another like method can be performed in an automated manner using a computer equipped with a processor and memory.

[0059] It is another advantage of the invention that genotyping using the method as well as select procedures of the method can be performed in an automated manner using a computer equipped with a processor and memory.

[0060] It is still another advantage of the invention that allele copy number determination using the method as well as one or more procedures that make up the method can be performed in an automated manner using a processor that interfaces with a memory.

[0061] It is a still further advantage of the invention that one or more of the procedures used in the genotyping method can be used separately and independently.

[0062] It is still another advantage of the invention that estimation of the size (e.g., molecular weight or molecular mobility) of each amplicon is more precise.

[0063] It is still another advantage of the invention that genotyping error is reduced.

[0064] It is still another advantage of the invention in that it does not require prior knowledge of pedigree information to assign genotypes.

[0065] It is still another advantage of the invention in that it enables the recovery of codominant genotypes from dominant phenotypic data.

[0066] These and other objects, features and advantages of this invention will be apparent in view of the following detailed description of the best mode, appended claims, and accompanying drawings in which:

[0067] FIG. 1 is flowchart outlining the steps of a preferred fluorescent-based genotyping method in accordance with the invention;

[0068] FIG. 2 is a diagram depicting a generalized fAFLP technique;

[0069] FIG. 3 is a diagram depicting plateau-phase effects on DNA concentration in PCR-amplified amplicons;

[0070] FIG. 4 is a graph showing raw fluorophore emissions of a relatively large set of scans derived from one-dimensional source data;

[0071] FIG. 5 is an enlarged portion of a graph depicting raw fluorophore emissions prior to baseline adjustment;

[0072] FIG. 6 is a graph depicting the fluorophore emissions shown in FIG. 5 after baseline adjustment;

[0073] FIG. 7 is a graph showing fluorophore emissions for a relatively large set of scans after application of a baseline adjustment to the data shown in FIG. 4;

[0074] FIG. 8 is a graph depicting emission spectra of four different fluorophores illustrating the phenomenon of spectral overlap;

[0075] FIG. 9 is a graph depicting fluorophore emissions after performing baseline adjustment and after performing a four-fluorophore multicomponent procedure;

[0076] FIG. 10 is a graph depicting baselined and multicomponented fluorophore emissions of five molecular-weight standards prior to noise reduction;

[0077] FIG. 11 depicts a frequency-domain graph of a Discrete Hartley Transform (DHT) of the time-domain spectral emission data illustrated in FIG. 10;

[0078] FIG. 12 is a graph illustrating a Gaussian apodization function;

[0079] FIG. 13 is a graph depicting truncated and Gaussian apodized data derived from the DHT illustrated in FIG. 11;

[0080] FIG. 14 is a graph illustrating baselined and multicomponented fluorophore emissions of the five molecular-weight standards illustrated in FIG. 10 after noise reduction;

[0081] FIG. 15 is a graph showing an example of baselined and multicomponented fluorophore emissions of a single noise-spike prior to noise reduction;

[0082] FIG. 16 is a graph depicting baselined and multicomponented fluorophore emissions of the noise spike after noise reduction;

[0083] FIG. 17A illustrates an example of raw spectral data taken from two-dimensional source data;

[0084] FIG. 17B illustrates the same data after performing baselining, multicomponenting and noise reduction procedures;

[0085] FIG. 18A is a graph illustrating a single baselined, multicomponented and noise reduced spectral peak derived from one-dimensional source data;

[0086] FIG. 18B is a graph illustrating a portion of the leading edge of the peak shown in FIG. 18A used in determining peak variance;

[0087] FIG. 18C illustrates the peak shown in FIG. 18A along with a Gaussian peak, shown in phantom, fitted thereto using orthogonal distance regression;

[0088] FIG. 19A is a graph depicting five baselined, multicomponented, and noise-reduced spectral peaks derived from one dimensional source data;

[0089] FIG. 19B is a graph displaying the peaks shown in FIG. 19A each fitted with a Gaussian peak shown in phantom;

[0090] FIG. 20 is an exemplary sizing curve generated from 22 DNA standards;

[0091] FIG. 21A is a contour plot illustrating the distribution of spectral intensities within a single amplicon derived from two-dimensional data;

[0092] FIG. 21B is a contour plot showing anomalous distribution of spectral intensities due to amplicon passage through an anisotropic portion of a gel matrix;

[0093] FIG. 22A illustrates a three-dimensional plot of baselined, multicomponented, and noise-reduced spectral intensity data of the amplicon shown in FIG. 21A;

[0094] FIG. 22B illustrates the plot of FIG. 22A after fitting a Gaussian function to it using orthogonal distance regression;

[0095] FIG. 23A is a contour plot of processed spectral data for a single amplicon;

[0096] FIG. 23B is a surface plot of a fitted TPS function to the data in FIG. 23A;

[0097] FIG. 23C is a contour plot of a fitted Gaussian function to the data in FIG. 23A;

[0098] FIG. 23D is a contour plot of the fitted TPS function to the data in FIG. 23A;

[0099] FIG. 24 is a diagram depicting genotypes located for nine autosomal (1-9) and two sex-linked loci (10 and 11) in Domestic Chicken (Gallus gallus) using a preferred embodiment of a method of the invention.

[0100] Before explaining embodiments of the invention in detail, it is to be understood that the invention is not limited in its application to the details of construction and the arrangement of the components set forth in the following description or illustrated in the drawings. The invention is capable of other embodiments or being practiced or carried out in various ways. Also, it is to be understood that the phraseology and terminology employed herein is for the purpose of description and should not be regarded as limiting.

[0101] The following definitions are intended to assist in providing a clear and consistent understanding of the scope and detail of the terms:

[0102] AFLP. Amplified Fragment-Length Polymorphism. A technology used to produce a specific type of DNA fingerprint. In its original form, the AFLP process used radioisotopes to visualize AFLP-generated DNA fragments. See, e.g., U.S. Pat. No. 6,045,994, which is expressly incorporated herein by reference, for a description of the process.

[0103] Allele. A fAFLP/AFLP-generated DNA template derived from corresponding loci on either homologous chromosome and categorized by size (molecular weight or molecular mobility). In one preferred embodiment of the method described herein, the chromosome from which the fAFLP/AFLP allele originates is not distinguished. See locus.

[0104] Amplicon. A quantity of PCR-amplified, identically-sized DNA molecules (fragments) derived from a specific DNA template (allele). The term amplicon is used herein to generally describe any DNA fragment or band produced by the AFLP/fAFLP process and suitably detected.

[0105] Aneuploid. An organism (or cell) in which the chromosome number is not an exact multiple of the haploid number. See genotype.

[0106] Channel. With particular reference to ABI-based detectors, a channel is a one-dimensional vector of data that refers to the number of read regions that are scanned across a gel or capillary (e.g. from left to right) during electrophoresis. In Applied Biosystems (ABI) slab-gel machines, the number of channels varies with options specific to each machine. For example, the number of channels in a 377, 377XL and 377-96 are 194, 388, and 480, respectively. For ABI capillary-based machines, e.g. ABI 3100, there are three channels. Similarly, with other capillary-based machines, the number of channels can also vary. See scan.

[0107] Chromatogram. A generalized reference to a computer file containing photometric data. For example, see gelfile or sample file.

[0108] Codominant marker. Any marker that allows discrimination of homozygote and heterozygote genotypes. See genotype.

[0109] DNA. Deoxyribonucleic acid. An organic polymer composed of four nitrogenous bases (guanine, adenine, thymine, and cytosine) linked via intervening units of phosphate and the pentose sugar deoxyribose.

[0110] DNA fingerprint. Any method that generates profiles of DNA based on the manipulation of the host's chromosomal material.

[0111] DNA fragment. A small segment of DNA, restricted here to those DNA fragments generated by the fAFLP/AFLP process. See allele.

[0112] Dominant marker. Any marker that is scored as either present or absent (null). Dominant markers cannot distinguish between homozygote and heterozygote genotypes. See genotype.

[0113] Electropherogram. A generalized reference to the physical display, in two dimensions, of photometric data, typically by scan number (X dimension) and spectral amplitude (Y dimension). For example, see sample file.

[0114] FAFLP. Fluorescent Amplified Fragment-Length Polymorphism. This procedure is identical to AFLP except that it uses fluorescent labels (e.g., fluorophores) and instruments, such as automated DNA sequencers or densitometers, to detect fAFLP-generated DNA fragments.

[0115] Filter. With particular reference to ABI-based CCD detectors, fluorophore emissions are collected as signal intensities recorded from specific locations on a charged-coupled device (CCD) instrument. These locations are sensitive only to particular wavelengths of light. Thus, the CCD instrument filters predetermined wavelengths of light similar to how physical filters separate individual wavelengths. However, since no physical filtering is actually performed, CCD filtering in this manner is referred to as a virtual filter. The majority of ABI-based CCD instruments possess four virtual filters each corresponding to a specific range of detectable wavelengths.

[0116] Fluorophore. A fluorescent moiety that can be attached covalently to a DNA fragment (typically the 5 position). Commonly used fluorophores include 5-FAM (5-carboxyfluorescein), 6-FAM (6-carboxyfluorescein), HEX (6-carboxy-1,4-dichloro-2,4,5,7-tetra-chlorofluorescein), JOE (6-carboxy-4,5-dichloro-2,7-dimethoxy-fluorescein), NED (proprietary), TAMRA (6-carboxytetramethylrhodamine), TET (6-carboxy-1,4-dichloro-2,7-dichloro-fluorescein), and ROX (6-carboxy-X-rhodamine).

[0117] Fragment binning. A method by which molecular weights or mobilities of fAFLP/AFLP generated amplicons are converted from continuous values (representing size) to discrete values indicating the presence or absence of a given amplicon of a specific size in a specific individual. See sizing curve.

[0118] Gel electrophoresis. The movement of charged molecules through a gel or polymer matrix in an electrical field. It is used herein to define the separation of DNA fragments based on their size (molecular weight or molecular mobility).

[0119] Gelfile. A computer file containing the raw fluorescence data (and possibly other information) recorded from all samples processed on a slab-based automated sequencer. Gelfiles are operationally two-dimensional i.e. one dimension contains the individual samples (e.g. X) which are electrophoresed in the second dimension (e.g. Y). For example, one type of gelfile is an ABI gelfile produced using an ABI 377XL automated DNA sequencer. This gelfile is generated with proprietary ABI Collection Software and organized in an unpublished format. See sample file.

[0120] Genotype. Agenetically distinct class of amplicon inferred from the number of DNA templates (alleles) present at a specific locus from which the amplicon was originally derived. In a model described herein, it is assumed that, but not limited to, at any single locus, there may be at most, two alleles that could be detected and amplified by the AFLP/fAFLP and PCR process, respectively. Operationally, it is assumed (but not required) that alleles from the identical locus are identical-in-state (i.e. contain the identical nucleotide sequence) if they have the same size as determined by, for example, their molecular weight or mobility in an electrophoretic gel. Therefore, for example in diploid systems, only three possible genotypes exist: a homozygote for the presence of the allele (two copies doses of the DNA template are present, one on each homologous chromosome), heterozygote for the presence of the allele (one copy dose of the DNA template is present on which homologue the template resides on is not distinguished) or a homozygote for the absence of the allele (neither copy dose of the DNA template is present on either homologue null). Although the model described herein makes reference only to genotypic configurations present in normal diploid organisms, the model can easily accommodate genotypes in aneuploid or polyploid individuals. In accordance with the invention described herein, the genotype of a locus is thus inferred by the relative amounts of PCR-amplified template DNA present in individual amplicons. That is, the intensity of an amplicon is used to infer allele dosage (2, 1, or 0 doses) which in turn is a direct indicator of genotype (homozygote present, heterozygote, or homozygote absent null).

[0121] With regard to classical genetic terminology, the results of PCR-based assays are actually amplification products (amplicons), observed in a gel and represent a phenotype as the genotype is typically indistinguishable (cf. dominant marker). This distinction is important because confusion with classical genetic terminology incurs definitional ambiguities. For example, amplicons do not have to be derived from a gene per se (and therefore need not be an allele sensu stricto). References to alleles contained herein are made without regard to functionality but rather to homologous chromosomal position (i.e. a locus).

[0122] To understand why a PCR-derived amplicon in a gel is a phenotype, rather than a gene or a genotype, first consider DNA templates that do or do not support PCR amplification. Those templates capable of supporting amplification are defined as 1 and those not capable of supporting amplification are defined as 0. Therefore, in diploid systems, three possible genotypic states exist: 1:1, 1:0, and 0:0. Homozygous 1:1 and heterozygous 1:0 templates both support amplification whereas the 0:0 template does not. Thus, the presence of an amplicon signifies a phenotype associated with both 1:1 and 1:0 templates; the absence of an amplicon represents the phenotype of 0:0. Therefore, the 1 template is dominant to the recessive 0 template because it is present in both homozygous 1:1 and heterozygous 1:0 phenotypes. These precise definitions provide a basis for using classical genetic terminology applied synonymously to amplicon patterns observed in electrophoretic gels:

	Classical allelic	Operational amplicon
	definition.	definition.
Phenotype	1	Dominant	Present
	0	Recessive	Absent (null)
Genotype	1:1	Homozygous dominant	Homozygous present
	1:0	Heterozygous	Heterozygous
	0:0	Homozygous recessive	Homozygous absent (null)

[0123] Heterozygous. Having only one template (allele) at the corresponding loci of either homologous chromosome. See genotype.

[0124] Homozygous absent. Having neither template (allele) at the corresponding loci of both homologous chromosomes. See genotype.

[0125] Homozygous present. Having both templates (alleles) at the corresponding loci of both homologous chromosomes. See genotype.

[0126] Information content. The number of alleles per marker in a sample of individuals.

[0127] Ligation. The process of joining two pieces of DNA (or RNA) with ATP-dependent enzymes called DNA (or RNA) ligases.

[0128] Locus. The classification of a fAFLP/AFLP-generated fragment derived from a specific combination of restriction enzymes and PCR amplification primers, defined by its size (molecular weight or molecular mobility). See allele.

[0129] Marker. See locus.

[0130] Molecular mobility. A physical characterization of molecular size, restricted here to include measurements of DNA fragments with units of length defined by nucleotide bases.

[0131] Molecular weight. A physical characterization of molecular size, restricted here to include measurements of DNA fragments with units of length defined by nucleotide bases.

[0132] Monomorphic. The existence of a single, genetically distinct class of genotype at a specific locus. See genotype.

[0133] Multicomponenting. The process of removing fluorescence overlap between different fluorophores in a multicomponent detection system.

[0134] Multiplex ratio. The number of markers detectable in a single reaction.

[0135] Nucleic acid. A family of molecules that includes RNA and DNA molecules. More specifically, nucleic acid is the phosphate ester polymeric form of ribonucleosides (adenosine, guanosine, uridine or cytidine; RNA molecules) or deoxyribonucleosides (deoxyadenosine, deoxyguanosine, deoxythymidine, or deoxycytidine; DNA molecules), or any phosphoester analogues thereof, such as phosphorothioates and thioesters, in either single stranded form, or a double-stranded helix. Double stranded DNA-DNA, DNA-RNA and RNA-RNA helices are possible.

[0136] PCR. Polymerase chain reaction (PCR) is an enzyme-catalyzed, in vitro copying process of specific DNA sequences that can use extremely small amounts of template DNA. The process includes cycles of denaturation, annealing with primer, and extension with DNA polymerase are used to amplify the number of copies of a target DNA sequence by >10⁶ times. The PCR allows the selective amplification of specific DNA sequences. The PCR process for amplifying nucleic acids is disclosed in U.S. Pat. Nos. 4,683,195, and 4,683,202, which are incorporated herein by reference for a description of the process.

[0137] Phenotype. The observed manifestation of a genotype. The phenotype may be expressed physically, biochemically, or physiologically. See genotype.

[0138] Photometric data. Refers specifically to data comprised related to light that need not be visible. Photometric data typically relates to fluorophore emission intensities where fluorophores are used. In the case of ABI-generated gelfiles, the photometric data is a matrix of spectral intensities from the number of channels across the gel (e.g. 194, 388, or 480) from four filters (0, 1, 2, and 4), multiplied by the number of scans (variable) collected. In the case of ABI-generated sample files, the photometric is organized as described above, but with only a single channel of data.

[0139] Plateau phase. Refers to the phenomenon of non-exponential, and typically linear, production of amplicons during the later cycle numbers of the PCR amplification process.

[0140] Polymerase chain reaction. See PCR.

[0141] Polymorphic. The existence of two or more genetically distinct classes of genotype at a specific locus. See genotype.

[0142] Polyploid. An organism, such as a cell, in which the chromosome number is an exact integer multiple of two or more complete sets of chromosomes. See genotype.

[0143] Primer. A relatively short single-stranded sequence that can bind or anneal to a complementary sequence and can serve as a starting point for DNA synthesis in a PCR reaction.

[0144] Putatively monomorphic. With reference to amplicon binning, the existence of a single phenotypic class (i.e. the presence of an amplicon) at a locus.

[0145] Quantitative trait loci (QTL). Regions of the genome affecting variation in quantitative (phenotypic) traits. QTLs are identified by (1) generating genotypic markers in individuals with a known pedigree, (2) creating a linkage map that shows the order of the markers and relative distance (in centimorgans, cM) between them, and (3) testing for statistical associations between markers (genotype) and phenotypic expression of the trait(s) of interest.

[0146] Restriction enzyme. Enzymes that recognize and cut (restrict) double-stranded DNA at specific short nucleotide sequences. The most commonly used restriction enzymes recognize sequences of four to eight nucleotides in length. Restriction enzymes used in the fAFLP/AFLP protocol produce ends with an overhang of one or more nucleotide bases.

[0147] Restriction site. The nucleotide base-pair sequence of a DNA molecule recognized by a specific restriction endonuclease.

[0148] RNA. Ribonucleic acid. An organic polymer composed of four nitrogenous bases (guanine, adenine, uridine, and cytosine) linked via intervening units of phosphate and the pentose sugar ribose.

[0149] Sample file. A computer file containing processed (and raw) photometric data associated with a single sample. Operationally, sample files are one-dimensional as they either (1) represent a subset of processed photometric data, often derived from slab-based gelfiles, or (2) represent samples derived from individual capillaries in capillary gel electrophoresis. For example, an ABI gelfile is processed with ABI proprietary software (e.g. Genescan Analysis Software) to produce a sample file. The sample file is organized in an unpublished format designed by ABI.

[0150] Scan. With particular reference to ABI-based detectors, a scan refers to the number of traverses across the gel made by the scanning laser. For an ABI-based detector, one complete scan corresponds to four traverses of the laser or exciting device. See traverse.

[0151] Sizing curve. A method to predict the size (e.g. in units of molecular weight or mobility) of size-unknown amplicons that is based on a mathematical relationship between fragment mobility and molecular size of a known DNA standard.

[0152] Spectral overlap. See multicomponenting.

[0153] Tracefile. See sample file.

[0154] Traverse. With particular reference to ABI-based detectors, one traverse is defined as the path the scanning laser makes from channel 0 to the largest channel supported by the specific machine. See scan.

[0155] Traditionally, analysis of fAFLP/AFLP data has been limited to scoring only for the presence or absence, i.e., phenotype, of a particular amplicon because no methodology currently exists to reliably distinguish among confounded homozygote present and heterozygote genotypes in such analyses. In other words, as traditionally-practiced, fAFLP/AFLP cannot produce an intrinsically quantitative result that can distinguish whether an observed amplicon was originally derived from one or two copies of template DNA (alleles). This characteristic results in substantial analytical drawbacks because it reduces total information content and poses peculiar statistical biases not always readily tractable. A method is described herein that can distinguish whether an AFLP-processed and PCR-generated amplicon was originally derived from one or two alleles of template DNA, thereby allowing discrimination among homozygote present and heterozygote genotypes. This method is advantageous as it renders the fAFLP/AFLP technology to serve as a codominant marker system.

[0156] A preferred implementation of a method of this invention is based on the precise quantitation of spectral emissions generated by induced fluorescence of fluorophore-labeled DNA molecules generated in the FAFLP process. Preferably, the spectral intensity of an amplicon is used to infer allelic dosage (2, 1, or 0 doses), which in turn provides a direct indicator of genotype (homozygote present, heterozygote, or homozygote absent null). A preferred implementation of the method assumes that the spectral intensity of PCR-amplified DNA fragments comprising an amplicon is directly proportional to the actual number of DNA molecules contained therein. Since only fluorophore-labeled DNA molecules are detected, the total spectral intensity of all fluorophore-labeled DNA molecules in an amplicon is therefore a good estimate of the actual number of DNA molecules contained within the amplicon. Assuming at most, but not limited to, two identically sized DNA templates (alleles) per locus, heterozygous amplicons are expected to exhibit approximately one-half the spectral intensity as a homozygous amplicon having been originally derived from two copies (doses) of the DNA template. The alternative homozygote produces no spectral intensity as neither copy of the DNA template is present on either homologous chromosome and therefore is not amplified. Thus, the difference in spectral intensity between genotypes is attributable directly to relative differences in the amount of starting template and the number of PCR cycles performed. Overall, the physical amount of DNA present in a fAFLP/AFLP-generated marker, as inferred by the amount of spectral intensity exhibited by the marker, can indicate genotype: heterozygote genotypes should exhibit an intermediate amount of spectral intensity compared to when neither or both DNA templates are present.

[0157] However, due to numerous sources of random variation, e.g., spectral noise or the relative efficiencies of individual PCR amplifications, measurements of spectral intensity differ greatly among amplicons, even at loci with known monomorphic genotypes, and therefore do not assort into the three expected genotypic classes. In slab-gel electrophoretic detection systems, the majority of spectral intensity variation results from: (1) the detection procedure, (2) the electrophoresis conditions, (3) the marker assay, and (4) the sample loading volume. While these same components contribute to variability in capillary-based detection systems, their relative contribution to the overall variability may be different. A preferred approach implemented in a preferred embodiment of a method of the invention to distinguish among genotypes relies on minimizing the variance attributed to these noise components by first applying to the original data a series of mathematical transformations and data-processing procedures. Once applied, the spectral intensity of an amplicon, and therefore the relative amount of DNA contained in the amplicon can be accurately quantified ultimately allowing a distinction between homozygote and heterozygote genotypes ascertained by relative proportions of spectral intensity.

[0158] The major processing steps of a preferred embodiment of a method 40 of the invention described herein, such as is preferably employed for one-dimensional or two-dimensional fluorescent detection systems, is shown in FIG. 1 and include at least some of the following steps: (1) collecting source data 42, (2) creating a spectral baseline using the source data and performing baseline adjustment 44, (3) removing spectral overlap by multicomponent adjustment 46, (4) performing spectral noise reduction 48 by attenuating low- and/or high-intensity frequency components of spectral noise, (5) performing lane tracking 50 to identify individual sample lanes, (6) performing amplicon identification 52, (7) performing amplicon sizing 54 to assign a molecular size to each amplicon or fragment, (8) performing spectral intensity estimation 56 for each amplicon or fragment, (9) performing spectral normalization 58 for all amplicons or fragments, and (10) genotyping 60 each amplicon or fragment. In some instances, the order can differ, one or more steps can be omitted, and, one or more steps not listed above can be performed in addition to one or more the above-recited steps.

[0159] A preferred embodiment of the method was validated in an insect system, namely Drosophila melanogaster, and an avian system, namely Gallus gallus. Each system had known pedigrees so that genotypes could be unambiguously ascertained for both parental and F1 generations. For example, in the avian system, two highly inbred lines of domestic chicken, Ancona hens and a Rhode-Island Red rooster, were crossed. One such mating produced ten healthy F1 individuals. DNA from each individual was processed, six replicates for each parent and two replicates for each F1, using FAFLP and an embodiment of the method described herein. Eleven loci were examined with nine of the loci being autosomal loci and the remaining two loci being sex-linked loci. The correct genotype was assigned correctly for all test loci examined.

[0160] Application of this method to natural populations of birds, namely White-bearded Manakins, Manacus manacus, has also allowed inferences of parentage and dispersal in the subject populations to be drawn that would have been otherwise difficult, if not impossible to achieve. Other applications of a method of the invention include, but are not limited to, systemic pathotyping, population and conservation genetics, disease and quantitative trait locus mapping in humans, plants, or livestock, determination of human identity and genealogical relationships, forensic DNA analysis, and basic population genetic studies on human or non-human biological entities.

[0161] Referring to an example of a fAFLP procedure schematically diagrammed in FIG. 2, a recent improvement in the AFLP process has been the incorporation of fluorescent moieties, commonly referred to as fluorophores, to label individual DNA molecules. This feature allows for convenient detection of labeled DNA fragments 62 and eliminates the need for potentially hazardous radioactive materials, as used in the original AFLP procedure. A fluorophore is attached covalently to the 5 position of one of the two selective PCR primers, typically the primer homologous to the recognition sequence of the less frequently cutting restriction enzyme (e.g. Eco RI). Following PCR amplification, fluorophore-labeled DNA fragments (amplicons) 62 can be separated and viewed on a machine, typically a detection system, such as an automated DNA sequencer, analyzer, or the like. Currently available automated DNA sequencers typically separate DNA fragments based on differences in size by either (1) slab-gel electrophoresis or (2) capillary electrophoresis. DNA fragment separation by size is useful because fragment size generally corresponds to molecular weight or molecular mobility.

[0162] During electrophoresis, amplicons of differing sizes are separated in a matrix that typically is comprised of a polymeric material. In slab-gel electrophoresis, amplicons are separated in a gel matrix comprised of polymerized acrylamide. Operationally, electrophoresis in slab-based gels is two-dimensional as it includes the simultaneous separation of two or more samples in a single gel. The sample being analyzed is applied to the gel across one dimension and is simultaneously separated in the second dimension. For example, as many as 48, 64, or 96, samples can be run simultaneously on ABI slab-gel machines such as the ABI 377, 377XL, or 377-96.

[0163] In capillary electrophoresis, DNA fragments are separated in a narrow-bored capillary containing a liquid polymer matrix. Compared to slab-based gel electrophoresis, capillary electrophoresis is operationally one-dimensional as only one sample at a time passes through the capillary. Multiple samples can be processed simultaneously, but each is done so in its own capillary. For example, ABI 310, 3100 and 3700 capillary-based machines respectively contain 1, 16, and 96 capillaries.

[0164] An excitation device, typically an argon-ion laser, illuminates the gel matrix at a fixed point, called a stage. When fluorophore-labeled DNA fragments in the gel matrix pass by the stage, the fluorophores are excited by the argon-ion laser and their emission spectra are collected by a detector of the detection system that is sensitive to specific emission wavelengths intrinsic to each fluorophore. The detector collects these emissions at time-dependent intervals generally referred to as scans. Electrophoresis continues for a specified amount of time, generally dictated by the maximum fragment size the operator wishes to detect. More than one scan can be performed during analysis of the matrix.

[0165] For clarity, it must also be noted that many automated detection systems can simultaneously detect multiple fluorophores each with emission spectra at distinctly different wavelengths. For example, ABI produces a number of such systems. In contrast to the generalized procedure described above, this allows multiple samples each labeled with a different fluorophore to be electrophoresed simultaneously on a per lane (slab-gel electrophoresis) or per capillary (capillary electrophoresis) basis. Therefore, following the above examples, commonly available ABI machines have the capacity to simultaneously electrophorese 92, 256, or 384 slab-gel samples and 4, 64, or 384 capillary samples.

[0166] Following electrophoresis, fluorophore emission data derived from slab-based gels are stored in computer files commonly called gelfiles or collection files. To meaningfully analyze the data, lane tracking is required to delimit individual samples. Lane tracking is accomplished by defining a single tracking spline, either by a human operator or using a software algorithm, that bisects each amplicon in a single sample lane. Once each sample and its associated amplicons are delimited by a tracking spline, data associated with each sample can be extracted from the gelfile. This procedure is typically implemented in one of two ways: (1) for each time interval or scan, the recorded value of spectral intensity nearest to the tracking spline is extracted, or (2) for each time interval or scan, the average value of recorded spectral intensity in a defined proximity to the tracking spline is computed and stored. This processing step is performed automatically in capillary electrophoresis, believed to be via software control, as emission spectra for each sample are collected independently.

[0167] The end result for conventional processing of photometric data from both slab-gel and capillary electrophoresis is the creation of a sample-specific computer file, commonly called a tracefile, chromatogram, or electropherogram, that, when graphically displayed, depicts the processed signal in two dimensions. Amplicon position is located along one axis, typically the abscissa or X-axis, and relative fluorophore intensity is located along another axis, typically the ordinate or Y-axis. Position corresponds to the scan number, which in turn is proportional to the size of an amplicon.

[0168] The presence at a particular position of a high-density of fluorophore-labeled DNA molecules is represented as an increase in spectral intensity above a defined background threshold. Such a region of increased spectral intensity is typically referred to as the spectral peak of the amplicon. The position of the peak, as a finction of scan number, provides a basis by which to estimate the size of the amplicon (gauged by a standard of known DNA fragment sizes) and ultimately the presence or absence of homologous or unique fragments in other samples.

[0169] Before any source data can be collected, the chromosomal material (total genomic DNA) for each sample must first be prepared for analysis. FIG. 2 depicts an example of a preferred method of preparing DNA for analysis using fAFLP. In the first step, identified in FIG. 2 by the label DNA Preparation, template DNA 64 is obtained by conventional isolation and purification methods. The template DNA 64, e.g., 200 nanograms of purified genomic DNA, is digested to completion with a rarely cutting restriction enzyme, e.g., Eco RI, and a frequently cutting restriction enzyme, e.g., Ase I. Double-stranded adaptors, e.g. 75 picomoles each,

	5-CTCGTAGACTGCGTACCCATCTGACGCATGGTTAA-5	(Eco RI)
	and
	5-GACGATGAGTCCTGAGTACTCAGGACTCAT-5	(Ase I)

[0170] with adaptor ends complementary to the restricted template DNA ends, are ligated with T4 DNA ligase 66 thereby generating a template 68 suitable for PCR, such as is shown in FIG. 2.

[0171] Depending on genome characteristics, the restriction-ligation procedure may generate thousands of adapted fragments. For sufficient resolution of each amplicon, i.e., discrimination of individual PCR-amplified amplicons, to occur in denaturing polyacrylamide gel electrophoresis (dPAGE) or capillary electrophoresis, the population (or complexity) of differing templates must first be reduced. This is accomplished in accordance with the original AFLP process by performing two rounds of PCR amplification of a subset of these templates. Each PCR amplification is designed to reduce the complexity of the template population by preferentially amplifying a specific subset of template DNA.

[0172] In round one, preselective PCR amplification 70 is performed as illustrated in FIG. 2. During preselective amplification, two primers, e.g., 15 picomoles each, 5-GACTGCGTACCAATTCC-3 (Eco RIC) and 5-GATGAGTCCTGAGTAATG-3 (Ase IG) initially reduce the template population complexity by extending the known adaptor sequence of each primer into the unknown portion of the adapted or ligated template by a single nucleotide base, i.e., Eco RIC, Ase IG. Theoretically, each additional nucleotide base extending into the unknown template sequence 70 reduces the complexity of the template population by a factor of 4ⁿ, with n being the number of selective nucleotides on both primers. In the present example, two additional nucleotides are used resulting in the reduction of the template complexity by a factor of 16. The resulting amplicons are diluted, for example, by a factor of 19 with 10 mM Tris-HCl buffer, pH 8.5.

[0173] A second round of PCR amplification, labeled Selective PCR Amplification with Fluorophore-Labeled Primer in FIG. 2, is performed, for example, using two additional selective nucleotides on the Eco RI primer 5-GACTGCGTACCAATTCCAT-3 (Eco RICAT, e.g., 1 picomole), and one additional nucleotide on the Ase I primer 5-GATGAGTCCTGAGTAATGC-3 (Ase IGC, e.g., 25 picomoles). This step further reduces the complexity of the template population by a factor of 64 and overall by a factor of 1024. Preferably, the goal in adding selective nucleotides is to reduce the template complexity to a point where non-homologous amplicons (i.e. amplicons of similar size but of different chromosomal origin) do not comigrate with each other. The choice of specific nucleotides to incorporate into the selective primers is made through routine testing and experimentation.

[0174] Contrary to the preselective PCR amplification procedure 70, the selective PCR amplification procedure 62 uses an Eco RICAT primer labeled with a fluorophore (e.g. 6-FAM) to enable detection of it on an ABI 377XL automated sequencer. Thus, only those amplicons 62 containing an Eco RI restriction site having a fluorophore label will be detected. Those fragments containing only Ase I sites will not be detected as the Ase I primer is not labeled with a fluorophore. In a preferred implementation, selective PCR amplification is performed under conditions of decreasing primer annealing stringency utilizing a procedure known as touch-down PCR amplification. The touch-down procedure begins with the fist few cycles of the PCR amplification procedure under conditions of extremely high annealing stringency, for example at 65 C. Such a high annealing temperature insures that only perfect matches of primer and template will result in template elongation. This greatly reduces the potential of spurious amplifications of imperfect primer/template matches. The annealing temperature is gradually reduced, for example by 1 C. each cycle, until the optimal annealing temperature of the primer is reached, for example, 56 C.

[0175] As typically practiced in the selective amplification phase of the fAFLP/AFLP process, the number of PCR cycles completed often creates conditions where (1) either the PCR primer is prematurely exhausted and/or (2) amplicon production enters a plateau phasewhere amplicon production deviates from exponential expectations and becomes linear. Typically, a plateau phase is entered after completing approximately 30 to 35 PCR cycles. The point at which any amplicon enters a plateau phase is dependent on: (1) the initial concentration of template DNA, (2) the concentration of both PCR primers, and (3) the number of selective PCR cycles.

[0176] Referring to FIG. 3, the plateau-effect phenomenon of PCR amplification may directly prevent discrimination of any relative differences in DNA concentration among amplicons. As a consequence, inferences concerning relative differences in the original amount of template DNA, as reflected in the PCR-amplified amplicon, cannot generally be made. For example, consider any amplicon generated from a 1 and 2 difference in initial template DNA concentration (cf. the difference in expected template concentration when comparing heterozygote and homozygote present genotypes), such as is depicted in FIG. 3. There is a theoretically approximate 7 cycle lag behind which an amplicon derived from the lesser concentration of DNA (e.g. 1) will enter the plateau phase compared to an amplicon derived from a greater concentration of DNA (e.g. 2). If both amplicons have already entered the plateau phase, any method that quantifies the DNA concentration in amplicons, either directly or indirectly, will be unable to infer relative differences in the original template DNA concentration. This is a critical problem as heterozygote genotypes can actually appear as homozygote genotypes and be misidentified as such (FIG. 3).

[0177] To help minimize this problem, the original AFLP protocol is modified in the present invention by (1) increasing the amount of both PCR primers in both amplification steps and by (2) adjusting the number of PCR cycles to restrict entry into the plateau phase preferably maximizing the relative difference in amplicon concentration between homozygote and heterozygote genotypes. In a preferred implementation, amplicons generated in the pre-selective PCR amplification phase, such as by using 15 picomoles each of unlabeled primer, are diluted with 10 mM Tris-HCl, pH 8.5, and adjusted to a concentration of 20 nanograms per microliter. Five microliters of this dilution provides, on average, 100 nanograms template for the selective PCR amplification step. Primer concentrations in the selective PCR amplification step are adjusted to 5 picomoles for the fluorophore-labeled primer and 25 picomoles for the unlabeled primer. Following initial touch-down PCR amplification protocols, the selective PCR amplification phase continues for an additional 17 cycles. Our empirical analyses have indicated that these conditions provide an environment where the vast majority of amplicons will not enter the plateau phase, thus providing a preferred physical criterion by which to accurately and precisely infer differences in initial template DNA concentrations. This ultimately helps to ensure that an accurate assignment of genotype can be made for each amplicon.

[0178] Using the aforementioned template preparation process, different primer pairs, which differ by, for example, a single nucleotide base, will amplify a different subset of adapted DNA templates. By using combinations of primers with different selective nucleotides, a series of fAFLP amplifications will amplify DNA templates (alleles) from a large number of loci in the genome. With the ability to control the number of selectively amplified templates, an optimal number of DNA amplicons may be generated, thereby avoiding complications intrinsic to the fAFLP process itself. For example, excessive DNA template amplifications are avoided. This problem usually produces uninterpretable smears of amplicons when electrophoresed in a gel. In cases where smears are not observed, excessive template amplication often results in unacceptable levels of amplicon co-migration, a phenomenon where amplicons derived from non-homologous loci migrate in a gel at the same molecular weight. Moreover, adjusting the number of PCR cycles in the selective amplification step restricts amplicon production to the log-linear phase, thus preserving a distinction of any differences in initial DNA concentrations.

[0179] In a preferred implementation, amplicons containing fluorophore-labeled DNA fragments generated during selective PCR amplification with the aforementioned modifications are purified over Sephadex columns, preferably using SEPHADEX G-75 marketed by Sigma-Aldrich Chemical Corporation of St. Louis, Mo. This step preferably removes unincorporated primers (both labeled and unlabeled) and nucleotides as well as contaminating salts that may interfere with the electorphoresis procedure. The amplicons are then separated on a dPAGE automated sequencer using a slab gel made with a 5% acrylamide mix, preferably using a 5% Long-Ranger acrylamide mix marketed by BioWhittaker Molecular Applications, Inc. of 191 Thomaston Street, Rockland, Md. A standard containing known sizes of DNA fragments, such as ILS-600, marketed by Promega Corporation of 2800 Woods Hollow Road, Madison, Wis., is included in each sample to control for any lane-to-lane amplicon migration variability and thus provide a basis from which to determine the sizes of the size-unknown amplicons. The DNA size-standard is labeled with a different fluorophore, e.g. ROX, than the fAFLP-generated sample amplicons and allows both to be easily distinguished.

[0180] In a preferred implementation, amplicons are separated and visualized on an automated DNA sequencer manufactured by ABI. A currently preferred model of a slab-based automated DNA sequencer is the ABI 377XL. The ABI 377XL separates and records the presence of fluorophore-labeled DNA fragments at a fixed point in a polyacrylamide gel using a charged coupled device (CCD) detection apparatus. Other instruments may also be used including those that use a complementary metal oxide device (CMOS) to detect spectral emissions.

[0181] In one preferred implementation using an ABI 377XL automated sequencer, the upper buffer chamber is loaded with 0.70TBE buffer and 1.0TBE buffer is loaded in the lower chamber. Heat denatured samples are loaded onto the gel and electrophoresed at reduced power, typically 10 Watts, for two minutes. A solution of 10TBE is then added to the upper buffer chamber to achieve a 1TBE concentration followed by electrophoresis at normal power, which typically is about 40 Watts.

[0182] This procedure, known as a water load enhances fragment stacking, which is a phenomenon where differently sized DNA fragments orient themselves according to size upon initial entry into the polyacrylamide gel matrix. This has a marked effect of reducing the apparent broadness of a group of identically sized DNA fragments (an amplicon) as they migrate through the gel, thereby increasing overall amplicon resolution. The electrophoresis power is set preferably to a limiting 40 Watts, which appears to optimize amplicon cohesiveness by minimizing excessive diffusion (observed at very low power) and by minimizing physical distortions due to sieving effects (observed at very high power). For a 36 cm gel plate, a seven-hour runtime detects amplicon sizes of approximately 20 to 1200 bases.

[0183] Other electrophoresis conditions may also be used. For example, electrophoresis power can be adjusted such that it differs from that discussed above. Samples can be loaded with 1.0TBE buffer in both chambers with or without a two-minute electrophoresis pre-run. Other buffers differing in composition and concentration may alter the electrical characteristics of the system and may beneficially affect the electrophoretic behavior of the DNA fragments. Routine testing and experimentation can be used to determine whether one or more of the aforementioned parameters affect optimal base stacking and amplicon migration.

[0184] Separated fluorophore-labeled DNA fragments (amplicons) are visualized by exciting the fluorophore with an external emission source and recording the intensity of fluorescence emitted therefrom. These data are collected and recorded in real time during the length of the electrophoresis session. Typically, the fluorophore emission data is stored in a computer file as detailed below. For example, ABI Collection Software v2.5, used on an ABI 377XL automated sequencer, creates a gelfile that contains a complete record of the raw fluorophore emissions for an entire electrophoresis session.

[0185] If desired, fluorophore-labeled DNA fragments (amplicons) can be separated and visualized on different types of imaging devices. Such devices include, but are not limited to: the NEN Global IR² DNA Sequencer System from LI-COR, Inc. of Lincoln, Nebr., which detects near-infrared fluorescent dyes, the ALFexpress II DNA Analysis System from Amersham Pharmacia Biotech, of Piscataway, N.J., and the FluorImager, which is a multi-purpose fluorescent sample imager for multiple-color analysis from Molecular Dynamics, a company of Amersham Biosciences of Piscataway, N.J. Fluorophore-labeled DNA fragments can also be separated and detected using instruments that employ capillary electrophoresis technology, such as is used in the ABI 310, ABI 3100 and ABI 3700 Genetic Analyzers. Each of these devices has their own unique operating procedures. Adjustments to the aforementioned slab-gel or capillary-based electrophoretic procedures may be necessary to accommodate the requirements of each detection technology.

[0186] To accurately and precisely distinguish among the three possible genotypic classes of amplicons generated by the fAFLP process, the resultant source data containing the actual fluorophore emission data, henceforth termed photometric data, must first be processed in accord with the invention detailed below. In its original form, the photometric data contains various sources of variability and error, typically in the form of noise, which negatively affects accurate quantitation of spectral intensity, and therefore the amount of DNA in any amplicon. In slab- and capillary-based detection systems, the majority of random variation observed in measurements of DNA concentration is believed to result from differences in: (1) the marker assay, (2) the sample loading volume, (3) the electrophoresis conditions, and (4) the detection procedure. These four sources of variability are set forth in Equation I below,

[0187] where

[0188] ²_Total is the total estimated variability in the source data;

[0189] ²_{Marker Assay} is the total estimated variability due to variations in the marker assay;

[0190] ²_{Loading Volume} is the total estimated variability due to variations in sample loading volume;

[0191] ²_{Eectrophoresis} is the total estimated variability due to variations in electrophoresis conditions;

[0192] ²_Detector is the total estimated variability due to variations in detector operation.

[0193] A present method of distinguishing among genotypes in amplicons generated from fAFLP assays in accordance with a method of the invention relies on first minimizing the variance attributed to these components by applying to the original data a series of mathematical transformations and image-processing procedures. Once applied, the relative amount of DNA in any amplicon, as reflected by spectral intensity, can be accurately estimated and normalized across multiple samples. This procedure enables an inference to the original number of DNA templates (alleles) originally present at any locus thus allowing a determination of the genotype of the amplicon.

[0194] Briefly, the procedure consists of applying baseline, multicomponent, and Fourier-based noise reduction procedures to reduce variation attributed to the detection and electrophoresis apparatus. Sample amplicon sizes are estimated using a statistical model designed to enhance the predictive accuracy and precision of sizing curves generated from size-standardized DNA fragments. Estimation of spectral intensity of each amplicon is based on a statistical model derived from characteristics of single-stranded DNA in denaturing electrophoretic gels. This model provides an accurate means by which to measure the amount of DNA in each amplicon. Normalization of these estimated spectral intensities is accomplished by using a monomorphic locus present in the population of samples. This type of normalization provides a basis for minimizing variation attributed to the marker assay itself and unequal sample loading volumes.

[0195] A purpose of this invention is to apply one or more of the procedures disclosed herein in a manner such that direct comparisons of spectral intensity, and therefore DNA abundance, can be made among amplicons present at any locus. In doing so, a distinction as to whether the amplicon was originally derived from one or two alleles can be inferred. With this information, a distinction between homozygous and heterozygous genotypes can be made for any amplicon. The method described herein is applicable to a broad class of genetic markers and marker-detection technologies and can be achieved using source data comprising a one- or two-dimensional vector of spectral intensities sampled at regular time intervals.

[0196] Although no reasons are known to exist why the method of the invention would not work with other marker assays, detection systems, aneuploid or polyploid organisms, the description set forth below has been limited to DNA-based, fAFLP-generated amplicon assays in diploid organisms visualized on an ABI 377XL slab-gel automated sequencer for the sake of describing at least one preferred implementation of the method. Where practical, it is identified herein what procedures are applicable to specific technologies and what algorithm is used to perform the particular procedure. The order in which each procedure is discussed herein reflects a preferred order of its application.

[0197] When the source data contains information in addition to that required for genotype discrimination, it may be necessary to first preprocess the source data using an algorithm preferably implemented in software that extracts the pertinent data. For example, when using ABI automated sequencers, ABI Collection Software, e.g., version 2.5, is used to record raw fluorophore emissions from the detection apparatus. The data comprising the raw fluorophore emissions is stored in a file that has a proprietary format either as: (1) a gelfile (or collection file), such as for slab-based machines or (2) a sample file (or tracefile), such as for capillary-based machines and for at least certain types of source data derived from some slab-based machines. Data from these files, or files such as these, contain the actual source data. If desired, the source data can also constitute data outputted directly from any such machine or other suitable device in real time and need not be pre-stored in a file.

[0198] In one method of preprocessing, source data generated from an automated sequencer or other like device, source data from a gelfile or sample file are used as input that is processed to locate, extract, and, if desired, manipulate at least some of the data contained therein. If desired, two-dimensional source data, such as photometric data from a gelfile, can be processed to obtain one-dimensional source data that resembles or is similar to data contained in a sample file. As previously mentioned, such source data can come from a slab-based DNA sequencer such as an ABI automated sequencer. ABI-derived data or like data from other types of sequencers or fluorescence imaging devices that produce or contain photometric data can also be used as source data. Hence, the term gelfile is used herein to generically indicate two-dimensional source data and sample file is used herein to generically indicate one-dimensional source data.

[0199] Source data are typically stored in a computer readable file. Photometric data are processed by first locating it within the source data, extracting it, such as by copying it, and, if desired, subsequently formatting it in a manner that preferably simplifies the manner in which it is stored. The extracted and/or processed photometric data can be written to another output file for subsequent use or can be used without being first written to an output file. Since not all of the previously stored data are extracted during processing, the resultant output file is often smaller than the original file containing the source data.

[0200] Depending on how the source data is organized, additional information other than the photometric data itself may be present. This information can be used in locating and identifying photometric data as well as other relevant photometric-related data. Once located, the photometric data are extracted along with any other related data desired or deemed necessary. The data can then be reorganized to a format that is preferably more easily manipulated and inspected or used as needed.

[0201] In one preferred implementation, photometric data stored in an ABI-generated gelfile or sample file, in a form representing the intensity of recorded fluorophore emission data as integer values between 0 and 8191, is extracted. This data may be rescaled into integer values ranging from 0 to 255. Although rescaled, this new range of integer values still represents the relative intensity of fluorophore emission data originally collected. An advantageous result of this procedure is that the rescaled data are easier to manipulate such that it can also be used by other algorithms, such as those found in graphics software, graphics related software, statistical analysis software, gene fragment analysis software, data visualization software, as well as in other types of software. For example, processed and resealed photometric data can be analyzed with freely available gene fragment analysis software, e.g., Cross Checker, which is available from the Department of Plant Science, University in Wageningen, the Netherlands.

[0202] In one preferred implementation where a PC or PC-compatible computer is used, it may be necessary to convert the source data into another byte order, e.g., LITTLE_ENDIAN, before performing any preprocessing procedures. It may also be necessary to convert the source data into another byte order, e.g., BIG_ENDIAN, when UNIX operating systems or Macintosh computers are used to preprocess gelfiles or sample files originally stored on a PC or PC-compatible computer. Software modules are available that easily convert data to and from LITTLE_ENDIAN and BIG_ENDIAN format. One such module is available from the Medical Research Council as part of the Staden Package at the MRC Laboratory of Molecular Biology, Hills Road, Cambridge, United Kingdom. Such a conversion may be needed when going from one type of platform to another. However, such conversion or reformatting may not always be needed.

[0203] The presently preferred source data preprocessing procedure locates, obtains and extracts the photometric information contained within the source data, e.g., either a gelfile or sample file. Header information included in the source data being preprocessed is examined to determine where this information resides.

[0204] Where the source data is an ABI gelfile, it is assumed that the gelfile begins with one or more specific file identifying characters. For example, ABI gelfiles derived from ABI Collection Software v2.5 begin with 11 concatenated header records, which vary in length depending on the number of channels scanned during an electrophoresis run. For example, depending on options specific to individual machines, the ABI 377 series has the capability to scan 194, 388 or 480 channels. Record lengths (LEN_RECORDS) for these scan options are known to be 407, 795 and 979 bytes respectively.

[0205] The initial eleven individual header records are designated by two single ASCII characters beginning at byte (LEN_RECORDS*i) and byte (LEN_RECORDS*i)11, where iheader record 0, 1, 2, . . . 10. Immediately following the ASCII character at header byte (LEN RECORDS*i) is a numerical value designating LEN_RECORDS after subtraction of a portion of the record header and checksum bytes (see below).

[0206] Any information stored within these initial eleven headers appears to begin at byte (LEN_RECORDS*i)17. At present, only 5 of the 11 header records have a known function, including header: (1) ABI collection software version, (3) ABI machine type, (5) Number of referenced data blocks (NUM_RECORDS), (6) offset to the primary photometric data block, and (7) length of the primary photometric data block.

[0207] In contrast to the organization of the initial eleven record headers, all data records associated with actual photometric data begin with the ASCII character > and are followed by a numerical value designating the amount of data, in bytes, stored in each data record after subtracting a portion of the record header, which currently is four bytes in length, and a data checksum, which currently is two bytes in length. The three possible initiators for 194, 388 and 480 channel ABI 377 data records are >401, >789, and >973 respectively. Recorded fluorophore emission data is stored in the primary photometric data block as binary data representing two-byte integers with corresponding decimal values ranging from 0 to 8191.

[0208] Fluorophore emission data is recorded in a channel by filter by scan fashion. Each channel can be considered as a column vector of data values, each representing divisions across the width of the gel or capillary, with scans represented as rows across the column vectors. Each filter comprises a distinct range of detectable wavelengths emitted by a particular set of fluorophores. Most ABI automated sequencers can concurrently detect 4 fluorophores. Thus, the photometric data is a matrix of two-byte integers representing fluorescence intensity at a particular channel by row position, sorted by filter number.

[0209] The first photometric data record is located at a static offset from byte 0 in 194, 388, and 480 channel gelfiles and is 4477, 8745, or 10769 bytes, respectively. Also, depending on the number of channels, the photometric data contains length-records (LEN_RECORD) of 407, 795, or 979 bytes that include a 17-byte sub-record and a two-byte checksum. Within each of these records, the portion of the record allocated to actual photometric data is 388, 776, or 960 bytes, respectively.

[0210] Each data record contained in the primary photometric data block contains 10 fields: (1) RECORD_START, (2) DATA_BYTES, (3) an ASCII space character, (4) GEL_VOLTAGE, (5) GEL_CURRENT, (6) GEL_TEMPERATURE, (7) FILTER_NUMBER, (8) TRAVERSE_COUNTER, (9) SCAN DATA, and (10) SCAN_CHECKSUM.

[0211] The RECORD_START field contains the initiator for a new scan line (>) followed by a three-byte field, DATA_BYTES, which indicates the total record bytes remaining after subtracting the number of data bytes associated with the RECORD_START, DATA_BYTES and SCAN_CHECKSUM fields. The two-byte fields of GEL_VOLTAGE, GEL_CURRENT, and GEL_TEMPERATURE provide values of gel voltage in units of volts/10, gel current in units of milliamperes, and gel temperature in units of degrees Celsius, respectively, at the initiation of each new scan.

[0212] The format in which ABI stores fields 4, 5, and 6 can be confusing. Recorded values are apparently stored as integer-coded two-byte values. However, it is the corresponding hexadecimal representations of these bytes that are actually interpreted. For example, GEL_VOLTAGE is stored as an integer-coded, two-byte value. An example of such a value representative of GEL_VOLTAGE is 61 38. This corresponds to the ASCII values a and 8. These two ASCII values are concatenated to form a hexadecimal value of 0xa8. The decimal equivalent of 0xa8168. As this value is in the units of Volts/10, the actual voltage in this example is 1680 Volts. Similarly, for example, GEL_CURRENT: 31 61ASCII 1a0x1a26 mAmps; for GEL_TEMPERATURE: 33 31ASCII 310x3149 C.

[0213] The field FILTER_NUMBER refers to the filter numbers assigned to detect emissions from particular fluorophores. Where four fluorophores are used, filter numbers 0, 1, 2 and 3 are used. The TRAVERSE_COUNTER field records the cumulative traverses of the laser across the gel. One scan corresponds to 4 traverses of the laser as one traverse is performed for each filter number. The TRAVERSE_COUNTER initializes at 1 and resets to 1 after 2{circumflex over ()}16 traverses of the laser. Immediately following is the SCAN_DATA field, which actually contains the photometric data stored as two-byte integers. Lastly, a two-byte field, SCAN_CHECKSUM, appears to be some sort of checksum resulting from modulus arithmetic on the SCAN_DATA. The form of the checksum is presently unknown.

[0214] In a preferred implementation of a procedure for extracting photometric data, two-dimensional source data, typically in the form of a two-dimensional matrix of raw digital data, is processed. Such is typically the case where slab-gel electrophoresis is used and a gelfile is created. Photometric data in the source input is processed by locating it and extracting it, preferably to another output file for subsequent use.

[0215] In a preferred implementation, an extraction algorithm processes any gelfile derived from an ABI 373 automated sequencer, including units having 194 and 388 channels, and an ABI 377 automated sequencer, including units having 194, 388, and 480 channels. First, the expected eleven header records are checked for the correct content and offsets. Depending on the particular model of automated sequencer, if the initial header does not exhibit the aforementioned features associated with its respective machine model, an error message preferably is generated. Next, the offset to the photometric data and the length of the photometric data is retrieved from the header records 6 and 7 described above. The photometric data is extracted and saved to a new file. For example, in one preferred implementation, a gelfile of length 76,901,936 bytes, derived from an ABI 377XL 388 channel automated sequencer, was processed. Beginning at an offset of 8745 bytes from byte 0 and extracting 16011 scans* 4 data rows per scan* 766 bytes of photometric data per row resulted in a new file, containing only photometric data, of length 49,698,144 bytes.

[0216] The structure of the new file containing the extracted photometric data is greatly simplified compared to the original format. In one preferred implementation, a single line of two, five digit numerals, separated by a space, indicate the number of channels and the number of scans. This information is immediately followed by the extracted photometric data written as rows of two-byte, binary-encoded integer data arranged in a scan by filter by row fashion. For example, following the above criteria, this organization would be represented as:

	00388 16011
	scan 1:filter 0:row 1
	scan 1:filter 1:row 2
	scan 1:filter 2:row 3
	scan 1:filter 3:row 4
	scan 2:filter 0:row 5
	. . .
	scan (n):filter 0:row (4n-3)
	scan (n):fllter 1:row (4n-2)
	scan (n):filter 2:row (4n-1)
	scan (n):filter 3:row (4n-0)

[0217] Note that in this example, the newly created file of extracted photometric data will be 12 bytes larger than expected from the previous formula (i.e. 49,698,156 bytes). This is due to the addition of the 12-byte channel number by scan number header information.

[0218] In addition to extracting the raw photometric data contained in a gelfile, the extraction method preferably also extracts the machine runtime parameters of voltage, current, and temperature, and outputs simple summary statistics of the photometric data including the minimum and maximum relative fluorescence intensity (RFI) in each filter along with a histogram of RFI values observed in each filter. This output can be and preferably is used in the optimization of fluorophore concentrations in FAFLP reactions and to detect any anomalous machine behavior.

[0219] Where the source data is an ABI sample file, it is assumed that the sample file begins with one or more specific file identifying characters. For example, ABI sample files begin with a single 128-byte file header. The first four bytes of this header (i.e., 0x41 0x42 0x49 0x46: ABIF) designate the file as an ABI processed sample file. In its current implementation, a single statically located 28-byte record, offset six bytes from the beginning of the sample file, provides information about the structure and location of dynamically located blocks of data. This table of contents or index contains a plurality of records that each provide information about the various types of data stored in the sample file. The index is made up of a group of records each referred to as a FLAG. Each FLAG record of an ABI sample file is also 28 bytes long and contains eight fields: (1) FLAG_HEADER, (2) FLAG TAG_VALUE, (3) DATA_TYPE, (4) RECORD_SIZE, (5) NUM_RECORDS, (6) LEN_RECORDS, (7) DATA_VALUE, and (8) an unknown field. The term FLAG is used herein in the context of being a type of record and therefore is used in a manner that is different from its conventional usage in computer parlance. Each FLAG contains identifiers that describe the type of data to which the FLAG points to or contains. Typically, the associated data record is stored elsewhere in the tracefile, the location of which is identified somewhere in the FLAG record.

[0220] The FLAG_TAG_VALUE consists of an integer value that distinguishes among multiple FLAG records that are identically named. For example, where a FLAG record identifies or points to photometric data, or to its location, the FLAG_TAG_VALUE field indicates which fluorophore was used or analyzed as well as whether the photometric data are raw or otherwise processed, for example, with proprietary ABI analysis software such as GeneScan. More specifically, because there are several records, typically twelve, that contain the FLAG DATA, associated FLAG_TAG_VALUEs 1-4 indicate that the associated FLAG respectively identifies or points to raw, e.g., unprocessed, spectral data for fluorophores 1-4, and associated FLAG_TAG_VALUE 9-12 indicate that the associated FLAG respectively identifies or points to ABI processed spectral data for fluorophores 1-4.

[0221] The DATA_TYPE field indicates the structure or type of data that is stored in the location to which the associated FLAG points to or identifies. For example, DATA_TYPE 0x0004 means that the data stored in the location to which the associated FLAG points to or identifies is a two-byte integer. In a current implementation, all raw or processed ABI photometric data is stored as DATA_TYPE 0x0004.

[0222] The RECORD_SIZE field indicates the length, in actual bytes, of the referenced DATA_TYPE stored in the location to which the associated FLAG points to or identifies. For example, where the data are of DATA_TYPE 0x0004 the RECORD_SIZE field will hold the value 0x0002, i.e., a two-byte integer).

[0223] The NUM_RECORDS field holds a value that indicates how many entries of a particular DATA_TYPE is present starting at the location to which the associated FLAG points to or identifies. For example, where an associated FLAG points to or identifies a location in the tracefile at or after which 100 (0x64) observations of fluorescent intensities are stored, the value of the associated NUM_RECORDS field will be (0x64).

[0224] The LEN_RECORDS field holds a value that denotes the total length, in bytes, of all NUM_RECORDS associated with a particular DATA_TYPE and FLAG. NUM_RECORDS can also be characterized as denoting the total length, in bytes, of NUM_RECORDS of referenced data. For example, where the particular data associated with a particular FLAG record has DATA_TYPE 0x0004 and NUM_RECORDS has value (0x64), the value of LEN_RECORDS will 0x00C8 (i.e., 2 bytes per record*100 records200 bytes).

[0225] In a preferred implementation, an extraction method is used to process ABI-generated sample files to extract the photometric data contained therein. First, the header content is checked for the correct content and offsets: byte 1-4 must contain the values 0x41 0x42 0x49 0x46 (ABIF) and must include a 28-byte FLAG index record located at an offset of six bytes from byte 0. If the initial header does not exhibit the aforementioned features, an error message preferably is generated. Next, each FLAG record is sorted and accessed to process record information. In a preferred implementation, only available FLAGS and associated records are sorted and displayed. The pointer index, which is a table of FLAG records, is consulted to determine the location of the photometric data contained in the sample file. For example, where the file is an ABI-generated sample file, each FLAG record that points to photometric or photometric related data, such as certain DATA FLAG records, is identified. Using the location associated with each such pointer, the DATA_TYPE, RECORD_SIZE, NUM_RECORDS and LEN_RECORDS information is used to specify and extract the photometric data to a new output file containing only the extracted data. This data is retrieved because it contains the actual photometric data in the form of preferably raw and/or processed ABI fluorophore emission data. Data contained in other sample files can also be selectively located and extracted in this same manner.

[0226] In describing the genotyping process, the forthcoming procedures are each presented with reference to a data vector as operationally defined by the dimensionality of the source data. For example, one-dimensional source data, such as photometric data obtained from a processed sample file, consists of a one-dimensional vector of photometric data. A two-dimensional data source, such as photometric data obtained from a gelfile, consists of a series of adjacent one-dimensional vectors that correspond to channels of photometric data.

[0227] Operationally, for some procedures described herein, the two-dimensional data is reduced to a series of one-dimensional vectors by applying a procedure to each data channel separately. For example, in the procedures of baseline adjustment, multicomponent correction, and spectral noise reduction, a one-dimensional form of the procedure is applied successively to each data channel, even though the original source data may have been two-dimensional. This does not mean two-dimensional applications of these procedures do not exist, or possibly even warrant merit. In fact, the use of a two-dimensional approach to estimate spectral intensity, e.g., source data from gelfiles, is well suited for being practiced using a preferred implementation of a method of the invention despite the fact that a one-dimensional approach, e.g. sample files, can be used. In a currently preferred implementation, a one-dimensional approach is preferred as it (1) provides excellent results and (2) costs less in terms computational time required to complete certain calculations.

[0228] During operation of an electrophoresis apparatus, a detector records the presence and, within the limitations of the detector, the relative intensity of emissions of fluorophore-labeled DNA fragments as they pass by a detection area. While the signal that the detector produces or outputs can be a digital signal, it typically is an analog signal that corresponds to the intensity of the fluorescence. The analog signal is subsequently converted to digital data. It is this digital data that comprises the source data.

[0229] Where an excitation device is used, its output ideally produces a response only from fluorophore-labeled DNA fragments in the polymer matrix that is read by the detector. For example, a laser-based excitation device causes fluorophores attached to DNA fragments to fluoresce and the detector is designed to detect that fluorescence. However, as is exemplified by FIG. 4, components other than fluorophores themselves may also fluoresce, therefore confounding the signal associated with the fluorophore-labeled DNA fragments. This phenomenon often leads to signal variability or error, typically in the form of spectral noise, which manifests itself as unwanted components present in the resultant detected signal.

[0230] More specifically, anything irradiated by the excitation device may produce some amount of response that can contribute to spectral noise. Where fluorophore excitation is intended to cause fluorescence, the polymer matrix that holds the DNA fragments itself can produce some aberrant spectral signals. Other components can also contribute additional spectral noise. For example, the glass-plate material holding the matrix, debris on the glass plates, capillary tubing, as well as other material in and around the matrix, can contribute noise. In most cases, the presence of this spectral noise manifests itself in an additive manner, i.e., an offset, to a signal that is only supposed to correspond to desired fluorophore emissions.

[0231] Such offsets, generally identified by reference numeral 72 in FIG. 4, have been shown to generally possess relatively uniform and constant amplitude. Where such an offset is an issue, a baseline adjustment preferably is performed to determine an offset correction that removes it completely or at least a substantial portion of it. The remaining component therefore consists primarily of the signal contributed by desired fluorophore emissions. Ideally, when baseline adjustment is complete, the signal intensity of all of the photometric data originates at about the same level and enables a better comparison of relative spectral signal intensity.

[0232] In a preferred implementation of a procedure for performing baseline adjustment, data from the electrophoresis apparatus are outputted for a certain amount of time and the data retained as source data. To the extent needed, photometric data is extracted and this data is processed to adjust its baseline.

[0233] Still referring to FIG. 4, each data point of a particular scan and filter that is processed to adjust its baseline is defined by a position that corresponds to its scan number along the X-axis, and a value indicating relative spectral intensity (RFI), that corresponds to its RFI along the Y-axis. For example, referring to FIG. 5, the baselining procedure first groups data points into windows, W₁, W₂, and W₃, of data points each having a fixed number of data points. In the preferred implementation depicted in FIG. 5, each window includes 50 data points, but each window can be larger or smaller depending on the degree of baseline removal desired.

[0234] For each window, a search is performed to determine the data point having the lowest value of spectral intensity or RFI. The slope of a line 74 between the lowest value of one window, W₂, and the lowest value of the next adjacent window, W₃, is determined, preferably by calculation. Next, the magnitude of baseline offset of each data point between the two window minima is computed. An offset value is then computed for each data point between the lowest values of spectral intensity of one window, W₂, and up to the data point with lowest value of spectral intensity of the next adjacent window, W₃. The interpolated baseline offset is subtracted from the value of spectral intensity associated with each data point. In a preferred implementation, a baseline offset is calculated for each data point using multipoint linear regression, and the calculated offset is subtracted from the value of spectral intensity of each data point.

[0235] Referring to FIG. 6, after the offset is subtracted, the values of spectral intensity preferably are shifted to a common baseline. Referring to FIG. 7, this is done for each scan and filter. An advantage of the present method of performing a baseline adjustment on the data is that the baseline is implemented dynamically. It is therefore capable of adjusting the baseline in response to changing spectral noise conditions. A preferred implementation of determining baseline-offset correction is presented in more detail below.

[0236] In single- or multi-component fluorophore-based detection systems, signal variability due to detection process often results in the accumulation of substantial spectral noise visually discernible as a signal offset from a baseline of zero spectral amplitude. To dynamically remove this offset, first consider a one-dimensional vector of spectral amplitudes (intensity responses) Y_i, i1, 2, . . . , n each indexed with corresponding predictor positions X_i, i1, 2, . . . , n . This vector is partitioned into Wn/l windows, where l represents the minimum level at which baseline-offset variability is removed (typically 50 to 100 predictor positions). In each window W_j, j1, 2, . . . , n/l, the minimum spectral intensity response W_j(Y_i^*) and its corresponding predictor position, W_j(X_i^*) is identified. The slope, _j, between minima in successively adjacent windows is then computed, $\begin{array}{cc}{W}_j\left({}_j\right)=\frac{W_{j+1}\left(Y_i^{*}\right)-W_j\left(Y_i^{*}\right)}{W_{j+1}\left(X_i^{*}\right)-W_j\left(X_i^{*}\right)}j=1,2,n/l.& \left(\mathrm{Equation}\mathrm{II}\right)\end{array}$

[0237] Next, the magnitude of baseline offset, _j^*, for each W_j(X_i^*, Y_i^*) is computed,

[0238] Using estimates of W_j(_j) and W_j(65_j^*), offset values, _i, for each response Y_i between W_j(X_i^*, Y_i^*) and W_j1(X_i^*, Y_i^*) are interpolated,

[0239] Each interpolated offset, W_j(_i) is then subtracted from each corresponding response Y_i between W_j(X_i^*, Y_i^*) and W_j1(X_i^*, Y_i^*) to shift the data to a common baseline. This procedure is performed piecewise for all W_j, j1, 2, . . . , n/l windows. It is desirable that the window l be larger than the width of any expected amplicon.

[0240] FIG. 4 illustrates an exemplary set of spectral data before baseline adjustment is performed and FIG. 6 depicts this same data after adjustment. FIG. 5 depicts source data of a single filter taken during a single scan before adjustment and FIG. 7 illustrates the same data after adjustment. The same procedure is repeated for all of the source data. The spectral intensity data will thereafter have a common baseline of zero amplitude.

[0241] It should be noted that it is preferred that baseline adjustment be performed where the photometric data are obtained from fluorophore emissions and the resulting fluorescence recorded by some style of detector. However, performing this algorithm may not be necessary where other methods or techniques of obtaining information from DNA fragments being analyzed are used.

[0242] If emission spectra from different types of fluorophores overlap, i.e., have some emission wavelengths in common, and are simultaneously recorded by the detector, it may be necessary to remove spectral overlap. This procedure is generally performed after any baseline adjustments are made. To remove spectral overlap, a multicomponent correction preferably is performed on the photometric data. An example of such overlap between the emission spectra from four different fluorophores 76, 78, 80, and 82 is shown in FIG. 8. If desired, removal of spectral overlap can be performed where fluorophore emission occurs at more than one wavelength and this multiple-wavelength emission is detected. This procedure help ensure that the photometric data subsequently analyzed consists primarily of emissions at the wavelength at which a particular fluorophore was intended to emit.

[0243] However, in cases where: (1) fluorophore emissions do not have overlapping wavelengths, (2) the fluorophore emits at only a single wavelength, or (3) the detector is capable of selectively detecting emissions at only a single wavelength with no overlap from other fluorophores, it may not be necessary to remove spectral overlap. Additionally, it may not be necessary to remove spectral overlap where the system is not a fluorescence-based detection system. A preferred implementation of a multicomponenting method that is well suited for being implemented in software is presented in more detail below.

[0244] In fluorescence-based systems, a suitable detector records excitation emissions by a plurality of different fluorophore types. Such a system is generally referred to as a multicomponent system and is generally capable of simultaneously detecting two or more distinct emission spectra. While each different type of fluorophore is chosen or designed to fluoresce maximally at a specific wavelength, fluorophores typically also fluoresce, to a lesser degree, at some range of wavelengths above and below their wavelength of maximal emission intensity. As a result, fluorophores typically have maximum emission intensity at their designated frequency and a lesser fluorescence intensity at other wavelengths, which can cause the spectrum of wavelengths in fluorophores of one type to overlap with those of another in the manner depicted in FIG. 8.

[0245] When detected, these additional emissions lead to a spectral intensity along the overlap that is greater than what it should actually be. This, of course, contributes to variability of intensity values, and it is therefore desirable to compensate for this overlap.

[0246] To compensate for spectral overlap, the fluorescence intensity components contributed by fluorophore emissions at wavelengths other than its designated frequency are removed typically by a n-fluorophore multicomponent procedure. In a multicomponent system containing n fluorophores, the fluorescence intensity s, obtained at the i^th wavelength, is the sum of the recorded fluorescence intensities at that wavelength from each fluorophore component $\begin{array}{cc}{s}_i=\underset{j=1}{\stackrel{n}{}}{m}_{\mathrm{ij}}{f}_j& \left(\mathrm{Equation}V\right)\end{array}$

[0247] where the m_ij are constant coefficients that define the degree of spectral overlap between recorded wavelengths and _j is the concentration of the j^th fluorophore. Equation V may be represented more compactly in standard matrix form,

[0248] Here, {right arrow over (s)} is a vector with n components, where n is the number of detectors at different wavelengths, and {right arrow over ()} is a vector with m components, where m is the number of fluorophore components. Thus, M is a nm matrix. In order for M to be inverted, note that n must equal m . For example, in the case of industry-standard 4-fluorophore chemistries, {right arrow over (s)} and {right arrow over ()} are vectors each with four elements and M is a 44 matrix, $\begin{array}{cc}\left[\begin{array}{c}{s}_1\\ s_2\\ s_3\\ s_4\end{array}\right]=\left[\begin{array}{cccc}{m}_{11}& m_{12}& m_{13}& m_{14}\\ m_{21}& m_{22}& m_{23}& m_{24}\\ m_{31}& m_{32}& m_{33}& m_{34}\\ m_{41}& m_{42}& m_{43}& m_{44}\end{array}\right]\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\\ f_4\end{array}\right].& \left(\mathrm{Equation}\mathrm{VII}\right)\end{array}$

[0249] The matrix coefficients vary depending on the detection instrument and the emission properties of each fluorophore.

[0250] Several published methods of determining the matrix coefficients exist. Most are based on an assumption of a linear relationship between the relative spectral intensities contained in the raw data. Once defined, the matrix M provides a transformation to obtain measures of fluorophore concentrations {right arrow over ()} from the actual signals s recorded by the detection apparatus. This process of 4-fluorophore multicomponent transformation uses the specified matrix that, when multiplied against each four-element vector of spectral data, results in a four-element vector of fluorophore concentration

[0251] The transformation performed by Equation VIII will function correctly only if the signal measurements contained in {right arrow over (s)} are approximately linear and correctly positioned to each other relative to a common baseline. Therefore, before multicomponent transformation can be applied, baseline adjustment of the spectral signal, where necessary, must first be performed.

[0252] An example of raw fluorophore emissions derived from an ABI 377XL automated sequencer is shown in FIG. 4. In a preferred implementation, using the method of Li and Speed (1999, Electrophoresis, 20, pp. 1433-1442), a 4-fluorophore multicomponent correction matrix was calculated and applied to the previously baselined data. The results of the baseline adjustment and multicomponent correction are shown in FIG. 9.

[0253] No matter what detection system is used, noise reduction, preferably spectral-based noise reduction, is performed to increase the signal to noise ratio of the photometric data. Such noise reduction is particularly desirable because it significantly reduces the impact of low amplitude background noise and can also reduce the impact of noise spikes. It does so by preferably substantially eliminating the majority of low-amplitude background noise not previously removed by any baselining or multicomponenting procedures that have been performed. Such noise reduction can also be tailored to remove or substantially reduce random noise spikes that sometimes can occur during the detection procedure.

[0254] Spectral responses recorded in single- and multi-component systems often contain considerable amounts of additional noise even after other baseline and spectral subtraction corrections have been made. For example, in fluorescent-based systems, a very common type of noise manifests itself as low-amplitude background noise comprised typically of 50 RFI or less.

[0255] Another type of noise sometimes observed in fluorescent-based systems manifests itself in the form of a random spike. A random spike typically has an unusually large RFI typically associated with only a single data point. Typically, the RFI of the spike is usually much higher in amplitude than the highest amplitude produced by any spectral emission derived from a real fluorophore present in the system.

[0256] It is believed that random noise spikes and other unwanted low-amplitude noise can also contribute to error in other types of detection systems that need not be fluorescence based. As a result of noise reduction, these spikes and other unwanted low-amplitude noise are removed or reduced significantly, thereby smoothing the spectral output of the photometric data and improving its signal-to-noise ratio.

[0257] To further enhance accurate measurement of the spectral intensity response, it is desirable to increase the signal-to-noise ratio of the data. This procedure is best performed in a manner designed to minimize the degradation of the true spectral signal while removing the majority of the noise. In its preferred implementation, a highly flexible Fourier-based transform, preferably a discrete Fast-Hartley Transform, is used to perform this noise reduction. A preferred implementation of a preferred noise reduction method is presented in more detail below.

[0258] Consider a one-dimensional vector h of real-valued spectral intensities sampled at N uniformly spaced positions (e.g. scans). Denote these samples h_n where n0, 1, . . . , N1. In this implementation, h_n, corresponds to the amplitude of spectral intensity at scan position n. The output of the Discrete Hartley Transform (DHT) is a real-valued sequence H_k given by $\begin{array}{cc}{H}_k=\frac{1}{N}\underset{n=0}{\stackrel{N-1}{}}{h}_n\cos \left(\frac{2\mathrm{nk}}{N}\right)+\sin \left(\frac{2\mathrm{nk}}{N}\right)& \left(\mathrm{Equation}\mathrm{IX}\right)\end{array}$

[0259] for k0, 1, . . . , N1. The resulting spectrum is a mapping of N discrete, time-domain spectral intensity values into N discrete frequency-domain signal components.

[0260] These discrete frequency domain signal components take the form of a series, H₀, H₁, H₂, . . . H_N1, where n is equal to the number of scans and k corresponds to the specific frequency. Thus, H₀, is the resultant component that has the lowest frequency, and H_N1 is the resultant component that has the highest frequency. The amplitude of each DHT component expresses the relative contribution of that frequency to the periodicity in the original spectral data.

[0261] After noise reduction is performed (described below), that data which remains is transformed back into the time domain using an inverse Fourier transform that preferably is an inverse Discrete Hartley Transform (IDHT). The preferred form of the inverse IDHT is given by $\begin{array}{cc}{h}_n=\underset{n=0}{\stackrel{N-1}{}}{H}_k\cos \left(\frac{2\mathrm{nk}}{N}\right)+\sin \left(\frac{2\mathrm{nk}}{N}\right)& \left(\mathrm{Equation}X\right)\end{array}$

[0262] for n0, 1, . . . , N1 and is used to convert the frequency-domain components of the DHT back into the time-domain values of spectral intensity.

[0263] After transforming the data into the frequency domain using equation IX, unwanted high-frequency noise, e.g., high-amplitude noise spikes, is removed by first computing H_k for any spectral signal sequence h and then truncating H_k at an empirically derived frequency position, T_i. The truncation position, T_i, is chosen to maximize the signal-to-noise enhancement while minimizing any distortions in h_n.

[0264] In the frequency domain provided by the DHT transformation, this preferred form of truncation is equivalent to multiplying a window function g(x) to H_k and results in smoothing of undesirably high amplitude spectral intensity data once transformed back into the time domain using the IDHT shown in Equation X. A truncation of this sort using the window function, g(x), is expressed by Equation XI below

[0265] where $g\left(x\right)=\{\begin{array}{c}1T_{0}T{T}_i\\ 0\mathrm{otherwise}\end{array}.$

[0266] Where this type of truncation is performed, I_k is substituted for H_k in Equation X above.

[0267] If it is also desired to remove unwanted low-amplitude noise, e.g., low-amplitude spectral scatter, that was present in the data before being transformed into the frequency domain, a partial truncation may be applied to the amplitude of H_k, preferably at or near it lowest frequency component, H₀, as the amplitude of this frequency component contains the vast majority of low-frequency spectral components. Performing partial truncation involves truncating a portion of the amplitude of the frequency component being so truncated.

[0268] Both types of truncation can be performed on the same set of data after transformation into the frequency domain. If desired, only frequency component truncation of one or more higher-frequency components can be performed where high amplitude noise reduction is of interest. Similarly, amplitude truncation of at least one lower frequency component can be performed where only low amplitude noise reduction is of interest. In one preferred implementation of a procedure for reducing noise in spectral data, frequency component truncation is performed before amplitude truncation. However, the order in which truncation is performed is not critical.

[0269] After truncation is performed, the remaining frequency components are transformed back into the time domain preferably by using the IDHT shown in Equation X. However, after transformation back into the time domain, either type of truncation can give rise to convolution artifacts that visibly manifest themselves as excessive broadening of the spectral intensity distribution in each distinct amplicon.

[0270] Multiplication of a Gaussian apodization function, a(x), scaled from (10) over the interval (T₀T_i) in the truncated sequence I_k, can be used to effectively attenuate these artifacts before transformation back into the time domain

[0271] where $a\left(x\right)=\{\begin{array}{c}{\exp }^{-\left(\frac{{}^2{x}^2}{4\ln 2}\right)}\\ 0\mathrm{otherwise}\end{array}{T}_{0}T{T}_i.$

[0272] Where Gaussian apodization is applied, the noise-reduced time-domain spectral signal is recovered by substituting J_k for H_k in Equation X.

[0273] Examples of this preferred noise reduction procedure is graphically depicted in FIGS. 10-17. FIG. 10 illustrates photometric data derived from a one-dimensional source before noise reduction. The X-axis depicts the scan number with the corresponding spectral intensity on the Y-axis. The data is displayed in the time domain with each data point having a value between 0 and 8191. That data has five distinct areas 84, 86, 88, 90, and 92 that exhibit spectral intensities much greater than zero amplitude. These areas contain greater concentrations of fluorophores, which correspond to the location of fluorophore-labeled DNA fragments or groupings of fluorophore-labeled DNA fragments. These locally dense concentrations of fluorophore-labeled DNA fragments are generally termed a peak (i.e., a response) of the amplicon. Between each of these peaks 84, 86, 88, 90, and 92, exists an amount of low-amplitude background noise.

[0274] In a preferred implementation of a procedure for reducing noise, the photometric data is transformed from the time domain into the frequency domain using a Fourier transform that preferably is a DHT. Equation IX above preferably is used to transform this photometric data from a time domain into a frequency domain, the graphical results of which are shown in FIG. 11.

[0275] The transformed photometric data is analyzed after transformation to determine a truncation position, T_i, that preferably will increase the signal-to-noise ratio when the transformed photometric data is retransformed back into the time domain. Where a truncation position, T_i, is determined by visual inspection, graphically displaying the DHT results can help make it easier to determine approximately where the signal-to-noise ratio approaches one. In one preferred implementation, T_i is selected at the frequency where the signal-to-noise ratio is approximately one.

[0276] Referring to FIG. 11, the transformed data visibly indicates several lower frequency components having amplitude much different than zero and many higher-frequency components having amplitude very similar to zero. In a preferred implementation, T_i is chosen to truncate at least 20% of the higher frequency components.

[0277] In the preferred implementation graphically depicted in FIG. 11, the truncation position, T_i, is set at a frequency of approximately 9.1438E-04 (scan number¹) as this has been empirically shown to provide good results where the range of each data point is between 0 and 8191. After truncation, the portion of the frequency components having a frequency greater than T_i are simply discarded. Although perhaps 9.1438E-04 is a relatively conservative truncation value, truncation at this position significantly reduces the amplitude of high-amplitude spectral noise. The result of such noise reduction after retransformation back into the time domain is graphically depicted in FIG. 14. The truncation position may differ depending upon the maximum theoretical amplitude of the spectral intensity, which may depend upon the type of detection system used.

[0278] If it is desired to reduce low-amplitude spectral noise, one or more of the lower frequency components having amplitude much different from zero can be truncated, preferably by using partial amplitude truncation. In one preferred implementation, H₀, the first frequency component, is truncated such that its amplitude is reduced by at least 25%. Preferably, it is the only frequency component that is partially truncated. In another preferred implementation, each lower frequency component sought to be truncated is multiplied by a truncation coefficient and replaced with its result. For example, in one preferred implementation, H₀ is multiplied by 0.6 and replaced with the result.

[0279] Before transformation back into the time domain, an additional operation can be performed that often provides further enhancement to the noise-reduction procedure. FIG. 12 graphically depicts an apodization finction scaled from 0 to 1 over the range T₀ to T_i. The apodization function is multiplied to the signal amplitudes remaining after any truncation has been performed on the frequency and/or amplitude components. Applying an apodization function to the truncated signal helps attenuate noise that may be introduced as convolution artifacts due to the physical truncation of the signal.

[0280] In a preferred implementation, the apodization function is of Gaussian form, which is multiplied to the remaining signal using Equation XII. In FIG. 11, the truncation position, T_i, corresponds to the specific component frequency, 9.1438E-04, at which the original signal was truncated. The apodization function is applied to the signal having a component frequency no greater than that of T_i. In the present example, FIG. 13 illustrates the truncated signal after application of the apodization finction shown in FIG. 12.

[0281] Once all manipulation of the frequency-domain components is made, the signal is transformed back into the time domain. In a preferred implementation, Equation X is used to transform the frequency components back into the time domain components of spectral amplitudes. An example of photometric data after noise reduction and retransformation back into the time domain is shown in FIG. 14. As is shown in FIG. 14, the amplitude of each peak has been slightly reduced and low-amplitude background noise between fragment peaks has been dramatically attenuated.

[0282] FIGS. 15 and 16 graphically depict the results of high-amplitude noise reduction using high-frequency component truncation. FIG. 15 shows time-domain data of a single baselined and multicomponented noise-spike 94 present at scan number 748. In the same manner previously described, the DHT was applied and the resulting signal was truncated at T_i9.7656250E-03. The Gaussian apodization function was scaled from 0 to 1 over the range T₀ to T_i and multiplied to the truncated frequency-domain data. Following transformation back into the time domain, the amplitude of the noise spike 94 is significantly reduced. This is graphically illustrated by FIG. 16, which shows the spike 94 of FIG. 15 after transformation using DHT, performing high-amplitude truncation, applying a Gaussian apodization function, and transforming the data back into the time domain using IDHT. As a result of this high-amplitude noise reduction, the impact of this noise spike 94 on the subsequent analysis discussed below preferably is minimal, if not completely negated.

[0283] This same noise reduction procedure may be applied to two-dimensional source data, e.g., a gelfile. FIG. 17A graphically illustrates the raw spectral data from a portion of a gelfile. FIG. 17B graphically illustrates the results of this noise reduction procedure applied to this same portion of the gelfile.

[0284] An advantage of performing noise reduction in this manner is that it later aids in amplicon analysis and that it is well suited for implementation in software. For example, locating the top or apex of each amplicon is made easier because the noise between amplicons is greatly reduced. The reduction in noise also makes it easier to detect the leading and trailing edges of each amplicon, thereby making it easier to more accurately estimate parameters of each amplicon, such as maximum spectral amplitude, length or width. Moreover, reducing high-amplitude noise spikes aids greatly in achieving a suitable range of spectral intensities in which to rescale the data. Taken together, all of these features greatly enhance the accurate estimation of total spectral intensity discussed herein.

[0285] In one preferred implementation, raw and/or processed photometric data are rescaled preferably to facilitate visual inspection of the data. Depending on the format in which the photometric data is initially presented, or what procedures are to be applied to the data, e.g. baseline, multicomponent, noise-reduction, etc., resealing may or may not be desirable or necessary. If rescaling is performed, photometric data are resealed into another digital representation that can be preferably more easily manipulated and used.

[0286] In one preferred implementation, photometric data are rescaled into a different digital format after applying baseline, multicomponent, and noise reduction procedures. For example, where an ABI gelfile is rescaled, two-byte (INT2) photometric data having equivalent decimal values ranging from 0 to 8191, i.e., 2¹³ values, is converted into two-byte (INT2) photometric data having equivalent decimal values ranging from 0 to 255, i.e., 2⁸ values.

[0287] Where resealing is performed, the rescaling of the photometric data into a more common and compact form is advantageous as it permits the photometric data to be more easily displayed and analyzed in visual manner. Additionally, rescaling into a more common form means that the photometric data can be used in any one of dozens of commercially available programs capable of accepting standardized image data.

[0288] It should be noted that rescaling could be performed wherever the photometric data comprises fluorescent data. However, performing this step may not be necessary where other methods or techniques of obtaining information from DNA fragments being analyzed are used.

[0289] In one preferred method of resealing, a histogram of spectral intensity values is first made of the photometric data. Thereafter, the data are scaled preferably to an industry-standard 8-bit scale such that the magnitude of the photometric data ranges from 0 to 255. To make the data easier to rescale, a portion of the largest and smallest histogram values can be removed. For example, in one preferred implementation, the largest and smallest three percent of the spectral amplitudes in the photometric data are first removed, preferably by truncation or magnitude reduction, to facilitate resealing the data into an 8-bit format.

[0290] Where the data is rescaled, it preferably is reformatted into a standard graphics format, such as bitmap (BMP), Portable Network Graphics (PNG), Portable Pix Map (PPM), Portable Gray Map (PGM), Tagged Image File Format (TIFF), Graphics Interchange Format (GIF), PC paintbrush (PCX), Joint Photographic Experts Group (JPEG), Flashpix (FPX), or another computer readable format. For example, in one preferred implementation, the rescaled photometric data is written to a Portable Gray Map (PGM) file format. In doing so, a header is added to the photometric data that denotes it as being in particular graphics data format and that includes certain dimensions of the file or other necessary information. When displayed on a screen, such as a computer driven display or the like (see below), a graphical image of the recorded spectral intensities can be depicted, such as shown in FIGS. 17A and 17B.

[0291] It is often desired to have the ability to display the photometric data at some point. For example, it is desirable to display the photometric data in order to inspect the results of the marker assay. Inspecting the results also provides a means to evaluate the progress of an experiment and also enables lane tracking to be performed, especially where two-dimensional photometric data is involved.

[0292] The majority of genetic marker-assay detection systems record the presence of a substance (response) as an ordinal response preferably on a uniformly defined time scale. Typically, these responses are measured as spectral intensities, indicating relative concentration, ranging from 0, i.e., no response, increasing to an upper value that typically is limited by the physical properties of the detector. These values may be translated into graphical representations suitable for human visualization.

[0293] Consider each response recorded by the detector as a single element of amplitude in a vector of integer-valued measurements ordered by time. Each element may then be defined as a single graphical pixel in a one- or two-dimensional array. These pixels may be displayed graphically by constructing a grayscale or false-color map that allows for visual discrimination among the spectral response data.

[0294] Due to the fact that commonly available graphical formats cannot directly display unprocessed photometric data, the data can be transformed into an appropriate form by first resealing it as described previously. For example, in the case of industry-standard four-fluorophore chemistries, after rescaling, one preferred implementation maps the four-component fluorophore intensity space into a three-component Red-Green-Blue (RGB) color space. This is accomplished by assigning four different RGB colors to each of the four fluorophoresred, green, blue, and yellow. Individual component spaces are defined by mixing the four colors according to the intensity of the corresponding pixel at each position in the pixel vector.

[0295] Where slab-gel electrophoresis is performed, multiple samples are typically loaded across one dimension of the gel (X) and simultaneously separated along the other dimension (Y). Compared to capillary electrophoresis, slab-gel electrophoresis devices do not have any intrinsic capability to precisely delimit individual samples on the gel itself. For example, standardized square-shaped loading combs provide an organization of sample wells that delimit, preferably passively, each sample with a small blank area, but only offer a simple visual means of distinguishing among adjacent sample lanes. This feature is insufficient for precisely defining DNA fragments (e.g. amplicons) present in individual sample lanes. Therefore, each sample lane, and therefore amplicons contained within that lane, must be defined by a human operator or by a computer implemented algorithm designed to perform lane tracking 50 (FIG. 1). This procedure is also necessary where a sharks-tooth comb is used, which results in almost no visually ascertainable gap between sample lanes.

[0296] In a present implementation, the task of defining the individual amplicons belonging to each sample lane is accomplished by manually constructing a tracking-spline that bisects each amplicon present in the marker assay for each sample lane. Spectral intensities are averaged laterally over a region of the sample lane, e.g., a region defined by the average width of a sample lane (cf. channel). The result of this procedure is a one-dimensional vector of averaged spectral intensities spaced over uniformly spaced time intervals (cf. sample file).

[0297] In another preferred implementation, where lane tracking is already performed and saved to a file, the tracking information is read directly from the file and later used to obtain the location of the lanes. For example, where the source data is an ABI-generated gelfile processed with ABI Genescan Analysis or Sequencing Analysis software, the lane tracking information can be extracted from the gelfile and subsequently used. In one preferred implementation, the lane tracking information is extracted and stored for later use, such as by saving it in a file. The location of each lane is obtained from the stored lane tracking information and used with the processed photometric data to show the location of the lanes when the photometric data is displayed. For example, processed and rescaled photometric data can be displayed and the lanes can be superimposed on the photometric data using the extracted lane tracking information.

[0298] Once defined, the lane-tracking data will serve as input data to enable extraction of a one-dimensional vector of photometric data from the gelfile (cf. sample file). This procedure is implemented in one of two preferred methods: (1) for each time interval or scan, the recorded value of spectral intensity (cf. channel) nearest to the tracking spline is extracted or (2) for each time interval (scan), an average value of recorded spectral intensity in a defined proximity to the tracking spline is extracted. In one preferred implementation, the proximity is defined as the number of channels defining the average width of an amplicon in the gel, typically 4 or 5 channels on an ABI 377XL automated sequencer.

[0299] Once extracted, the photometric data will be used to: (1) identify each amplicon in a sample, preferably including molecular weight standards, (2) estimate amplicon parameters such as maximum spectral intensity (amplitude), width and position, and (3) generate a sizing curve based on data derived from included molecular-weight standards. Information gathered from each procedure allows the assignment of a molecular size, either as a molecular weight or molecular mobility, to amplicons of previously unknown size.

[0300] The photometric data shown in FIG. 14 depict five points of maximally increased spectral intensity, i.e. a spectral peak, surrounded on either side by regions of reduced spectral intensity. The purpose behind identification of these areas 52 (FIG. 1) is that they indicate the presence of locally-dense concentrations of fluorophore-labeled DNA fragments. Taken together, these regions define, in FIG. 14 for example, five differently-sized amplicons. The location of the spectral peak relative to its position on the X-axis (i.e. scan number) within any single amplicon is defined to be the absolute center of that amplicon. The distinction of the peak maximum is important because even if it is assumed that all DNA fragments comprising a single amplicon are identically-sized, diffusion of these fragments does occur. Therefore, instead of occupying a single discrete location in the gel, the DNA fragments comprising an amplicon are actually spread out over a region and confound the true position of the peak. An example of the spectral distribution of a single amplicon, in one dimension, is shown in FIG. 18A. Each amplicon has associated with it a maximum height or amplitude (peak) with an inflection point or apex, and a width or variance. In the case of fAFLP, as is exemplified here, these localized regions are indicative of an amplicon, provided that the peak maximum is greater than a predefined threshold.

[0301] In a preferred implementation of a procedure for locating amplicons, values of spectral intensity are sequentially examined and compared to spectral intensities immediately preceding and following the current location. At the leading edge of a putative amplicon, i.e., where spectral intensity begins to increase steadily, each successive value of spectral intensity increases until a maximum value is reached before the spectral intensity begins to decrease. At the position (scan number) where spectral intensity first decreases, the scan number immediately prior to that position is identified as being the location of maximum spectral intensity (i.e. the inflection point or peak apex) of that putative amplicon. An advantage of this preferred implementation is that it is capable of not only quickly and simply identifying the apex of any amplicon, but it also is capable of identifying the majority of peak apexes should adjacent amplicons overlap.

[0302] This procedure also determines a putative position, typically as a scan number, of each amplicon relative to all other amplicons present in a given sample. The search algorithm also compares the value of spectral intensity at each putative position of a putative peak apex against a predetermined threshold value of spectral intensity. When the threshold is exceeded, the position associated therewith is noted and the spectral peak of the amplicon is identified. Those putative positions that do not equal or exceed the threshold are discarded.

[0303] After the peak apex has been identified, such as in the manner discussed immediately above, the trailing edge of the peak can be identified. Knowing the location of the trailing edge of an amplicon and the location of its leading edge, the amplicon full-width at half-maximum amplitude, i.e., amplicon variance, may be estimated. More specifically, in one preferred implementation, the amplicon variance preferably is calculated as the distance between the leading and trailing edges of an amplicon at the point of one-half the distance to the peak apex.

[0304] However, in a preferred implementation of a method of this invention, only the leading edge and the peak apex need be determined to provide suitable information for determination of amplicon variance. In another preferred implementation, only the trailing edge and peak apex need be determined. Preferred implementations of a procedure for locating an amplicon, its peak apex, and its variance, are discussed below.

[0305] The presence of all amplicons or DNA fragments (e.g. DNA standards) in each marker assay must first be ascertained. This is accomplished by scanning the data vector defined by the tracking-spline for amplicons containing spectral intensities greater than a defined minimum-response value C. Let be a one-dimensional vector of spectral intensities arranged by scan number. Assume that is continuous and differentiable on the interval (a, b) defined to be the beginning and end, respectively, of . Define to be increasing in (a, b) if $\begin{array}{cc}\frac{f\left(x\right)}{x}>0x\left(a,b\right)& \left(\mathrm{Equation}\mathrm{XIII}\right)\end{array}$

[0306] and decreasing in (a, b) if $\begin{array}{cc}\frac{f\left(x\right)}{x}<0x\left(a,b\right)& \left(\mathrm{Equation}\mathrm{XIV}\right)\end{array}$

[0307] Suppose that contains possibly many values c_i, such that $\begin{array}{cc}\frac{f\left(c_i\right)}{c}=0\mathrm{If}& \left(\mathrm{Equation}\mathrm{XV}\right)\\ \frac{f\left(x\right)}{x}>0x\left(a,c_i\right)\mathrm{and}& \left(\mathrm{Equation}\mathrm{XVI}\right)\\ \frac{f\left(x\right)}{x}<0x\left(c_i,b\right)& \left(\mathrm{Equation}\mathrm{XVII}\right)\end{array}$

[0308] then has a local maximum at xc_i. Provided that c_iC, each c_i defines a single spectral peak with an apex located at the corresponding position x_i. As a result of knowing the peak apex position, its maximum amplitude can also be ascertained.

[0309] Referring to FIG. 18, FIG. 18B illustrates a portion of the leading edge of a portion of the amplicon shown in FIG. 18A beginning at scan numbers 623 and 635, respectively. Although the peak apex is not shown in FIG. 18B, the vertical line centered at position X_LM corresponds to the location of the peak apex, i.e., the position of the peak maximum amplitude, with magnitude Y_LM. The solid circles indicate observed values of spectral intensity at their corresponding scan numbers, X_(L) and X_(L1). The open circle shows the intersection of an imaginary dashed line drawn perpendicular to the peak apex position centered at one-half the maximum amplitude of the peak apex, Y_(LM/2), located at position X.

[0310] To obtain an estimate of the amplicon variance (full-width at half-maximum amplitude), first define the quantity Y_LM as the local maximum amplitude at position c_iX_LM of the amplicon. An approximation of the position X where a line originating at (X_LM, Y_LM/2) perpendicular to (X_LM, Y_LM) crosses the leading edge of the amplicon is given by: $\begin{array}{cc}{\hat{X}}_{}=X_{\left(L+1\right)}-\left(\frac{\left(Y_{X_{\left(L+1\right)}}-\left(Y_{\mathrm{LM}}/2\right)\right)}{\left(Y_{X_{\left(L+1\right)}}-Y_{X_L}\right)}\right)& \left(\mathrm{EquationVIII}\right)\end{array}$

[0311] where Y_XL(Y_LM/2)Y_X(L1) is constrained by the condition (X_L1X_L)1. Assuming roughly symmetrical amplicon proportions, a good estimate of the amplicon full-width at half-maximum amplitude (FWHM) is given by: 2(X_LMX).(Equation XIX)

[0312] The calculated estimates of peak apex position, X_LM, maximum amplitude, Y_LM, amplicon full-width at half-maximum amplitude, , provide starting values required by the nonlinear amplicon-parameterization procedure described below.

[0313] A preferred implementation of a procedure for detecting a leading edge of an amplicon involves scanning through the data vector and determining the presence of a plurality of consecutively increasing spectral intensities. Where the presence of a putative amplicon is identified, it is accepted as an amplicon if the peak apex is above a threshold spectral intensity value. Otherwise, it is rejected. Such a preferred implementation helps minimizes the occurrence of false negatives and false positives. In cases where this procedure is implemented in software, a message preferably is generated.

[0314] In one preferred implementation, a leading edge is confirmed when five spectral intensities are observed to consecutively increase. An amplicon is accepted if the maximum spectral intensity of the series of increasing spectral intensities detected is greater than a peak spectral intensity threshold, for example, of 50 RFI. If desired, a greater number of consecutively increasing spectral intensities can be used to detect the presence of a leading edge where the risk of false positives is desired to be very small. If desired, the threshold spectral intensity value can be higher or lower depending on the needs of the user.

[0315] In another preferred implementation, a check is made two determine whether the spectral peaks of two amplicons are located relatively close together. For example, where two spectral peaks are detected within 10 datapoints or scans of each other, the location (scan number) of each peak is identified and a message preferably is generated. Identifying incidences of two closely spaced spectral peaks enables a more detailed manual examination or automated analysis of this region to be made as such an occurrence is indicative of amplicon overlap. A second check preferably is also performed to determine the presence of flat maxima, that is, a region where spectral intensities are the same as the maximum spectral amplitude that the detector is capable of detecting or outputting. A message preferably is also generated where a flat maxima occurs so that this region of data can be more closely examined and disregarded, if necessary, as these values are prone to extreme amounts of error.

[0316] Despite having attenuated significant amounts of noise and possibly having performed other procedures to enhance the quality of the photometric data, there may still be appreciable amounts of error in the location, amplitude, and width of a DNA fragment or amplicon. For instance, slight variations in the polyacrylamide gel matrix may cause identically-sized DNA fragments or amplicons to travel through the gel at different rates. This can cause location error, which can ultimately lead to error in assigning a molecular size, either, e.g., as a molecular weight or molecular mobility, even when the DNA fragment or amplicon is compared to a known standard.

[0317] Other types of error may also exist. For example, error in spectral intensity measurement manifesting itself as amplitude error may lead to error in amplicon variance, i.e. full-width at half-maximum amplitude. Fragment location error can also occur because the detector has limits to the amount of data it can record. As a result, for example, the apex of an amplicon can be completely missed by the detector, causing its location to be erroneously estimated. Other limitations in the detection process, including scanning error and frequency, can also result in incomplete signal capture and lead to such errors.

[0318] To help correct such error, it is desirable to fit an idealized function or curve to each amplicon identified in the source data. For example, where the system is a fluorescence-based system, the data comprises photometric data representing spectral intensities. The fitted function essentially replaces or supplants each amplicon defined by the observed data. This helps recapture and/or restore lost signal and can help correct location error.

[0319] In one preferred implementation of an algorithm for providing additional error correction, it is assumed that the distribution of spectral intensity in each identified amplicon can be approximated as being Gaussian in shape, and a Gaussian function is therefore fitted to each amplicon. Where the apex of an identified amplicon has somehow been missed or skewed, fitting such an ideal and continuous curve to data defining the amplicon can recover lost signal at the apex, while also helping to ascertain more accurately its true location. In fact, depending upon the resolution of the detector, the fitting of a continuous curve or function to an identified amplicon that is defined by discrete digital actual data can help more accurately define how the amplicon should have ideally been detected. The result is a continuous curve or function that adjusts itself to the observed data in a manner that more accurately reflects its expected signal amplitude value at each location of the amplicon.

[0320] In one preferred implementation, estimates of apex amplitude, apex position, and variance of an amplicon derived from the implementation of the amplicon response identification procedure disclosed in the preceding section, are used as starting parameters for the Gaussian-function fitting procedure. Preferably, a nonlinear, orthogonal distance regression (nODR) fitting algorithm is used. Such an algorithm assumes that there is error not only in the observed value of spectral intensity (i.e. the response variable) but also the position of the response (i.e. the predictor variable). This methodology differs markedly compared to traditional regression procedures which assume no error in the predictor variable. Since it as an objective of the invention to more accurately and precisely define the location of an amplicon, the use of nODR is clearly advantageous. In a preferred implementation, the location of the spectral peak corresponds to a scan number and is defined as the X variable that is also referred to as being the predictor variable. The amplitude at the apex position corresponds to spectral intensity and is defined as the Y variable that is also referred to as being the response variable.

[0321] Such a nonlinear fitting algorithm helps correct for or compensate for peak apex location, amplicon variance, and spectral amplitude errors. In one preferred implementation, a suitable nonlinear least-squares orthogonal distance regression approach, such as that set forth in Equations XX to XXIV below, preferably is implemented using ODRPACK. ODRPACK is a configurable subroutine library that is publicly available from the repository at NetLib at www.netlib.org. A preferred implementation of an algorithm that fits a Gaussian curve to each peak using orthogonal distance regression is discussed in more detail below.

[0322] All amplicons derived from the marker assay must be assigned a unique value so that homology, e.g., as defined by its size, may be assigned unambiguously across multiple sample assays. This is accomplished by defining a relationship between the molecular weight and/or molecular mobility of a DNA fragment (or an amplicon) in a gel and its size. Size-standardized DNA fragments included in each lane with the products of the marker assay form the basis by which size-unknown amplicons are estimated.

[0323] Let (X_i, Y_i), i1, 2, . . . , n be a one-dimensional vector of observed spectral intensity data with true, i.e. unknown, values (x_i, y_i), i1, 2, . . . , n. Assume that _i and _i is the random error associated with (x_i, y_i) is

[0324] Assume that the values of the spectral intensity response y_i are a nonlinear function of the position predictors (e.g. scan numbers) x_i and a set of unknown parameters . Then the observed value of y_i satisfies y_i(x_i; ) or

[0325] for some actual but unknown values, *, *(₁, ₂, . . . , _p) where p is the number of parameters to be solved. The explicit orthogonal distance regression (ODR) procedure approximates * by defining an orthogonal distance r_i from the point (X_i, Y_i) to the curve ({tilde over (x)}, {tilde over ()}) where indicates an unspecified value of the specific variable. Given that the errors are normally distributed and have zero-mean, i.e. _iN(0, _i²) and _iN(0, _i²), the minimization of the sum of squares of r_i²{tilde over ()}_i²{tilde over ()}_i² is the maximum likelihood estimate of *. Thus, * is the solution of $\begin{array}{cc}\underset{{}^{*},\stackrel{~}{},\stackrel{~}{}}{\min }\underset{i=1}{\overset{n}{}}\left[{\stackrel{~}{}}_i^2+{\stackrel{~}{}}_i^2\right]& \left(\mathrm{Equation}\mathrm{XXII}\right)\end{array}$

[0326] subject to the constraint

[0327] Since _i in Equation XXIII is linear, _i can be eliminated thereby resulting in an unconstrained minimization problem $\begin{array}{cc}\underset{{}^{*},\stackrel{~}{}}{\min }\underset{i=1}{\overset{n}{}}\left({\left[f\left(X_i+{\stackrel{~}{}}_i;{}^{*}\right)-Y_i\right]}^2+{\stackrel{~}{}}_i^2\right).& \left(\mathrm{Equation}\mathrm{XXIV}\right)\end{array}$

[0328] In a preferred implementation, it is assumed that the distribution of fluorophore intensity in each fully resolvable amplicon is closely approximated by a Gaussian function $\begin{array}{cc}{Y}_{i}f\left(X_i;\right){}_1\exp \left[-\left(\frac{{\left(X_i-{}_2\right)}^2}{2{}_3^2}\right)\right]& \left(\mathrm{Equation}\mathrm{XXV}\right)\end{array}$

[0329] where ₁ describes the observed amplitude at the peak apex, ₂ is a location parameter corresponding to the peak apex position, and ₃ describes the amplicon variance that preferably is full-width at half-maximum apex amplitude. For multiple amplicons in the same data vector, Equation XXV may be summed and solved simultaneously over all p amplicons of interest, $\begin{array}{cc}{Y}_{i}f\left(X_i;\right)\underset{=1}{\overset{3_p}{}}{}_1\exp \left[-\left(\frac{{\left(X_i-{}_2\right)}^2}{2{}_3^2}\right)\right]++{}_{3p-2}\exp \left[-\left(\frac{{\left(X_i-{}_{3p-1}\right)}^2}{2{}_{3p}^2}\right)\right].& \left(\mathrm{Equation}\mathrm{XXVI}\right)\end{array}$

[0330] Estimates of the required starting values for each of the 3p parameters in Equation XXVI, i.e., ₁, ₂, . . . , _3p, are obtained by preferably first fitting Equations XXV or XXVI to the data using an ordinary nonlinear least-squares (nOLS) procedure (i.e. assuming no error is observed in the predictor variable). Estimates of starting values for each of the 3p parameters required for this fit preferably are provided by the amplicon parameter estimates derived preferably from the initial amplicon identification step described above. Upon convergence of the nOLS algorithm, such estimates of are then used as starting values for fitting Equations XXV or XXVI to the data using a nODR procedure (i.e. assuming the predictor variable contains error). By first estimating the 3p parameters with the nOLS procedure, the amount of time necessary for convergence of the nODR algorithm is reduced significantly. The estimates of from the nODR procedure are also used to preferably generate a molecular sizing curve in a manner such as described below.

[0331] The above methodology can also preferably be implemented for the two-dimensional case except that two additional parameters are required. The first describes amplicon position in the second dimension and the other describes amplicon variance in the second dimension. Where such is desired, estimates of peak apex height, ₁, peak position in the first dimension (e.g. X), ₂, amplicon variance in the first dimension (e.g. X), ₃, peak position in the second dimension (e.g. Y), ₄, and amplicon variance in the second dimension (e.g. Y), ₅, can be obtained using a version of Equation XXV or XXVI modified for the two-dimensional case. For example, Equation XXV modified for the two-dimensional case becomes $\begin{array}{cc}{}_1\exp \left[-\left(\frac{{\left(X_i-{}_2\right)}^2}{2{}_3^2}\right)\right]\exp \left[-\left(\frac{{\left(Y_i-{}_4\right)}^2}{2{}_5^2}\right)\right].& \left(\mathrm{Equation}\mathrm{XXVII}\right)\end{array}$

[0332] An example depicting a preferred one-dimensional implementation is graphically represented in FIGS. 18-19. FIG. 18A illustrates a baselined, multicomponented, and noise-reduced spectral peak 96. FIG. 18B illustrates a portion of the amplicon used to determine a preferred variance, namely, the full-width at half-maximum. FIG. 18C illustrates an nODR fit of a Gaussian function 98 (in phantom) overlaid on top of the amplicon shown in FIG. 18A. Also displayed in FIG. 18C are the estimates of the three fitted parameters of peak amplitude, AMPLITUDE, position, X POSITION, and amplicon full-width at half-maximum amplitude, VARIANCE. For subsequent analyses, the fitted Gaussian amplicon 98 may or may not replace the original amplicon 96, depending on the needs of the user. For example, when constructing an amplicon sizing curve, in replacing the original peak 96, the fitted values of Y_i generated using Equations XXV or XXVI as a function of X_i preferably replace the observed (X_i, Y_i) data values.

[0333] In another preferred implementation, parameters of multiple amplicons may be simultaneously estimated using a sum of multiple Gaussian functions, preferably using Equation XXVI. FIG. 19A illustrates five baselined, multicomponented, and noise-reduced amplicons 100, 102, 104, 106 and 108. FIG. 19B illustrates a nODR fit of five Gaussian functions 110, 112, 114, 116 and 118 (in phantom) overlaid on top of the amplicons shown FIG. 19A. Also displayed in FIG. 19B are the estimates of the fitted parameters of peak apex amplitude, position and variance for each amplicon.

[0334] In an alternately preferred implementation, a starting value for amplicon variance (e.g. with reference to Equation XXV, parameter ₃) is provided from a table of estimated variances based on peaks observed in a DNA molecular-weight standard. The table includes the computed variance for each standard used. For example, where 22 standards are used such as with DNA standard ILS-600, the table contains a computed variance for each DNA standard. As such, the approximate location of the amplicon being fitted can be used to select an estimated variance using the DNA standard closest to that location. In a preferred implementation, linear regression is used to obtain a variance estimate where the peak apex location of the amplicon being fitted lies between two of the DNA standards. Thereafter, this amplicon variance estimate used as a starting value at the initiation of either minimization routine. In the preferred implementation set forth above, this variance estimate is input as variable ₃ into, for example, Equation XXV above.

[0335] Although the use of an estimated variance value that is based upon a standard rather than an actual amplicon may comprise an estimate that is not exact (especially as the variance is dependent on the amount of DNA deposited onto the gel), it appears accurate enough to allow the nOLS minimization routine to more quickly converge to an accurate estimate of variance than by simply guessing. This is due to the fact that estimates of amplicon peak apex location and amplitude are generally much more accurate, thereby reducing the necessity of an extremely accurate starting value of estimated amplicon variance. The amplicon variance estimate obtained by using estimates obtained from DNA standards is simply a reasonable estimate that helps minimize computational time of the minimization routine implemented in software, while still ensuring convergence to an accurate estimate of variance for the amplicon being fitted.

[0336] In another preferred implementation, estimates of amplicon variance used for starting values in the nOLS and/or nODR least-squares minimization routines are computed directly from the full-width at half-maximum amplitude value in the amplicon response identification procedure discussed in the preceding section. Preferably, the estimate of amplicon variance used as a starting value in the nODR least-squares minimization routine is the full-width at half-maximum amplitude value obtained from the initial nOLS fit or the amplicon peak identification procedure.

[0337] In another preferred implementation, the nODR procedure is used to estimate amplicon parameters of position, variance, and amplitude for each amplicon individually. Starting values for amplicon peak amplitude, position, and variance for the minimization procedure are obtained from the amplicon response procedure described above. Once all fitted parameters for each individual amplicon is obtained, these values can be used for starting values in the minimization procedure where parameters of all amplicons taken together are estimated.

[0338] This preferred fitting procedure advantageously helps minimize variability when binning amplicons based on their estimated sizes. This increases not only the certainty of ascertaining an accurate bin for each amplicon, but does so by increasing the accuracy of the molecular size obtained. For example, as a result of using this preferred fitting algorithm, molecular sizes of fragments associated with each amplicon can typically be estimated to a precision of at least 0.333 nucleotide bases as discussed in more detail below.

[0339] After nonlinear orthogonal distance regression (nODR) has been performed to obtain the best-fit parameters of peak amplitude, location, and variance for each amplicon, a sizing function or curve preferably is generated with DNA standards of known size. This preferably is done using DNA size standards that are electorphoresed simultaneously with amplicons from each sample. The function or curve is used to estimate the size of size-unknown amplicons as either, but not limited to, (1) a molecular weight, e.g. in units of nucleotide bases or (2) a molecular mobility, e.g. in units of scan numbers. By using molecular weight or mobility data from DNA size standards electrophoresed concurrently with size-unknown sample amplicons, the molecular weights or mobilities of such amplicons may be estimated 54 (FIG. 1).

[0340] In a preferred implementation that uses molecular weights to estimate the sizes of amplicons, 22 DNA fragment standards ranging in size from 60 to 600 nucleotide bases, such as with ILS-600, comprise the DNA-standard sample. Prior research has demonstrated that an approximately linear relationship exists between the molecular weight and the distance traveled (mobility) of single-stranded DNA fragments in a gel matrix. This approximate linear relationship has been used in the past in the form of a sizing-curve that is in turn used to extrapolate a molecular weight for any amplicon or DNA fragment of unknown size, within the size limits of the DNA standard used.

[0341] However, while this assumed linear relationship may indeed apply to DNA fragments having molecular weights within the range of approximately 100 to 400 nucleotide bases in length, it is believed that this assumption does not apply to DNA fragments having very large or very small molecular weights. As a result, estimates of molecular weights of DNA fragments outside of this range that are made assuming that such a linear relationship exists can be erroneous. Moreover, recent experimental evidence suggests that the common belief of an approximately linear relationship among molecular weight and mobility of mid-sized DNA fragments may not be very accurate.

[0342] To remedy such error, a preferred procedure for generating a sizing function of this invention posits that a locally-quadratic relationship applies to a small region surrounding each DNA-standard fragment. More specifically, in a preferred implementation, a locally-quadratic, weighted regression approach is used to determine a sizing function that accounts for any nonlinear migration of DNA fragments and or amplicons in a gel. In this preferred approach, locally-quadratic regression is performed on subsets of DNA standards smaller than the entire complement of all DNA standards. The resultant function produces a curve that helps to more accurately and more precisely estimate the molecular weights of size-unknown amplicons.

[0343] In a preferred implementation, for each DNA fragment standard, locally-quadratic, weighted regression is performed piecewise on successive subsets of all DNA fragment standards. Preferably, the locally-quadratic fit is applied to each DNA standard such that the weighted regression is performed using less than the entire complement of DNA standards. Preferably, approximately eight or nine DNA standards are used for each locally-quadratic fit performed.

[0344] In a preferred implementation, approximately eight standards are used for each locally-quadratic fit. For example, where 22 standards are used such as with ILS-600, a locally-quadratic fit is applied to each individual DNA standard using the nearest eight DNA standards. This pattern is repeated piecewise until such a locally-quadratic fit is generated for all 22 DNA standards. As is depicted by FIG. 20, the resultant fitted data can then be used to generate a nonlinear mobility-compensated sizing curve 120.

[0345] By performing such an analysis for the DNA size standards, the resultant sizing curve generated advantageously accommodates nonlinear or anomalous DNA fragment migration while helping to achieve a high degree of molecular weight prediction for sample amplicons with precision typically less than 0.333 nucleotide bases. In fact, the level of precision achieved is typically between 0.05 and 0.25 nucleotide bases. A preferred implementation of an algorithm for generating a sizing curve is discussed in more detail below.

[0346] To create a sizing function suitable for assigning molecular weights to amplicons of unknown sizes, a weighted, locally-defined quadratic weighted regression model is applied to defined neighborhoods of size-standardized DNA fragments. For example, at each DNA fragment standard position in the set of DNA standards, a polynomial of degree two (i.e. a quadratic polynomial) is fit to a subset of the entire complement of DNA standards near the fragment standard whose response is being estimated. A quadratic polynomial is fit, using weighted least-squares procedures, giving more weight to standard fragments near the DNA standard fragment whose response is being estimated and less weight to standard fragments further away. The value of the local regression function for the DNA standard being fitted is then obtained by evaluating the local polynomial using the explanatory variables computed for that standard. The regression procedure is complete after regression functions have been computed for each of the standards in the entire set of DNA standards.

[0347] The sizing function or curve is constructed by first choosing a smoothing parameter that denotes the fraction of DNA standards used in the computation of each locally-defined quadratic fit. The parameter is called a smoothing parameter because it controls the behavior of the regression function. For example, large values of produce very smooth functions. Small values of will more closely conform to fluctuations in the data. Using too large or too small a value of the smoothing parameter is not desirable either because the regression function will not reflect accurately subtle differences in fragment mobility (i.e. overfit the data) or begin to model the random error in the data, respectively. Generally, useful values of include those values of between 0.25 and 0.50.

[0348] Let qN where q is the number of DNA standards that define the neighborhood-size of the weighting function and N is the total number of DNA standards. Preferably, is a real number between (d1)/N and 1.0, where d is the degree of the local polynomial. The value of , and therefore the value of q, is dependent on the number of DNA standards used in defining the sizing curve. For a particular number of DNA standards, preferably is chosen to minimize the residual standard error (RSE) of the weighted least-squares fit of Equation XXX below as applied piecewise to the entire set of DNA standards. For example, where the DNA standard is an ILS-600 standard, preferably is between 0.136 and 1.0. Using the criteria of RSE minimization as described above, a value of of approximately 0.4 results in a RSE of about 0.094 which is significantly less compared to any RSE resulting from an such that 0.136<0.4 or 0.4<1. Where twenty-two standards are used, q is 8.8, meaning that, after rounding down to the nearest integer, eight DNA standards will define the size of the local neighborhood in which the weighted quadratic regression is performed. Where a different DNA standard is used, routine testing and experimentation can be used to determine a suitable for that standard.

[0349] As indicated above, the weighting function gives the most weight to the DNA standards nearest to the standard position currently being estimated and reduces the weight of standards more distant from the center of the neighborhood. The use of the weighting function is based on the idea that DNA standards nearer to each other are more likely to be related to each other (e.g. in terms of mobility) in a simpler way than DNA standards much further away. Thus, standards that are most likely to follow the local model best influence the local model parameters the most. Standards that are less likely to follow the local model (because they are further away from the center of the neighborhood) influence the local parameter estimates less.

[0350] Define x_i as the i^th DNA standard fragment where i1, 2, . . . , N. Let d_i be the distance (e,g, in units of scan number) from the estimated peak apex location of a standard fragment being fitted to its q^th nearest neighbor. The estimated peak apex location of the DNA standard fragment preferably is obtained using ODR as described above. Define the weighting function W(u) as $\begin{array}{cc}W\left(u\right)=\{\begin{array}{cc}{\left(1-|u{|}^3\right)}^3& 0u<1\\ 0& \mathrm{otherwise}\end{array}.& \left(\mathrm{Equation}\mathrm{XXVIII}\right)\end{array}$

[0351] The function W(u) imposes a weight to each standard fragment in the neighborhood as a decreasing function of its distance from the center of the neighborhood. The center of the neighborhood, in this instance, is the location, x_i, of the standard currently being fitted. The weight contributed by the nearby standards (x_k, y_k), where k1, 2, . . . , q, to that of the standard at location x_i being fitted is thus defined as $\begin{array}{cc}{w}_i\left(x_k\right)=W\left(\frac{\left(x_i-x_k\right)}{d_i}\right)& \left(\mathrm{Equation}\mathrm{XXIX}\right)\end{array}$

[0352] where x_k is the apex location of the nearby standard and y_k is its known (standardized) molecular weight. The weight for a specific point in any neighborhood is obtained by evaluating the weight function at the distance between the point x_k and the point of estimation, x_i, after scaling the distance so that the absolute distance over all the neighborhood points is exactly one. The fitted value for the DNA standard located at x_i is computed by the ordinary least-squares minimization of $\begin{array}{cc}\underset{k=1}{\overset{N}{}}{w_i\left(x_k\right)\left[y_k-\left({}_1+{}_2{x}_k+{}_3^2{x}_k\right)\right]}^2& \left(\mathrm{Equation}\mathrm{XXX}\right)\end{array}$

[0353] where w_i(x_k) is the weighting function of Equation XXIX. Once the best-fit estimate of parameters ₁, ₂ and ₃ are obtained (note that the parameters ₁, ₂ and ₃ of Equation XXX have no direct relationship to similarly denoted parameters in the ODR section), the molecular weight of size-unknown amplicons in each defined neighborhood are estimated by

[0354] where ₁ is the estimated molecular weight of the i^th size-unknown amplicon, x_i is the peak location estimate of that same amplicon that preferably is obtained by ODR as described above, and {circumflex over ()}₁, {circumflex over ()}₂ and {circumflex over ()}₃ denote the least-squares estimates of ₁, ₂, and ₃ obtained from minimizing Equation XXX.

[0355] FIG. 20 illustrates an example of one such sizing curve generated using this procedure. The curve shown in FIG. 20 is produced from the results of plotting the results of Equation XXX for each of the standards in the ILS-600 DNA standard.

[0356] After estimating a molecular weight for each size-unknown sample amplicon, a binning operation to group amplicons of similar size is preferably performed. This binning operation preferably is done because, although the aforementioned methods produce a highly accurate estimate of molecular weight for each size-unknown amplicon, some size variability nonetheless remains. This variability is manifested usually as very small differences in the estimated molecular weight of amplicons derived from the same locus thus preventing direct assignment of homology across samples. It is therefore desirable to bin or group similarly-sized amplicons and thus define groups of amplicons based on homology. It is posited that amplicons exhibiting estimated molecular weights that differ only slightly are in fact identically sized, (i.e. homologous), across samples. Preferably, the binning procedure eliminates any remaining variability attributable to small differences, typically less than 0.333 nucleotide bases, in the estimated molecular weight. This procedure also has the advantage of providing a means by which to assign amplicons to discrete and homologous groups across multiple samples. A preferred implementation of a procedure for binning the molecular weights of amplicons is discussed in more detail below. Other binning procedures can be used.

[0357] Although each amplicon now has an estimated size assigned to it, these molecular weights are numerically continuously valued. In other words, due to variability in assigning a molecular weight to any amplicon, the estimated size of the amplicon is unlikely to be an exact integer reflecting the precise size of the amplicon, for example, in nucleotide bases. Therefore, the estimated amplicon sizes must be transformed into a discrete representation of an amplicon at a specific locus. This procedure is necessary to make assign homology to homologous amplicons across other samples.

[0358] In a preferred implementation, for each processed sample, S_i, i1, 2, . . . , N, the estimated molecular weights of the j^th amplicon S_ij, j1, 2, . . . , n_i^*, are arranged in ascending order $\begin{array}{cc}{S}_{\mathrm{ij}}=\{\begin{array}{c}{S}_{11}<S_{12}<<S_{1j}\\ S_{21}<S_{22}<<S_{2j}\\ \\ S_{\mathrm{i1}}<S_{\mathrm{i2}}<<S_{\mathrm{ij}}\end{array}& \left(\mathrm{Equation}\mathrm{XXXII}\right)\end{array}$

[0359] where N is the total number of samples and n_i^* is the total number of amplicons in each sample S_i. Due to the possibility of polymorphism between samples, for example, in the case of a homozygote null in one or more samples, the value of n* for a specific sample S_i does not have to be equal over all N samples. For each amplicon represented as a molecular weight in S_ij, any two amplicons with molecular weights, S_ij, S_ij^*, are defined to be identically-sized if

[0360] where is the largest residual of all molecular weight estimates derived from Equation XXX above. To transform these molecular weights into discrete representations of homologous amplicons, a new matrix M is created and mapped onto the matrix S_ij where each element in M_ij is defined as $\begin{array}{cc}{M}_{\mathrm{ij}}=\{\begin{array}{cc}1& \left|S_{\mathrm{ij}}-S_{\mathrm{ij}}^{*}\right|\\ 0& \mathrm{otherwise}\end{array}.& \left(\mathrm{Equation}\mathrm{XXXIV}\right)\end{array}$

[0361] The end result of this procedure is a matrix of ones and zeros indicating whether or not any amplicon represented by its molecular weight in S_ij has a homologous representation in a given sample i at position j in matrix M_ij. Because the presences of all homologous amplicons across all samples are now accounted for, the value of n_i^* must be equal across all N samples in matrix M_ij.

[0362] The j_th column, j1, 2, . . . , n*, in matrix M_ij may be conveniently referred to as a locus using classical genetic terminology. Defined in this manner, matrix M_ij is equivalent to an indicator of amplicon phenotype (i.e. presence or absence) at each locus. Because the presence of an amplicon is indicative of either a homozygote present or heterozygote genotype, matrix M_ij is therefore a dominant representation for each amplicon genotype.

[0363] In M_ij, any locus j,j1, 2, . . . , n* is defined putatively monomorphic, pMM_k, k1, 2, . . . , m, if at all i, M_ij1. If at any i, M_ij1, then locus j is defined polymorphic, PM₁, l1, 2, . . . , p. In a preferred implementation, each locus is uniquely defined by computing the average molecular weight of all amplicons present at that locus. This procedure preferably is accomplished by first referencing the presence of each amplicon in M_ij(i.e. each M_ij1) back to the original molecular weight table, S_ij. When computing the average molecular weight for each locus, null homozygotes in S_ij (i.e. each M_ij0) are excluded from the calculation.

[0364] As described above, the matrix M_ij provides information concerning only the presence and/or absence of amplicons at a given locus in a given sample. This information is limited to identifying only the phenotype (i.e. presence and/or absence) of the amplicon and is therefore dominant data as genotypes for each amplicon remain unknown. The procedures described below provide a preferred methodology to enable the recovery of genotypic information previously confounded as an expression of dominant-only phenotypic data. The overall result of this procedure is a transformation of the matrix of dominant-only phenotypic data into a matrix of co-dominant genotypic data where each amplicon has assigned to it a specific genotype.

[0365] Referring additionally to FIGS. 21A, 21B, 22A and 22B, during gel electrophoresis, individual DNA molecules in a gel or polymer matrix undergo random movements over time due to both concentration and electrically-mediated displacement effects. Because of random statistical fluctuations in these complicated processes, the DNA molecules experience directional and displacement asymmetry, so that they spread outward from a uniformly distributed starting point. In this way, an initially very narrow band of identically-sized DNA molecules (e.g. an amplicon) widens into a characteristic Gaussian-shaped profile. See, e.g., FIGS. 17A-19B and FIGS. 21A-22B.

[0366] This physical phenomenon presents a convenient situation in which to mathematically model the distribution of spectral intensity in an amplicon and therefore estimate the amount of DNA molecules present in the amplicon. Because each DNA molecule is labeled only with a single fluorophore, and only fluorophore-labeled DNA molecules are detected, the total spectral intensity will be directly proportional to the actual number of fluorophore-labeled DNA molecules in the amplicon. Thus, the total amount of spectral intensity observed in an amplicon provides a robust measure of the total number of fluorophore-labeled DNA molecules in that amplicon. The ability to estimate the relative number of DNA molecules in an amplicon allows inferences concerning the number of alleles from which the PCR-generated amplicon was originally derived. With this information, genotypes may be accurately assigned to individual amplicons.

[0367] It should be noted that others have tried to estimate the total spectral intensity of amplicons by simply adding the recorded signal intensities of fluorophore-labeled DNA molecules. This approach is not always preferred because: (1) numerous sources of spectral noise exist and contribute significantly to spectral-intensity measurement error, (2) most fluorescence detectors have relatively low sampling rates and provide only a coarse sample of spectral emissions at any given position in the gel, such as is depicted in FIG. 22A, (3) the data available for use may have been processed using unknown procedures, and (4) no intrinsically-defined method to delimit the boundaries of an amplicon exist other than perhaps an arbitrary cutoff value.

[0368] A procedure that preferably is used in a method of the invention, described below, can be used to overcome these limitations. First, a mathematical model is provided that describes the effects on an amplicon of (1) passive diffusion and (2) the additional additive effects of electrically-mediated displacement. Based on this model, a parametric statistical model is then presented that preferably is suitable for use as a procedure implemented in software that estimates the total amount of spectral intensity present in an amplicon 56 (FIG. 1). This procedure is ultimately well suited for use in estimating the spectral intensity of an amplicon and therefore the relative number of fluorophore-labeled DNA molecules in that amplicon. Ultimately, the relative number of DNA molecules in an amplicon is used to infer allele dosage (2, 1, or 0 copies), which in turn is a direct indicator of genotype (homozygote present, heterozygote, or homozygote absent null).

[0369] It is assumed that the gel or polymer matrix in which passive diffusion of a single-stranded DNA fragment occurs is isotropic and effectively either one-dimensional, such as where capillary electrophoresis is used, or two-dimensional, such as where a slab gel electrophoresis is used. The change in concentration from an initially very narrowly deposited source of identically-sized DNA molecules (e.g. an amplicon deposited into a sample well) is described in one dimension (x) by $\begin{array}{cc}\frac{C}{t}=D_p\left(\frac{{}^{2}C}{x^2}\right)& \left(\mathrm{Equation}\mathrm{XXXV}\right)\end{array}$

[0370] and in two dimensions (x and y) by $\begin{array}{cc}\frac{C}{t}=D_p\left(\frac{{}^{2}C}{x^2}+\frac{{}^{2}C}{y^2}\right).& \left(\mathrm{Equation}\mathrm{XXXVI}\right)\end{array}$

[0371] For the above equations, C is the concentration of the diffusing DNA molecules, t refers to time, and D_P is a constant describing the passive diffusion of single-stranded DNA molecules in the defined isotropic matrix. For the one-dimensional diffusion case of Equation XXXV, C is a function of scan number or location x (e.g. the direction of electrophoresis), and x, y (e.g. the dimension containing the sample wells and the dimension of electrophoresis, respectively), for the two-dimensional diffusion case of Equation XXXVI. For the two-dimensional case, x specifically refers to the channel numbers spanned by each sample well, and y refers to scan number or location along a particular channel in the direction of electrophoresis. Preferably, D_P, defines the passive diffusion rate of single-stranded DNA molecules in a polyacrylamide gel or other defined gel matrix. For discussion contained herein, D_P is assumed to be constant for a given length of ssDNA. Methods of determining D_P are known in the literature and are dependent upon factors that include: (1) the porosity of the gel matrix, (2) the strength of the electric field used, and (3) the length of the DNA molecule. Routine experimentation can be performed to determine D_P and estimates are also available in the literature.

[0372] For the one-dimensional case, the initial deposition of an amount of identically-sized DNA molecules onto a polyacrylamide gel or other gel matrix is defined as an idealized unit impulse function, (xa). This function describes the properties of applying a unit amount of identically-sized DNA molecules onto the gel matrix at position xa (e.g. the top portion of the gel matrix) as unity but being zero for all other values of x

[0373] Therefore, at time t0, the origin of sample deposition at xa is characterized as an infinitely narrow band of identically-sized DNA molecules of unit area $\begin{array}{cc}{}_{-}^{+}\left(x-a\right)x=1.& \left(\mathrm{Equation}\mathrm{XXXVIII}\right)\end{array}$

[0374] Given these initial conditions, the fragment band will evolve by diffusion processes into a Gaussian shape.

[0375] Ordinary differentiation of Equation XXXV provides a solution for diffusion of DNA molecules in one dimension $\begin{array}{cc}C=\frac{\mathrm{const}}{\sqrt{t}}\exp \left(-x^2/4D_{p}t\right)& \left(\mathrm{Equation}\mathrm{XXXIX}\right)\end{array}$

[0376] where const is a constant related to the amount of DNA originally deposited onto the gel matrix. Equation XXXIX is symmetrical with respect to x0 and tends toward 0 as x either positively or negatively for t>0. For t0, x0 everywhere except at x0 where the function is described by the delta function of Equation XXXVIII. To evaluate the constant, const, consider the total amount of identically-sized DNA molecules deposited onto the gel matrix is defined as $\begin{array}{cc}M={}_{-}^{+}Cx.& \left(\mathrm{Equation}\mathrm{XL}\right)\end{array}$

[0377] Substituting x²/4D_Pt², dx{square root}{square root over (4D_Pt)}d provides a means for calculating M by $\begin{array}{cc}\begin{array}{c}M=\mathrm{const}\sqrt{4D_p}{}_{-}^{+}\exp \left(-{}^2\right)\\ =\mathrm{const}\sqrt{4D_p}.\end{array}& \left(\mathrm{Equation}\mathrm{XLI}\right)\end{array}$

[0378] Equation XLI shows that the amount M of DNA molecules diffusing remains constant and equal to the originally deposited amount at time t0 at position xa. In other words, Equation XXXIX normalizes to unit area (i.e. it applies to, for example, one DNA molecule per unit area) when Mconst{square root}{square root over (4D_P)}. Therefore, on substituting for const from Equation XLI in Equation XXXIX, the solution to the change in concentration in an amount M of identically-sized DNA fragments at a given time t and position x in one dimension is $\begin{array}{cc}C=\frac{M}{\sqrt{4D_{p}t}}{\exp }^{-\left(\frac{{\left(x-x_0\right)}^2}{4D_{p}t}\right)}& \left(\mathrm{Equation}\mathrm{XLII}\right)\end{array}$

[0379] and positions x, y in two dimensions is $\begin{array}{cc}C=\frac{M}{\sqrt{4D_{p}t}}{\exp }^{-\left(\frac{{\left(x-x_0\right)}^2}{4D_{p}t}-\frac{{\left(y-y_0\right)}^2}{4D_{p}t}\right)}.& \left(\mathrm{Equation}\mathrm{XLIII}\right)\end{array}$

[0380] Note that radial symmetry in the two-dimensional diffusion process is not assumed.

[0381] Diffusion of the DNA molecules occurs because they (1) drift under random kinetic motion and (2) more of them move from a high concentration region, located at the center of the amplicon, to a low concentration area that is located away from the center of the amplicon. The mean distance in which an individual DNA molecule diffuses is zero because diffusion in the x and x directions for the one-dimensional case is equally probable. However, the mean square distance is not zero: $\begin{array}{cc}{}_x^2={}_{-}^{+}{x}^2\frac{M}{\sqrt{4D_{p}t}}{\exp }^{\left(-\frac{{\left(x-x_0\right)}^2}{4D_{p}t}\right)}x& \left(\mathrm{Equation}\mathrm{XLIV}\right)\\ =2D_{p}t.& \left(\mathrm{Equation}\mathrm{XLV}\right)\end{array}$

[0382] Equation XLV shows that for passive diffusion, _x², which corresponds to the peak variance or width, is directly proportional to time, t, provided D_P remains constant. Similar arguments are also directly applicable to the two-dimensional case of diffusion.

[0383] Random events other than passive diffusion of DNA molecules also contribute substantially to the effective diffusion rate of a DNA molecule in an isotropic gel matrix. In gel electrophoresis systems, DNA molecules are also actively transported by an electrical field through a porous gel or polymer matrix. The molecular conductance pathways change randomly in response to random changes in the composition of the gel or polymer matrix; DNA molecules may be forced to stop and diffuse laterally before resuming forward migration. Hence, a molecule of DNA will undergo many changes in velocity, and thus position, thereby contributing to an increase in the effective diffusion rate. Considering only one-dimensional molecular displacement, incorporating the effects of electrically-mediated displacement velocities into Equation XXXV above leads to $\begin{array}{cc}\frac{C}{t}=-V^{}\left(\frac{C}{x}\right)+D_P\left(\frac{{}^{2}C}{x^2}\right)& \left(\mathrm{Equation}\mathrm{XLVI}\right)\end{array}$

[0384] where V is the displacement velocity of the DNA molecule in the gel matrix. Note that V will also experience random velocity fluctuations and therefore contribute an additional variance component to the effective diffusion rate. To account for these random velocity fluctuations, the displacement velocity V can be broken into two parts

[0385] where V is the average displacement velocity (assumed to be constant) and V_A incorporates the random component of velocity fluctuations. Therefore, $\begin{array}{cc}\frac{C}{t}=-V\left(\frac{C}{x}\right)-V_A\left(\frac{C}{x}\right)+D_P\left(\frac{{}^{2}C}{x^2}\right).& \left(\mathrm{Equation}\mathrm{XLVIII}\right)\end{array}$

[0386] The random nature of V_A thus contributes an additional diffusion term $\begin{array}{cc}-V_A\left(\frac{C}{x}\right)=D_A\left(\frac{{}^{2}C}{x^2}\right)& \left(\mathrm{Equation}\mathrm{XLIX}\right)\end{array}$

[0387] where D_A is a constant describing the diffusion of a DNA molecule attributable to active (e.g. electrically mediated) displacement. Substituting Equation XLVIII into Equation XLVII leads to $\begin{array}{cc}\frac{C}{t}=-V\left(\frac{C}{x}\right)+\left(D_A+D_P\right)\left(\frac{{}^{2}C}{x^2}\right).& \left(\mathrm{Equation}L\right)\end{array}$

[0388] Equation L is of the same form of Equations XXXV and XXXVI and yields similar Gaussian solutions as presented by Equations XLII and XLIII. The effective diffusion rate is enhanced because the coefficient ²C/x² of is now the sum of two diffusion terms

[0389] One term, D_P, represent passive diffusion of the DNA molecule and the other term, D_A, represents diffusion due to active (e.g. electrically mediated) displacement of the DNA molecule. Similar to Equation XLV, the total effective diffusion coefficient, D_Total, of each DNA amplicon acquires a variance in the form of

[0390] that is also proportional to time, t, provided that D_Total remains constant. For example, at constant V, the distance, x, traversed by the amplicon in time, t, is simply xVt. Substituting tx/V into Equation LII leads to $\begin{array}{cc}{}_x^2=\left(\frac{2D_{\mathrm{Total}}}{V}\right)x.& \left(\mathrm{Equation}\mathrm{LIII}\right)\end{array}$

[0391] Equation LIII shows that _x², based on the combined effects of passive and active diffusion, is also directly proportional to the amplicon migration distance, x. Identical arguments can easily be made for the same passive and active diffusion effects in two-dimensional systems.

[0392] It is a well-known that the separation of DNA molecules in a gel matrix occurs because of different migration velocities attributable to DNA molecules of distinct sizes. For example, small DNA molecules move through a gel matrix at greater velocities than do larger DNA molecules. This phenomenon creates a third source of apparent diffusion due to differences in the velocity by which a DNA molecule passes by the detection region of the detector. For example, many automated DNA sequencers (e.g. ABI capillary and slab-gel sequencers) operate in what is known as a finish-line mode where differently-sized DNA molecules migrate the same distance but at different times.

[0393] When applying finish-line electrophoresis, small-sized amplicons will move across the detector region relatively quickly while larger-sized amplicons will pass by the detector region more slowly. Consequently, due to the size-dependent migration velocities of DNA molecules, the apparent magnitude of diffusion will be larger for amplicons of larger size. In other words, the magnitude of _x² and _y², as predicted by the aforementioned one- and -two dimensional diffusion models, respectively, will typically be an underestimate of the observed values of _x² and _y² because of the time a particular DNA molecule takes to pass the detection region.

[0394] Provided the gel matrix and electrophoresis conditions remain constant, the observed time-delays for identically-sized DNA molecules should be identical or nearly identical across samples and across electrophoresis runs. However, in situations where this assumption cannot be made, it is desirable to compensate for this effect. A preferred approach is to convert the observed values of _x² and/or _y² on the original time-delayed scale to actual values of amplicon width (in units of length) by multiplying with the velocity of the migrating amplicon. For example, the velocity at which an amplicon moves through a gel matrix may be computed by the simple relationship, L/t, where is the velocity of the DNA molecule (e.g. amplicon), L is the distance from the sample-well (or injection position) to the center of the detection region, and t is the time it takes for the amplicon to reach the center of the detection region. In ABI-based detection systems, the distance from the sample-well to the detection window is known. The time it takes the amplicon to reach the detection window may also be determined by simply converting the unit of scan number to a unit of time.

[0395] It may be desirable to calculate the true diffusion coefficients, e.g. to compare the observed results to theoretical expectations. In a preferred implementation, the total variance in amplicon width is rewritten as a sum of the variance in amplicon width due to electrophoresis, _{EIectrophoresis}², and variance in amplicon width as a result of additional detector delay, _Detector², leading to ²_{ElectrophoresisDetector}². By assuming that the additional increase in observed amplicon width is entirely due to diffusion, the true diffusion coefficient is related to the additional amplicon width as _Detector²2D_Totalt_delay, where t_delay is the delay time of an amplicon in the detector region. Thus, the true diffusion coefficient may be estimated, for example, by plotting the amplicon-width variances as a function of the delay times; the diffusion coefficient is equal to 0.5 multiplied to the slope of a linear, least-squares fit to the amplicon widths.

[0396] The diffusion model described above provides a means by which to accurately model the distribution of DNA molecules comprising an amplicon in an electrophoretic system. In turn, the same model provides the mathematical criteria necessary to compute the total amount of spectral intensity in an amplicon from photometric data. This procedure not only provides a relative estimate of the number of DNA molecules in an amplicon, but also provides a means to infer the number of DNA templates from which an amplicon was originally derived. Together, these capabilities ultimately make available a means to assign genotypes to specific amplicons where none could be designated previously. A preferred implementation of a procedure for estimating spectral intensity 56 (FIG. 1) is described herein further detail.

[0397] Dividing Equation LII (and its two-dimensional form) into Equations XLII and XLIII, respectively, and replacing the amplitude variable M/{square root}{square root over (4D_Totalt)} with the variable I, solvable parametric models of DNA molecules comprising an amplicon, C(x), in an isotropic medium are obtained for the one-dimensional case $\begin{array}{cc}C\left(x\right)=I{\exp }^{-\left(\frac{{\left(x-x_0\right)}^2}{2{}_x^2}\right)}& \left(\mathrm{Equation}\mathrm{LIV}\right)\end{array}$

[0398] and for the two-dimensional case $\begin{array}{cc}C\left(x,y\right)=I{\exp }^{-\left(\frac{\left(x-x_0\right)}{2{}_x^2}-\frac{\left(y-y_0\right)}{2{}_y^2}\right).}& \left(\mathrm{Equation}\mathrm{LV}\right)\end{array}$

[0399] Equations LIV and LV may be fitted to photometric data preferably using nODR or nOLS procedures as described above. Preferably, estimates of peak amplitude, peak location, and amplicon variance that have been previously obtained with nODR or nOLS techniques or other methods described herein are used as starting values for the preferred nODR minimization procedure.

[0400] In the one-dimensional case, for example, the value of I corresponds to peak amplitude (e.g. in units of RFI), x₀ is the true peak location (e.g. in units of scan number), x is an observed location in the set of photometric data (e.g. in units of scan number), and _x² is the amplicon variance (e.g. in units of scan number, either corrected or uncorrected for detector time-delay effects). In a preferred implementation, it is desired to use estimates of these variables, namely parameters I₁, x₀₂, and _x²₃, previously generated using Equation XXV above. In that case, Equation LIV is rewritten to $\begin{array}{cc}C\left(x\right)={}_1{\exp }^{-\left(\frac{{\left(x-{}_2\right)}^2}{2{}_3^2}\right).}& \left(\mathrm{Equation}\mathrm{LVI}\right)\end{array}$

[0401] For the two-dimensional case, the definition and units of the amplitude I remain the same as in the one dimensional case. However, the definitions and units of x₀, x and _x² must be adjusted to reflect changes due to the addition of another dimension. In two dimensions, x₀ is the true peak location in the first dimension (e.g. in units of channel number), x is an observed location in the first dimension in the set of photometric data (e.g. in units of channel number), and _x² is the amplicon variance in the first dimension (e.g. in units of channel number). Note that there is no delay response associated with the detector region in the x dimension of channel number.

[0402] To mathematically accommodate the second dimension, three additional variables are necessary: y₀ is the true peak location (e.g. in units of scan number), y is an observed location in the set of photometric data (e.g. in units of scan number), and _y² is the amplicon variance (e.g. in units of scan number either corrected or uncorrected for detector time-delay effects). In a preferred implementation, it is desired to use estimates of parameters ₁, ₂, and ₃ estimated from Equation XXV as starting values for parameters I₁, y₀₄, and _x²₅in the two dimensional case, respectively. Note that the x dimension coordinates and parameters were simply transposed to the y dimension coordinates for parameters ₄ and ₅.

[0403] The remaining variables in the second dimension, x₀ and _x² are estimated by parameters ₂ and ₃, respectively. Estimates for starting values for these two additional parameters are obtained as follows: In a preferred implementation, ₂ preferably is initially estimated as the channel number nearest to the intersection of the tracking spline that bisects the amplicon being fitted. In a preferred implementation, ₅ is preferably estimated by directly measuring the width of the sample well, in units of channel numbers. For example, 64-lane combs (part no. CAG64-020) manufactured by The Gel Company (665 Third Street, Suite 240, San Francisco, Calif. 94107) have a well-width of approximately 4-5 channels. If other sample-well combs are used, and/or used on other detectors with differing numbers of channels, routine experimentation can be performed to determine the width of the sample wells in units of channels. Preferably, as discussed above, estimates of all parameters in either one or two dimensions are determined prior to their use in this part of the procedure as they are highly accurate starting values and will allow quick convergence of the fitting algorithm.

[0404] In a preferred implementation of the two-dimensional case, estimates of these variables, as reflected in parameters ₁ . . . ₅, are substituted into Equation LV above. In that case, Equation LV is rewritten as $\begin{array}{cc}C\left(x,y\right)={}_1{\exp }^{-\left(\frac{{\left(x-{}_2\right)}^2}{2{}_3^2}-\frac{{\left(y-{}_4\right)}^2}{2{}_5^2}\right).}& \left(\mathrm{Equation}\mathrm{LVII}\right)\end{array}$

[0405] Estimates of the parameters in either the one- or two-dimensional case are made by fitting Equations LVI or LVII to the processed photometric data, respectively, preferably using nonlinear ODR or nonlinear OLS procedures as described above.

[0406] To compute the total amount of spectral intensity present in an individual amplicon, Equation LIV is integrated from to yield an analytical solution for the one-dimensional case $\begin{array}{cc}{}_{-}^{+}C\left(x\right)x=I{}_{-}^{+}{\exp }^{-\left(\frac{{\left(x-x_0\right)}^2}{2{}_x^2}\right)}x& \left(\mathrm{Equation}\mathrm{LVIII}\right)\\ =I\sqrt{2{}_x^2}& \left(\mathrm{Equation}\mathrm{LIX}\right)\end{array}$

[0407] Substituting the estimated values for peak amplitude, {circumflex over ()}₁, and amplicon variance, {circumflex over ()}₃, into Equation LVIII yields the following $\begin{array}{cc}{}_{-}^{+}C\left(x\right)x={\hat{}}_1\sqrt{2{\hat{}}_3^2}.& \left(\mathrm{Equation}\mathrm{LX}\right)\end{array}$

[0408] Equation LX allows direct calculation of the spectral intensity of the entire amplicon as reflected by the parameters which minimize the error of the fitting procedure described above. Integration of Equation LV yields a similar analytical solution for the two-dimensional case $\begin{array}{cc}{}_{-}^{+}{}_{-}^{+}C\left(x,y\right)xy=I{}_{-}^{+}{}_{-}^{+}{\exp }^{-\left(\frac{\left(x-x_0\right)}{2{}_x^2}-\frac{\left(y-y_0\right)}{2{}_y^2}\right)}xy& \left(\mathrm{Equation}\mathrm{LXI}\right)\\ =I\sqrt{4{}^{2}2{}_x^{2}2{}_y^2}.& \left(\mathrm{Equation}\mathrm{LXII}\right)\end{array}$

[0409] Substituting the estimated values for peak amplitude, {circumflex over ()}₁, peak variance in the x dimension, {circumflex over ()}₃, and peak variance in the y dimension, {circumflex over ()}₅, into Equation LXI yields the following equation $\begin{array}{cc}{}_{-}^{+}{}_{-}^{+}C\left(x,y\right)xy={\hat{}}_1\sqrt{4{}^2{\hat{}}_3^2{\hat{}}_5^2}.& \left(\mathrm{Equation}\mathrm{LXIII}\right)\end{array}$

[0410] Similarly to the one-dimensional case, Equation LXIII allows direct calculation of the spectral intensity of the entire amplicon as reflected by the parameters which minimize the error of the fitting procedure described above.

[0411] Once spectral intensities have been estimated for all amplicons of interest, all elements in the putative monomorphic/polymorphic matrix where M_ij1 are substituted with the value of integrated spectral intensity corresponding to the mapped amplicons in matrix S_ij.

[0412] For the vast majority of amplicons, such as those depicted in FIG. 21A, the preceding implementation of spectral intensity estimation provides extremely accurate estimates of total spectral intensity. However, situations do arise, especially in two-dimensional electrophoretic systems, where DNA molecules do not always follow the aforementioned diffusion processes. Typically, this occurs where the amplicon is an anomalously-shaped amplicon, such as is depicted in FIG. 21B. This phenomenon most commonly results from: (1) deformations in the sample loading well, (2) deformation of the gel or polymer matrix at the point of initial fragment deposition, (3) anisotropy of the gel or polymer matrix, and/or (4) physical debris within the gel matrix including air bubbles. In such cases, the model based on characteristics of ssDNA diffusion in a gel matrix, as described above, does not always yield a good fit to the data. This is often observed as an uncharacteristically large residual sum-of-squares of the model fit. In such cases, it is desirable to estimate the spectral intensity of an amplicon using a less restrictive (i.e. more general) surface-fitting model. A preferred implementation of such a model is described below.

[0413] In another preferred procedure for estimating the spectral intensity of amplicons, a highly generalized surface-fitting procedure employing thin-plate splines (TPS) can be implemented to accommodate anomalously-shaped DNA fragments. This type of surface-fitting methodology has an important advantage over the Gaussian model described above in that it does not depend on any intrinsic characteristics of (1) the electrophoretic device, (2) the separation matrix, or (3) the general behavior of DNA molecules in the separation matrix. Therefore, the TPS surface-fitting procedure can accommodate not only patterns resulting from the diffusion processes described above, but also almost any possible deformation that may occur during electrophoresis.

[0414] Briefly, TPS surface-fitting procedures can be considered as generalizations of standard multivariate linear regression procedures, except that the parametric model of Gaussian diffusion described above is replaced by a non-parametric function that preferably is a spline. Such TPS functions have certain desirable features, such as control over how smooth the function will fit to the observed data. The degree of smoothness of the fitted TPS function is determined by minimizing a measure of the predictive error of the fitted surface. This measure is preferably determined by generalized cross-validation techniques or by routine experimentation. Once the TPS is fitted to the observed spectral intensity data, the total spectral intensity may be computed by applying standard integration procedures to the fitted surface. Application of TPS functions fitted to photometric data as used herein has consistently yielded smaller variances in estimated spectral intensity than the preferred one- and two-dimensional Gaussian procedures discussed above. A more detailed description of a preferred implementation of the TPS surface-fitting procedure is given below. Although we restrict our discussion of TPS fitting to two dimensions, its use in one dimension is also possible and preferred when the one-dimensional Gaussian procedure does not yield acceptable results.

[0415] A TPS function is a generalized smoothing function that interpolates the surface of regularly or irregularly spaced data of the form

[0416] Restricting the TPS function to two dimensions, consider a two-dimensional independent variable X(x₁, x₂) and assume observations are from the additive model

[0417] where Y_i is the observed response at the i^th combination variable x_i, is the function of interest, and n is the number of observations. In practice, the true function is unknown, but is assumed to be smooth. The random error components, _i, are assumed to be independent and normally distributed as N(0, _i²). The goal of applying a TPS is to build a surface function (x_i) whose graph passes near or through all the observed landmark points x₁, x₂ while minimizing the function $\begin{array}{cc}{J}_2\left(f\right)={}_{-}^{+}{}_{-}^{+}\left({\left(\left[\frac{{}^{2}f}{x_1^2}\right]\right)}^2+{2\left[\frac{{}^{2}f}{x_1{x}_2}\right]}^2+\left[\frac{{}^{2}f}{x_2^2}\right]\right)x_1{x}_2.& \left(\mathrm{Equation}\mathrm{LXVI}\right)\end{array}$

[0418] Equation LXVI is often referred to as the bending energy of the TPS function as it reflects the overall degree of curvature imposed by the fit of the function to the given surface. With reference to the accompanying invention, x_i may be interpreted as pixel coordinates x_i, y_i and (x_i) as relative fluorescence intensity z_i preferably in two-dimensional gel systems. For tabulated coordinates, x₁ , x₂, (x_i), describing a single amplicon derived from a two-dimensional gel system, the TPS minimizing function is of the form $\begin{array}{cc}{S}_{}\left(f\right)=\frac{1}{n}\underset{i=1}{\overset{n}{}}{\left(y_i-f\left(x_i\right)\right)}^2+J_2\left(f\right)& \left(\mathrm{Equation}\mathrm{LXVII}\right)\end{array}$

[0419] where >0 is a smoothing parameter describing the roughness (cf. precision) of the TPS surface fit. In other words, the parameter controls the balance between the goodness of fit and the smoothness of the final TPS surface-fitting function. For example, when 0, the estimate of (x_i) is equivalent to exact interpolation of the data when the graph of the TPS function passes through every point in (x_i). Given the presence of noise in the spectral measurements of fluorescence intensity, exact interpolation is usually not desired. Increasing the value of allows greater flexibility of the TPS fit and provides an effective means by which to accommodate spectral noise.

[0420] Choosing the magnitude of the smoothing parameter greatly affects the final character of the TPS surface fit and thus requires careful selection. An objective methodology to select relies on generalized cross-validation (GCV) procedures, which approximate the predicted mean square error of the TPS-fitted function. Let A() denote a nn smoothing matrix that maps the observed data vector y onto the predicted values of the TPS fit

[0421] where T indicates the transpose of the matrix. The trace of A() is interpreted as a measure of the number of effective parameters present in the TPS representation of the surface. The residuals of the TPS fit are given by (IA())y with ntraceA() reflecting the degrees of freedom associated with the residual error. The GCV function is thus defined as $\begin{array}{cc}V\left(\right)=\frac{\left(1/n\right)y^T(I-{A\left(\right)}^{T}W\left(I-A\left(\right)\right)y}{{\left(1-\mathrm{trace}\left(A\left(\right)\right)/n\right)}^2}.& \left(\mathrm{Equation}\mathrm{LXVIX}\right)\end{array}$

[0422] The value of that minimizes the GCV function yields the TPS model fit that is least affected by any single data point. Once {circumflex over ()} is found, the variance of the surface fit, ², may be estimated by $\begin{array}{cc}{\hat{}}_{\hat{}}^2=\frac{{y^T\left(I-A\left(\hat{}\right)\right)}^{T}W\left(I-A\left(\hat{}\right)\right)y}{n-\mathrm{trace}\left(A\left(\hat{}\right)\right)}& \left(\mathrm{Equation}\mathrm{LXX}\right)\end{array}$

[0423] and is analogous to the estimate of the residual variance from classical linear regression. The use of the GCV function to select is based on asymptotic theory. Therefore, it is important that sample sizes contain enough observations to provide a sufficient amount of information so as to be able to discriminate between actual signal and background noise. In practice, experience suggests that the signal-to-noise ratio of data actually obtained is more than adequate to provide reliable estimates of using the GCV procedure. Once the TPS function is fitted to the data, the spectral intensity is estimated using standard integration techniques.

[0424] An example of this preferred implementation of a TPS-fitted surface to an amplicon in two dimensions is shown in FIG. 23. A contour plot of spectral intensity for a distorted amplicon (also see FIG. 21B) is shown in FIG. 23A. A TPS-fitted surface was computed as described above, using a value of 0.0001, and is shown in FIG. 23B. Compared to the Gaussian fit of a non-distorted amplicon shown in FIG. 22B, the TPS-fitted surface is not symmetrical and conforms much more closely to the original distribution of spectral intensity exhibited by the amplicon. This enhanced fit is graphically depicted in contour plots shown in panels C and D of FIG. 23. Compared to the surface fitted by the Gaussian model, the estimates of the residual sum-of-squares, R_SS, and the standard deviation of the fit, SD, is significantly less in the TPS-fitted surface.

[0425] After calculating the total spectral intensity, the Gaussian-fitted surface yielded an estimate of approximately 6 percent less RFI than the TPS-fitted surface. While 6 percent difference in RFI may not seem too significant, it is nevertheless additional variability that contributes to spectral noise and may prevent the correct assignment of genotype. Moreover, the magnitude of error is expected to increase depending on (1) the resulting distortion of the amplicon and (2) the magnitude of RFI in the amplicon. Therefore, it is desirable to minimize surface-fitting errors whenever possible.

[0426] The preceding transformations applied to the original photometric data have attenuated most of the random variation in signal intensity attributed primarily to electrophoresis- and detector-derived noise. However, variation in (1) sample-loading volumes and (2) concentration differences among amplicons in individual marker assays (i.e. the actual number of DNA molecules comprising a particular amplicon) and other variability also require correction 58 (FIG. 1). Once this procedure is done, preferably by the use of an internal normalization standard, the observed values of spectral intensity for each amplicon may be compared directly to each other to ascertain whether the amplicon was originally derived from one or two DNA templates (alleles). Recall that in the case of where no amplicon is present (null), no alleles were originally present and thus no amplicon is detected. As this comparison is relative across amplicons at an individual locus, a standard by which to gauge differences in spectral intensity among amplicons is mandatory.

[0427] Where pedigree information is known and provides certainty in knowledge of one or more monomorphic loci, a known monomorphic locus can be used to normalize the spectral intensities of all other amplicons using the normalization procedure set forth below. Where a monomorphic locus is not known a priori, e.g. when no pedigree or other suitable generic information is available, identification of a monomorphic locus must be performed before the normalization process can proceed.

[0428] In cases where monomorphic loci cannot be identified a priori, a minimum-variance search strategy preferably is implemented among putatively monomorphic loci, pMM_k, k1, 2, . . . , m (i.e. those loci j in M_ij where for all i, M_ij1), to identify the locus with the least amount of variability in spectral intensity. The putatively monomorphic locus having the least amount of variability preferably is designated as a true monomorphic locus and used to normalize the spectral intensities of all other amplicons. The rationale behind this strategy is that if a locus is actually monomorphic, the variance in spectral intensity among amplicons at that locus should be smaller than if there was at least one amplicon with a polymorphic genotype present at that locus. Therefore, any locus j in M_ij having at least one null homozygote (i.e. those loci j in M_ij where for all i, at least one M_ij0) cannot be monomorphic and therefore are designated polymorphic, PM_l, l1, 2, . . . , p.

[0429] Having now omitted polymorphic sites from consideration of a candidate monomorphic locus, the selection of an actual monomorphic locus rests on the ability to distinguish among those amplicons having either homozygote present (preferably) or heterozygote genotypes. A difficulty arises in that it is the primary goal of this invention to distinguish among these genotypes, but to do so first requires the identification of a monomorphic locus. The circularity of this problem can be resolved by applying the rationale of the minimum-variance search strategy described above to each of the putatively monomorphic loci. Thus, it is expected that if a locus is actually monomorphic (i.e. comprised only of amplicons with preferably homozygote present genotypes), the variance in spectral intensity at this locus will be smaller than if there was at least one amplicon with a heterozygote genotype present at that locus.

[0430] In a preferred procedure for selecting a monomorphic locus, a coefficient of variation, , in spectral intensity is first computed for all amplicons at each putative monomorphic locus. The coefficient of variation for spectral intensity at the k^th putatively monomorphic locus pMM_k, k1, 2, . . . , m, is given by $\begin{array}{cc}{}_k=\frac{\sqrt{S_{{\mathrm{pMM}}_k}^2}}{\stackrel{\_}{{\mathrm{pMM}}_k}}100.& \left(\mathrm{Equation}\mathrm{LXXI}\right)\end{array}$

[0431] In Equation LXXI, S_pMMk² is the k^th observed sample variance of spectral intensity in amplicons across all S_i, i1, 2, . . ., N samples and {overscore (pMM)}_k is the k^th average observed spectral intensity in amplicons across all S_i, i1, 2, . . . , N samples. Computed values of for all k putatively monomorphic loci preferably are arranged in ascending order. Where this procedure is implemented in software, ordering all of the values of _k may not be required.

[0432] The locus k exhibiting the smallest value of _k is designated to be the monomorphic reference locus, rMM, provided that _k is below an empirically derived magnitude, typically a value of 12.000. Cross-correlation analyses of spectral intensities of amplicons at rMM to amplicons at loci with successively larger values of may be performed to gauge the reliability of the original reference locus as a useful monomorphic locus. If a monomorphic reference locus cannot be defined, either by a priori knowledge or by using this procedure, it may not be possible to perform genotyping to discriminate among homozygote present and heterozygote genotypes. Where this procedure is implemented in software and a reference monomorphic locus cannot be defined, a message preferably is generated.

[0433] If a monomorphic reference locus is identified, preferably, the average spectral intensity, {overscore (rMM)}, of all amplicons at the reference locus is computed. This value is used to compute a normalizing coefficient, _i, for the amplicon in the i^th sample at the reference locus rMM $\begin{array}{cc}{}_i=\frac{{\mathrm{rMM}}_i}{\mathrm{rMM}}& \left(\mathrm{Equation}\mathrm{LXXII}\right)\end{array}$

[0434] where i1, 2, . . . N. This procedure generates a set of normalizing coefficients which, when used to adjust the spectral intensity of corresponding amplicons at other loci, will remove the majority of variability in spectral intensity due to differences in (1) the individual marker assay, _{Marker assay}², and (2) the sample loading volume, _{Loading volume}². In a preferred implementation, each coefficient _i is divided into corresponding value of spectral intensity observed at other amplicon positions i for any other locus j. For example, the normalizing coefficients may be applied to any locus in pMM_k, k1, 2, . . . , m or PM_l, l1, 2, . . . p. The resulting normalized spectral intensity will more closely reflect values of expected spectral intensity if no variability in the marker assays or loading volumes existed. Thus, the normalizing coefficients, _i, computed from the i^th amplicon in rMM serve directly as the source of normalizing coefficients for all other loci in the entire sample. These normalizing coefficients form the basis by which genotypes are assigned to individual amplicons at other loci. This procedure is described in further detail below.

[0435] Once the spectral intensities for amplicons at other loci are normalized, as described above, the resulting numerical values are prepared for genotyping 60 (FIG. 1). On a per locus basis, the spectral intensities of amplicons are expected to vary, possibly considerably, due to differences in the relative efficiency of the PCR reaction. It is therefore desirable to perform an additional normalization step to transform the normalized intra-locus spectral intensities to a common inter-locus scale.

[0436] In a preferred implementation, the normalized values of spectral intensity are scaled to unity. For example, at a given locus, the amplicon exhibiting the maximum normalized spectral intensity is identified and divided into the spectral intensities of the remaining amplicons at that locus. Thus, the spectral intensity for each amplicon at that locus will be uniformly scaled to an RFI of 1.0. An advantage of scaling in this manner is that it allows simultaneous graphical visualization of any number of loci, e.g. FIG. 24, all of which are normalized to a common scale. Moreover, this procedure directly allows amplicon genotype assignments to be made based only on the relative proportion of normalized spectral intensity. FIG. 24 graphically illustrates the proportion of normalized spectral intensity computed for eleven loci using a preferred implementation of the normalization procedure described above. At all eleven loci, the three genotypic classes, homozygote present, heterozygote, and homozygote absent (null), are clearly separated from each other. Thus, visual inspection alone, without the reliance on additional statistical models, is sufficient enough to accurately assign genotypes (but see below).

[0437] In practice, it appears uncommon not to be able to unambiguously assign amplicons to a specific genotypic class using a simple graphical display such as illustrated in FIG. 24. In such cases where a genotype could not be unambiguously assigned using this procedure, it has been attributed to either (1) amplicon entry into the plateau-phase or (2) a poor fit of the preferred Gaussian model (due to amplicon distortion). In the latter case, refitting the photometric data of the amplicon in question with the TPS model and recalculating the spectral intensity resulted in an unambiguous assignment of genotype to that amplicon. In the former case, reducing the number of PCR cycles in the selective-phase of the fAFLP process resulted in an unambiguous assignment of genotype to those amplicons.

[0438] An additional advantage of this method is that it allows a user to gauge whether or not amplicons at a particular locus are entering the plateau phase. For example, referring to FIG. 24, the vast majority of heterozygous amplicons have not yet passed the inflection point (the proportion of normalized spectral intensity of 0.50 or greater) of the log-linear amplification phase illustrated in FIG. 3. Thus, presenting the data in this manner is useful in determining the optimum number of PCR cycles for a particular primer pair.

[0439] The invention also appears robust when inferring genotypes that are not derived exclusively from 0, 1 or 2 copies of an allele. In fact, the invention appears capable of distinguishing among any number of alleles. Such an instance has been observed in certain individuals for sex-linked loci in Domestic Chicken (Gallus gallus). For example, where the pedigree of parents and offspring is known, as was the case in a validation study, the origin and distribution of sex-chromosomes in offspring is predictable. At locus 559.562, we obtained genotypes in certain individuals consistent with 2, 3, and 4 allelic doses.

[0440] Nevertheless, an indeterminable genotype may arise using the preferred procedures of this invention. If it is determined that other identifiable processes are not responsible, genotype assignments may also be made statistically as briefly described below.

[0441] If desired, a statistically-based model can be used to optionally gauge or predict classification probabilities of genotypes assigned to specific classes. Many such models already exist. For example, both (1) genetic models based on expected allele distributions or (2) purely mathematical models, e.g., discriminant models, are suitable for gauging or predicting classification probabilities into a particular genotypic class. Thus, either aforementioned model can be used to statistically assign a genotype to an amplicon that could not be ascertained by visual inspection after normalized spectral intensities were scaled to unity and plotted in the manner depicted in FIG. 24. Either of these models can also assign a probability of correct assignment into a specific genotypic class for any spectral intensity representing a specific amplicon. Note that reference to the following two classification models is not inclusive of all available models to accomplish this task.

[0442] For example, the model discussed by Piepho and Koch (2000), Genetics, 155(3):1459-68, is a genetic model based on Bayesian expectations in normal mixtures of genotypes. Alternatively, a purely probabilistic discriminant function may also be used, e.g. the three-class discriminant function described by Fisher, 1938. Although both types of models are implemented in very different ways, each would be used to predictively classify the normalized spectral intensity of amplicons at loci into a defined number of classes, for example, homozygous present, heterozygous, and homozygous absent (null).

[0443] Nevertheless, the performance of either model to classify genotypes using spectral intensity data is almost completely unknown and has not been applied to data processed in the manner described herein. For example, genotype assignments presented in Piepho and Koch (2000) show a strong overlap of the probability distribution between the classes of bomozygous present and heterozygote, a phenomenon we generally do not observe. It is asserted that this is likely due to the failure to accommodate into their model important issues related to (1) spectral noise minimization, (2) spectral intensity estimation, and (3) spectral intensity normalization. Thus, it is believed that their method would be more accurate using data processed in accordance with methods of this invention.

[0444] The following Examples are provided for illustrative purposes only. The Examples are included herein solely to aid in a more complete understanding of the presently described invention. The Examples do not limit the scope of the invention described or claimed herein in any fashion.

[0445] Sample validation. To validate the FAFLP codominant scoring procedure, a system where precise knowledge of parental and offspring genotypes was used. Two highly inbred lines of domestic chicken (Gallus gallus), a single Rhode-Island Red male and 8 Ancona females were crossed. A total of 48 F1 individuals were obtained from six hens. Nine autosomal and two sex-linked loci were used to validate the codominant fAFLP procedure. DNA samples of parents P:M3966 and P:F343 and 10 F1 individuals (Fl:M/F7163-7173) were randomized blindly with replacement and processed using the fAFLP genotyping procedure, such as that preferably described above. The results of the genotyping procedure are graphically illustrated in FIG. 24.

[0446] fAFLP Analysis. For fAFLP analysis of avian genomic DNA (e.g. G. gallus and Manacus manacus) two to three selective bases were sufficient to produce on average, 120 resolvable fragments 60 to 600 bases in length. For example, in M. manacus, each of twenty primer pairs, applied to individuals sampled at a small lek (7 males, 5 females), always demonstrated significant levels of polymorphism. The percentage of polymorphic loci for each of twenty primer pairs ranged from 15 to 47 percent, with an average of 23 percent.

[0447] It is understood that the various preferred embodiments shown and described above illustrate different possible features of the invention and the varying ways in which these features may be combined. Apart from combining the different features of the above embodiments in varying ways, other modifications are also considered to be within the scope of the invention. For example, the method can be used to analyze AFLP-generated fragments labeled with fluorescent dyes other than those specifically mentioned and labeled with dyes other than fluorescent dyes, e.g., radiolabels. The method can also be used to analyze AFLP fragments generated from nucleic acids other than DNA. For example, the fragments can be generated from mRNA. As noted above, the amplification products can be analyzed with a variety of instruments, including, but not limited to, automated DNA sequencers (both slab-gel and capillary-based), and densitometers.