Plotting DNA chromatogram trace data

Sanger DNA sequencing produces a chromatogram trace which can be visualized with a number of programs, including FinchTV or ChromasLite. The raw data consists of co-ordinates for each of the four DNA bases (A, C, G and T). However, plots are displayed as smoothed peaks, as shown in the image previously. Many programs allow the x- and y-axis of the plot to be scaled larger or smaller to change the shape of the plot. What mathematical method is used to plot smooth curves such as these from a (small) number of raw data points?

• convolution? What are the small number of raw data points? What does the smoothness tell you? I'm confused, I though the smooth signal was what came out of the sequencer, and you're interested in the correct letter that corresponds to. – Henry Gomersall Jul 9 '12 at 16:10
• @HenryGomersall: Yes and no :) What comes out of the sequencer are points which, when plotted, lie on that nice smoothed line. I want to plot the smoothed line. The letter corresponding to the peak is a different, but related, problem as I understand it. – SabreWolfy Jul 9 '12 at 17:48
• Oh, so your question is about interpolation for a pretty rendering? To be clear, this isn't about processing the data per se? – Henry Gomersall Jul 9 '12 at 18:30
• @HenryGomersall: Ultimately I wanted to do some processing of the data, but I wanted to start by figuring out how to plot the curves. – SabreWolfy Jul 9 '12 at 21:33

Implementation

Assuming you already have a line drawing routine, you just need to supplement that with some kind of interpolation. The curves are created by drawing enough interpolated short lines to make the result look smooth. A good starting point would be to use an existing interpolation routine, like the ones given by Paul Bourke here.

I'll illustrate this using the cubic routines he supplies, since those are some of the simplest that will still give reasonable results. Here's the first one (translated to python) for reference:

def cubic(mu,y0,y1,y2,y3):
mu2 = mu*mu
a0 = y3 - y2 - y0 + y1
a1 = y0 - y1 - a0
a2 = y2 - y0
a3 = y1
return a0*mu*mu2 + a1*mu2 + a2*mu + a3

Each routine has a parameter mu which represents the fractional part of the index you want to interpolate. Depending on the routine, the other parameters will be some number of samples surrounding the index in question. In the cubic case, you'll need four samples. For example, if your data is y[n], and you want the value at 10.3, mu would be .3, and you'd pass in y[9], y[10], y[11], and y[12].

Then instead of drawing a single line with endpoints, say, $(10,y_{10}) \rightarrow (11,y_{11})$, you'd draw a bunch of shorter ones using the interpolated values (e.g. $(10,y_{10}) \rightarrow (10.1,\text{cubic}(.1,y_9,y_{10},y_{11},y_{12})) \dots$). Obviously those points would need to be scaled to the $x$ and $y$ dimensions of the image to be rendered.

Theory

Now, since the page/routine I referenced doesn't cite any sources, it's worth explaining where those cubic routines come from (and how they work). Both the one I reproduced above, and the very similar Catmull-Rom spline he mentions just below it, are two specific cases of using the following cubic convolution kernel:

$$\psi(x) = \begin{cases} (\alpha + 2)|x|^3 - (\alpha + 3)|x|^2 + 1, & \text{ if } 0 \le |x| \lt 1 \\ \alpha|x|^3 - 5\alpha|x|^2 + 8\alpha|x| - 4\alpha, & \text{ if } 1 \le |x| \lt 2 \\ 0, & \text{ if } 2 \le |x| \end{cases}$$

The routine listed above corresponds to a value of $\alpha = -1$, and the Catmull-Rom spline corresponds to $\alpha = -1/2$. I won't go into too much detail about how the general form of the kernel is derived, but it involves various constraints like making sure that $\psi(x)$ is one at zero, and zero at all other integers.

This is what it looks like:

The two choices for the value of $\alpha$ come from attempts to match various aspects of the sinc function, the ideal reconstruction kernel. Setting $\alpha = -1$ makes the derivative of $\psi$ match the derivative of the sinc function at $x = 1$, and making it equal to $-1/2$ provides the best low-frequency approximation instead. By all accounts, a value of $\alpha = -1/2$ has much better properties overall, so it's probably the best value to use in practice. A much more extensive discussion can be found in the following paper, starting on page 328:

Meijering, Erik. “A Chronology of Interpolation: From Ancient Astronomy to Modern Signal and Image Processing.” Proceedings of the IEEE. vol. 90, no. 3, pp. 319-42. March 2002.

Insight

Now, just looking at that kernel in relation to the actual implementation of the interpolation code, it may not be clear how the two are related. Basically, the interpolation process can be thought of as adding together shifted copies of the kernel, which are scaled by the samples of the data, like so:

In fact, if you have an implementation of the kernel, you can use it directly to do the interpolation, as follows:

def kernel(x, a=-1.0):
x = abs(x)
if x >= 0.0 and x < 1.0:
return (a + 2.0)*x**3.0 - (a + 3.0)*x**2.0 + 1
elif x >= 1.0 and x < 2.0:
return a*x**3.0 - 5.0*a*x**2.0 + 8.0*a*x - 4.0*a
else:
return 0.0

def cubic(mu,y0,y1,y2,y3):
a = -1.0
result  = y0 * kernel(mu + 1, a)
result += y1 * kernel(mu, a)
result += y2 * kernel(mu - 1, a)
result += y3 * kernel(mu - 2, a)
return result

However, it's much less computationally efficient to do it that way. As a bridge from the direct kernel approach to the more streamlined one above, consider that with a little bit of algebraic manipulation, the first implementation can be put in the following form:

def cubic(mu,y0,y1,y2,y3):
mu2 = mu*mu
mu3 = mu*mu2
c0 = -mu3 + 2*mu2 - mu
c1 =  mu3 - 2*mu2 + 1
c2 = -mu3 + mu2 + mu
c3 =  mu3 - mu2
return c0*y0 + c1*y1 + c2*y2 + c3*y3

In this formulation, the c0...c3 values can be thought of as the coefficients of an FIR filter which is applied to the sample values. Now it's much easier to see how to derive the routine from the kernel. Consider the kernel with $\alpha = -1$, like so:

$$\psi(x) = \begin{cases} |x|^3 - 2|x|^2 + 1, & \text{ if } 0 \le |x| \lt 1 \\ -|x|^3 + 5|x|^2 - 8|x| + 4, & \text{ if } 1 \le |x| \lt 2 \\ 0, & \text{ if } 2 \le |x| \end{cases}$$

Now evaluate that kernel symbolically at various shifted offsets, keeping in mind that mu ($\mu$) ranges from 0 to 1:

$$\begin{eqnarray} \psi(\mu + 1) &=& -(\mu + 1)^3 + 5(\mu + 1)^2 - 8(\mu + 1) + 4 &=& −\mu^3 + 2\mu^2 - \mu & & (c_0) \\ \psi(\mu) &=& \mu^3 - 2\mu^2 + 1 &=& \mu^3 - 2\mu^2 + 1 & & (c_1)\\ \psi(\mu - 1) &=& (1 - \mu)^3 - 2(1 - \mu)^2 + 1 &=& -\mu^3 + \mu^2 + \mu & & (c_2)\\ \psi(\mu - 2) &=& -(2 - \mu)^3 + 5(2 - \mu)^2 - 8(2 - \mu) + 4 &=& \mu^3 - \mu^2 & & (c_3)\\ \end{eqnarray}$$

Note that $\mu - 1, \mu - 2$ get "flipped" to $1 - \mu, 2 - \mu$ respectively because of the absolute value on $x$ in the kernel definition. Now we have the exact polynomials that are used in the "FIR version" of the interpolation routine. The evaluation of those polynomials can then be made more efficient through standard techniques (e.g. Horner's method). Similar things can be done with other kernels, and there are other ways of constructing efficient implementations as well (cf. Digital Audio Resampling Home Page).