Traçando dados de rastreamento do cromatograma de DNA

O sequenciamento de DNA da Sanger produz um traço de cromatograma que pode ser visualizado com vários programas, incluindo FinchTV ou ChromasLite. Os dados brutos consistem em coordenadas para cada uma das quatro bases de DNA (A, C, G e T). No entanto, os gráficos são exibidos como picos suavizados, conforme mostrado na imagem anteriormente. Muitos programas permitem que os eixos x e y da plotagem sejam redimensionados em tamanho maior ou menor para alterar a forma da plotagem. Que método matemático é usado para plotar curvas suaves como essas de um número (pequeno) de pontos de dados brutos?

Implementação

Supondo que você já tenha uma rotina de desenho de linha, basta suplementar isso com algum tipo de interpolação. As curvas são criadas desenhando linhas curtas interpoladas suficientes para tornar o resultado mais suave. Um bom ponto de partida seria usar uma rotina de interpolação existente, como as fornecidas por Paul Bourke aqui .

Ilustrarei isso usando as rotinas cúbicas que ele fornece, pois essas são algumas das mais simples que ainda fornecerão resultados razoáveis. Aqui está o primeiro (traduzido para python) para referência:

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

Cada rotina possui um parâmetro muque representa a parte fracionária do índice que você deseja interpolar. Dependendo da rotina, os outros parâmetros serão um número de amostras em torno do índice em questão. No caso cúbico, você precisará de quatro amostras. Por exemplo, se seus dados estão y[n], e você quer o valor em 10.3, museria .3, e você passar em y[9], y[10], y[11], e y[12].

Em vez de desenhar uma única linha com pontos finais, digamos, $(10,y_{10}) \rightarrow (11,y_{11})$, você desenharia vários mais curtos usando os valores interpolados (por exemplo, $(10,y_{10}) \rightarrow (10.1,\text{cubic}(.1,y_9,y_{10},y_{11},y_{12})) \dots$) Obviamente, esses pontos precisariam ser escalados para o$x$ e $y$ dimensões da imagem a ser renderizada.

Teoria

Agora, como a página / rotina que referi não cita nenhuma fonte, vale a pena explicar de onde vêm essas rotinas cúbicas (e como elas funcionam). Tanto o que reproduzi acima, quanto o spline Catmull-Rom muito semelhante que ele menciona logo abaixo, são dois casos específicos de uso do seguinte núcleo de convolução cúbica:

$$\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}$$

A rotina listada acima corresponde a um valor de $\alpha = -1$e o spline Catmull-Rom corresponde a $\alpha = -1/2$. Não entrarei em muitos detalhes sobre como a forma geral do kernel é derivada, mas envolve várias restrições, como garantir que$\psi(x)$ é um em zero e zero em todos os outros números inteiros.

Isto é o que parece:

As duas opções para o valor de $\alpha$provêm de tentativas de combinar vários aspectos da função sinc , o núcleo de reconstrução ideal. Configuração$\alpha = -1$ faz a derivada de $\psi$ coincidir com a derivada da função sinc em $x = 1$e tornando-o igual a $-1/2$fornece a melhor aproximação de baixa frequência. Em todas as contas, um valor de$\alpha = -1/2$possui propriedades muito melhores no geral, portanto é provavelmente o melhor valor para usar na prática. Uma discussão muito mais extensa pode ser encontrada no documento a seguir, começando na página 328:

Meijering, Erik. "Uma cronologia da interpolação: da astronomia antiga ao processamento moderno de sinais e imagens". Anais do IEEE. vol. 90, n. 3, pp. 319-42. Março de 2002.

Discernimento

Agora, apenas olhando para o kernel em relação à implementação real do código de interpolação, pode não estar claro como os dois estão relacionados. Basicamente, o processo de interpolação pode ser considerado como a adição de cópias deslocadas do kernel, que são dimensionadas pelas amostras dos dados, da seguinte forma:

De fato, se você tiver uma implementação do kernel, poderá usá-lo diretamente para fazer a interpolação, da seguinte maneira:

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

No entanto, é muito menos eficiente computacionalmente fazê-lo dessa maneira. Como uma ponte da abordagem direta do kernel para a mais simplificada acima, considere que, com um pouco de manipulação algébrica, a primeira implementação pode ser colocada da seguinte forma:

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

Nesta formulação, os c0...c3valores podem ser considerados os coeficientes de um filtro FIR que é aplicado aos valores da amostra. Agora é muito mais fácil ver como derivar a rotina do kernel. Considere o kernel com$\alpha = -1$, igual a:

$$\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}$$

Agora avalie esse kernel simbolicamente em várias compensações deslocadas, tendo em mente que mu($\mu$) varia de 0a 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}$$

Observe que $\mu - 1, \mu - 2$ seja "virado" para $1 - \mu, 2 - \mu$ respectivamente, devido ao valor absoluto no $x$na definição do kernel. Agora, temos os polinômios exatos que são usados ​​na "versão FIR" da rotina de interpolação. A avaliação desses polinômios pode então ser mais eficiente através de técnicas padrão (por exemplo , método de Horner ). Coisas semelhantes podem ser feitas com outros kernels e também existem outras maneiras de construir implementações eficientes ( consulte a Home Page da Digital Audio Resampling Home ).