The present teachings pertain to systems and methods of aligning and analyzing point sets.
The skilled artisan will understand that the drawings, described below, are for illustration purposes only. The drawings are not intended to limit the scope of the present teachings in any way.
Protein expression studies can be invaluable in understanding basic biology and disease. Information gained from protein expression compares favorably to that obtained from mRNA expression. However, the technology behind protein expression analysis is less mature than that used for mRNA analysis, and many technical difficulties exist.
Recently developed liquid chromatography and mass spectrometry (LC/MS) protocols permit the study of high-throughput protein expression. In such protocols, protein mixtures are often digested into peptides and run through one or more chromatography stages. The output from the final chromatographic column can be sampled regularly, typically over a period of one to four hours, and the samples analyzed by a mass spectrometer. Each peptide produces a characteristic signal in the spectrometer, and these signals can be identified and gathered into features with a representative mass/charge (m/z) ratio, column retention time (rt), and signal intensity. A well-designed experiment will generally result in each peptide possessing a unique (m/z, rt) pair.
Early LC/MS protocols attempted to identify each peptide as it appeared in a mass spectrometer via tandem mass spectrometry (MS/MS). In many applications it can be more efficient to run many samples, identify peaks corresponding to peptides of interest, and later determine their identities. This requires the ability to compare the output of multiple LC/MS experiments, often called LC/MS maps, and identify which features correspond to the same peptide using only the (m/z, rt) pair information. This process is similar to comparing and associating spots on different two-dimensional gel images and is often referred to as the feature identification procedure. Once this procedure is complete, the data can be transformed into a matrix format where, for example, the rows correspond to peptides, columns correspond to experiments and matrix entries are the peptide intensities. Accuracy in peptide feature identification is often important to the quality of further analysis.
Two-dimensional gel electrophoresis (2DE) is another technology for conducting protein expression studies. For a comparison of 2DE and LC/MS see the following reference. (S. P. Gygi, G. L. Corthals, Y. Zhang, Y. Rochon, and R. Aebersold. Evaluation of the two-dimensional gel electrophoresis-based proteome analysis technology. Proceedings of the National Academy of Sciences, 97(17):9390-9395, 2000.) Despite differences in experimental technique, the problem of comparing results is similar to the LC/MS case. In 2DE, proteins form a set of spots on a gel. Software can be used to identify and quantify the spots, giving each spot a characteristic mass, isoelectric point, and intensity. These are analogous to the m/z, rt, and intensity values found in LC/MS maps. Typical software packages perform pair-wise analysis of gel images by comparing a new gel to a reference or a set of database gels. Pair-wise feature map comparison can be extended to a multi-feature map comparison by wrapping the pair-wise comparison into a hierarchical method like single-link clustering. (R. Sibson. SLINK: an optimally efficent algorithm for a complete link method. The Computer Journal, 16:30-34, 1973.) However, greedy hierarchical methods often perform poorly in the presence of significant noise. This is often due to an inherent difficulty in avoiding incorrect matchings at the pair-wise level.
While pair-wise comparison methods could be applied to LC/MS data, the present teachings provide several advantages over such methods. For example, the present teachings do not require a reference map. As well, the present teachings work well in noisy environments. LC/MS has the potential to be a high-throughput approach and a typical large-scale LC/MS experiment might involve running dozens of healthy samples, each with several replicates, and as many diseased samples, each with several replicates. The replicates can be important because LC/MS spectra are typically noisy and a large percentage of false features often appear in each map—potentially more false features than real ones. This can occur particularly when preprocessing software is tuned to be sensitive to low-abundance peptides. From this information, the researcher wishes to identify which features occur in which samples, and in particular, which are differentially expressed between the healthy and diseased classes.
The present teachings are useful for analyzing many LC/MS maps simultaneously, with or without a reference, and in the presence of significant amounts of noise.
Some embodiments of the present teachings associate with each feature a rectangular area of uncertainty. These rectangles can be selected so that features are more likely to be identical if their rectangles overlap. The rectangles can be constructed with prior knowledge of the measurement error, which can be obtained via simple experiments or estimated. The rectangles from all of the LC/MS maps are overlaid, and sets of mutually intersecting rectangles are identified. These represent the sets of features that originate from the same peptide; however, rectangles may appear in multiple sets. A heuristic can be used to assign each rectangle to the correct set. Since the heuristic has simultaneous access to all of the sets and their contents, the potential for more accurate results as compared to a greedy pair-wise merging strategy exists.
Prior to comparison, LC/MS maps are acquired and processed. Processing, is generally handled by signal processing software associated with the instrument. The present teachings can be used with a variety of LC/MS machines such as, but not limited to, the LCQ LC-MS from ThermoFinnigan, 1100 Series Trap/SL LC-MS from Agilent, MicroQuattro LC-MS from Waters or the API 2000™ LC/MS/MS System from Applied Biosystems. The instrument identifies features and provides the mass/charge ratio, elution time, intensity, and charge state of each feature. Generally, the features of interest are mass/charge ratio and elution time. In general, the accuracy of mass/charge ratios is higher than elution time, which can have significant variability.
In some embodiments, the input is a set of points (x1, y1), (x2, y2), . . . , (xN, yN) where the x-dimension represents the mass/charge measurements, and the y-dimension the elution time. The points correspond to the set of features identified across all maps. For each feature (xi, yj), εx(xi)=(xil,xiu) represents the uncertainty interval surrounding xi with lower bound xil and upper bound xiu. For example, if xi is the observed mass/charge ratio of a peptide in an LC/MS experiment, then the actual mass/charge ratio of the peptide lies within the interval (xil, xiu). Likewise, εy(yi) can represent the uncertainty interval around yi and the set (εx(xi), εy(yi)) represents a rectangular region, Ri, in the plane. The set of rectangles derived from the input points can be represented by R={R1, R2, . . . , RN}. Two features in the input set can be said to be compatible if their corresponding rectangles intersect. This indicates that both features might be observations of the same peptide. One stage of the feature identification procedure is to identify maximal sets (cliques) of mutually compatible features.
Some embodiments construct rectangular overlap graphs. In such embodiments, R can represent a set of iso-oriented rectangles in the plane, and G can represent an undirected graph such that there is a vertex corresponding to each rectangle in R and an edge between every pair of overlapping rectangles. Such a graph is called a rectangle overlap graph and is illustrated in
To implement such a method some embodiments employ an extension of the sweepline approach for finding maximal cliques among a set of intervals on a line. In the sweepline approach, the interval end points are the events, which the sweepline must process. In the sweepline method, a priority queue can be used to keep track of events to be processed by the sweepline. The queue tracks all of the intervals crossing the sweepline, and sorts them by their termination point. First, the left-hand end points of the intervals are sorted and the sweepline is placed at the leftmost point. At each step, the next event to be processed is either the next point from the list of left-hand end points (called a
The two-dimensional case is very similar.
Emission of sub-cliques of previously emitted cliques can be prevented via a simple observation. In
In some embodiments, the rectangles are projected to the x-axis, yielding a set of intervals. A sweepline processes these intervals as described above. However, the sweepline structure is itself a set of intervals in the y-dimension. When a left endpoint is encountered in the x-axis, the interval of the associated rectangle in the y-dimension is placed into the sweepline. When a right endpoint is encountered for the rectangle R, then the one-dimensional algorithm is run upon the intervals in the sweepline, and those cliques containing R are considered. Each clique is compared against the list of previously-emitted centroids as described above. If no centroids are contained within the clique's area of intersection, the clique is output and its centroid is added to the set of centroids.
A practical consideration is the assignment of rectangles to multiple cliques. If the input maps are sufficiently dense, this can be problematic. As each clique ideally represents a distinct peptide, resolving the assignment of ambiguous rectangles to cliques is important for the final result. Some embodiments assign rectangles to cliques using global context information thus performing a global optimization. Some embodiments employ randomized assignment, models of clique size distributions, models of intensity and charge distributions, or models of noise distribution. One skilled in the art will appreciate that a variety of assignment methodologies can be fashioned. Some embodiments employ a heuristic patterned on the expectation-maximization algorithm. For example, if area(R1, R2) denotes the area of intersection between rectangles R1 and R2, a score, s(R, C), for rectangle R with respect to clique C can be computed as follows.
The rectangle is then assigned to the clique that yields the greatest score. This scoring scheme generally avoids assigning rectangles to singleton cliques, and otherwise assigns them to cliques with which the average fit is best. Initially, each ambiguous rectangle is assigned to a clique. Some embodiments assign it to the clique with the most members. Then, the assignment procedure is iterated. During each iteration, each ambiguous rectangle R is moved to the clique C that maximizes s(R, C). Iteration can be terminated when there is no change or a maximum number of iterations is reached. Implementation of the scoring mechanism is illustrated in
Computer Implemented System
Computer system 500 may be coupled via bus 502 to a display 512, such as a cathode ray tube (CRT) or liquid crystal display (LCD), for displaying information to a computer user. An input device 514, including alphanumeric and other keys, is coupled to bus 502 for communicating information and command selections to processor 504. Another type of user input device is cursor control 516, such as a mouse, a trackball or cursor direction keys for communicating direction information and command selections to processor 504 and for controlling cursor movement on display 512. This input device typically has two degrees of freedom in two axes, a first axis (e.g., x) and a second axis (e.g., y), that allows the device to specify positions in a plane.
A computer system 500 performs the alignment method. Consistent with certain implementations of the invention, matched data points are provided by computer system 500 in response to processor 504 executing one or more sequences of one or more instructions contained in memory 506. Such instructions may be read into memory 506 from another computer-readable medium, such as storage device 510. Execution of the sequences of instructions contained in memory 506 causes processor 504 to perform the process described herein. Alternatively hard-wired circuitry may be used in place of or in combination with software instructions to implement the invention. Thus implementations of the invention are not limited to any specific combination of hardware circuitry and software.
The term “computer-readable medium” as used herein refers to any media that participates in providing instructions to processor 504 for execution. Such a medium may take many forms, including but not limited to, non-volatile media, volatile media, and transmission media. Non-volatile media includes, for example, optical or magnetic disks, such as storage device 510. Volatile media includes dynamic memory, such as memory 506. Transmission media includes coaxial cables, copper wire, and fiber optics, including the wires that comprise bus 502. Transmission media can also take the form of acoustic or light waves, such as those generated during radio-wave and infra-red data communications.
Common forms of computer-readable media include, for example, a floppy disk, a flexible disk, hard disk, magnetic tape, or any other magnetic medium, a CD-ROM, any other optical medium, punch cards, papertape, any other physical medium with patterns of holes, a RAM, PROM, and EPROM, a FLASH-EPROM, any other memory chip or cartridge, a carrier wave as described hereinafter, or any other medium from which a computer can read.
Various forms of computer readable media may be involved in carrying one or more sequences of one or more instructions to processor 504 for execution. For example, the instructions may initially be carried on the magnetic disk of a remote computer. The remote computer can load the instructions into its dynamic memory and send the instructions over a telephone line using a modem. A modem local to computer system 500 can receive the data on the telephone line and use an infra-red transmitter to convert the data to an infra-red signal. An infra-red detector coupled to bus 502 can receive the data carried in the infra-red signal and place the data on bus 502. Bus 502 carries the data to memory 506, from which processor 504 retrieves and executes the instructions. The instructions received by memory 506 may optionally be stored on storage device 510 either before or after execution by processor 504.
The foregoing description of an implementation of the invention has been presented for purposes of illustration and description. It is not exhaustive and does not limit the invention to the precise form disclosed. Modifications and variations are possible in light of the above teachings or may be acquired from practicing of the invention. Additionally, the described implementation includes software but the present invention may be implemented as a combination of hardware and software or in hardware alone. The invention may be implemented with both object-oriented and non-object-oriented programming systems.
Some of the embodiments herein are implemented in a software package called Rectangle Aggregator (RAG). RAG's performance has been compared to a program called the Aggregator that was developed to solve the feature identification problem using a hierarchical merging method. Results from RAG and the Aggregator on both simulated and real data are presented.
Simulated data has the advantage of being controllable and results from it are often more easily evaluated than when using real data. To test RAG, a large number of data sets were created. For each data set, a reference map was created with a number of features. The reference map was used only to assess the results. Each input map was populated with one observation per reference feature, normally distributed around the reference feature's position in both dimensions (m/z, rt). A number of uniformly distributed noise points were added to each map. The number of true features, noise features, and maps were varied to create many data scenarios. Ten replicates were performed for maps with 500 reference features; five were performed for maps with 2000 reference features. In each replicate, the input maps were provided to RAG and the Aggregator; neither program made use of the reference map. Each program was given a maximum matching distance in the mass/charge and retention time dimensions based on the parameters of the simulation.
The table shows that RAG and the Aggregator perform similarly in the ideal case where there is no noise, or where only two maps are being compared. However, as the amount of noise and number of maps grow, RAG is more accurate and the computational time is much lower. Extremely large test cases were not run with Aggregator due to its resource requirements.
RAG was tested on a set of 107 purified human serum albumin (HSA) samples. These samples were digested with bovine trypsin and sent through the LC/MS process. This resulted in 107 maps possessing between 504 and 715 features. To evaluate the results, several computational tryptic digestions of HSA and bovine trypsin—which remains in the sample after digestion—were performed. Only peptides likely to fall within the detectable mass/charge ratio range of the instrument were considered. If both ends of the peptides are constrained to be tryptic and no mis-cleavages occur, the result is 81 peptides. The number grows to 166 peptides with 1 mis-cleavage and 225 with 2. Of course, many variables affect how many features are detected; for example, a peptide may ionize to multiple charge states or, conversely, may fail to ionize at all. Nevertheless, the number of features called per map far exceeds the number expected; most of these extra features are much more likely to be noise than signal.
RAG detected 55 fully-populated groups across the samples, and 176 groups present in 80% of the groups. These numbers agree well with the theoretical digestions. In addition, a hand-curated set of 15 reliable, well-isolated reference features was developed in order to judge the quality of each individual map. RAG was able to successfully generate a clique for each reference feature and populate it with the correct representative from each map.
One skilled in the art will appreciate that the teachings herein are applicable to almost any protocol which produces sets of 1D or 2D points with error as output. Other potential applications include 1D and 2D gels, and single mass spectra, and the finding of biomarkers in serum which involves comparing single spectra from a large number of samples. (E. Petricoin, A. Ardekani, B. Hitt, P. Levine, V. Fusaro, S. Steinberg, G. Mills, C. Simone, D. Fishman, E. Kohn, and L. Liotta. Use of proteomic patterns in serum to identify ovarian cancer. Lancet, 359:572-577, 2002.)
Some embodiments of the teachings analyze samples in order to identify unknown peptides. In such cases, a reference sample containing known peptides is first analyzed producing a 2-dimensional feature map containing m/z and retention time information. Then, a plurality of samples, which can contain unknown peptides, are analyzed producing a plurality of 2-dimensional feature maps. Using embodiments of the teachings disclosed herein, the maps are aligned via global optimization. Unknown peptides that align with the known peptides contained in the reference map can be identified. One skilled in the art will appreciate that in addition to identification of the peptides, this technique can also be used to perform relative quantification by computing a ratio of the unknown to known peptides. One will also appreciate that absolute quantification can be reported if one of the amounts in the ratio is known.
This application claims a priority benefit under 35 U.S.C. § 119(e) from U.S. patent application Ser. No. 60/540,817, filed Jan. 30, 2004, which is incorporated herein by reference.
Number | Date | Country | |
---|---|---|---|
60540817 | Jan 2004 | US |