Encuadre de la distribución binomial negativa para la secuenciación del ADN

Question

La distribución binomial negativa se ha convertido en un modelo popular para los datos de recuento (en concreto, el número esperado de lecturas de secuenciación dentro de una región determinada del genoma a partir de un experimento determinado) en bioinformática. Las explicaciones varían:

• Algunos la explican como algo que funciona como la distribución de Poisson de Poisson, pero con un parámetro adicional, lo que permite modelar la verdadera distribución, con una varianza no necesariamente igual a la media
• Algunos la explican como una mezcla ponderada de distribuciones de Poisson (con una distribución de mezcla gamma en el parámetro de Poisson)

¿Existe alguna forma de conciliar estos razonamientos con la tradicional definición de una distribución binomial negativa como modelo del número de éxitos de los ensayos Bernoulli antes de ver un cierto número de fracasos? ¿O debería pensar que es una feliz coincidencia que un mezcla ponderada de distribuciones Poisson con una distribución de mezcla gamma tenga la misma función de masa de probabilidad que la distribución binomial negativa?

Answer 1

IMOH, realmente creo que la distribución binomial negativa se utiliza por conveniencia.

Así que en RNA Seq hay una suposición común de que si se toma un número infinito de medidas del mismo gen en un número infinito de réplicas entonces la distribución verdadera sería lognormal. Esta distribución se muestrea a través de un proceso de Poisson (con un recuento), por lo que la distribución real de lecturas por gen a través de réplicas sería una distribución Poisson-Lognormal.

Pero en paquetes que utilizamos como EdgeR y DESeq esta distribución se modela como una distribución binomial negativa. Esto no es porque los chicos que lo escribieron no sabían acerca de una distribución de Poisson Lognormal.

El motivo es que la distribución Lognormal de Poisson es terrible para trabajar con ella porque requiere integración numérica para realizar los ajustes, etc., de modo que cuando intentas utilizarla, a veces el rendimiento es realmente malo.

Una distribución binomial negativa tiene una forma cerrada, por lo que es mucho más fácil trabajar con ella, y la distribución gamma (la distribución subyacente) se parece mucho a una distribución lognormal, en el sentido de que a veces parece normal y a veces tiene cola.

Pero en este ejemplo (si te crees el supuesto) es imposible que sea teóricamente correcta porque la distribución teóricamente correcta es la lognormal de Poisson y las dos distribuciones son aproximaciones razonables entre sí pero no son equivalentes.

Pero sigo pensando que la distribución binomial negativa "incorrecta" es a menudo la mejor opción porque empíricamente dará mejores resultados, ya que la integración funciona lentamente y los ajustes pueden funcionar mal, especialmente con distribuciones con colas largas.

Answer 2

He buscado en algunas páginas web y no he encontrado una explicación, pero se me ha ocurrido una para valores enteros de $r$ . Supongamos que tenemos dos fuentes radiactivas que generan independientemente partículas alfa y beta a las velocidades $\alpha$ y $\beta$ respectivamente.

¿Cuál es la distribución del número de partículas alfa antes del $r$ ¿la partícula beta?

1. Considera las partículas alfa como éxitos y las beta como fracasos. Cuando se detecta una partícula, la probabilidad de que sea una partícula alfa es $\frac{\alpha}{\alpha+\beta}$ . Entonces, esta es la distribución binomial negativa $\text{NB}(r,\frac{\alpha}{\alpha+\beta})$ .

2. Tenga en cuenta el tiempo $t_r$ de la $r$ ª partícula beta. Esta sigue una distribución gamma $\Gamma(r,1/\beta).$ Si condiciona en $t_r = \lambda/\alpha$ entonces el número de partículas alfa antes del tiempo $t_r$ sigue una distribución de Poisson $\text{Pois}(\lambda).$ Por lo tanto, la distribución del número de partículas alfa antes del $r$ a partícula beta es una distribución Gamma-mixta de Poisson.

Eso explica por qué estas distribuciones son iguales.

Answer 3

• Algunos la explican como algo que funciona como la distribución de Poisson pero tiene un parámetro adicional, lo que permite más libertad para modelar la distribución verdadera, con una varianza no necesariamente igual a la media
• Algunos lo explican como una mezcla ponderada de distribuciones de Poisson (con una distribución de mezcla gamma en el parámetro de Poisson)

Matemáticamente se obtiene la binomial negativa integrando la distribución de Poisson sobre pesos distribuidos Gamma, véase Mezcla Gamma-Poisson . Este hecho matemático se mantiene independientemente de que lo aceptemos o no como justificación para utilizar la distribución.

La distribución de Poisson es una elección bastante natural cuando se habla de contar las lecturas derivadas de la secuenciación del ADN (se podría utilizar la binomial, pero dado que sólo se secuencia una pequeña fracción de las lecturas/ADN obtenidas de la muestra, la diferencia es insignificante, y podemos utilizar la que nos parezca más conveniente). También estamos seguros de que el parámetro de esta distribución de Poisson varía, aunque la razón de esta variación depende de la naturaleza exacta del experimento - por ejemplo, puede ser una variación debida a

• repetir varias veces el mismo experimento
• las lecturas procedentes de células diferentes con propiedades algo distintas
• comparar el número de lecturas correspondientes a diferentes genes
• genes que tienen una estructura química diferente y, por tanto, se amplifican de forma diferente por PCR o algunas lecturas tienen más probabilidades de llegar a la máquina de secuenciación
• el estudiante/postdoc que prepara las bibliotecas no es muy cuidadoso/consistente
• etc.

En otras palabras, estamos seguros de que la variación existe (y la observamos experimentalmente), pero no sabemos exactamente de dónde procede, y no podemos saber directamente qué distribución de probabilidad la describe. No podríamos modelizarla utilizando la distribución normal, ya que el parámetro de Poisson debería ser positivo, así que utilizamos la distribución Gamma, porque es "casi como la normal", pero con soporte no negativo... pero también podríamos haber utilizado la log-normal u otra cosa. Mientras no estemos buscando los efectos biológicos finos que podrían resultar ser artefactos de la distribución particular que usamos, cualquier cosa que sea computacionalmente conveniente es buena.

Observe que, además de la flexibilidad que proporciona un parámetro adicional, la binomial negativa tiene una cola más gruesa que la distribución de Poisson, lo que la hace menos sensible a los valores atípicos. Esto proporciona una motivación adicional para utilizar esta distribución: permite inferencia más sólida .

Answer 4

Sólo puedo ofrecer una intuición, pero la propia distribución gamma describe tiempos de espera (continuos) (cuánto tarda en producirse un suceso raro). Así que el hecho de que una mezcla distribuida gamma de distribuciones poisson discretas daría lugar a un tiempo de espera discreto (ensayos hasta N fallos) no parece demasiado sorprendente. Espero que alguien tenga una respuesta más formal.

Edición: Siempre he justificado la dist. binomial negativa para la secuenciación de la siguiente manera: El paso real de secuenciación es simplemente muestreo de lecturas de una gran biblioteca de moléculas (poisson). Sin embargo esa biblioteca se hace a partir de la muestra original por PCR. Esto significa que las moléculas originales se amplifican exponencialmente. Y la distribución gamma describe la suma de k variables aleatorias independientes distribuidas exponencialmente, es decir, cuántas moléculas hay en la biblioteca después de amplificar k moléculas de la muestra durante el mismo número de ciclos de PCR.

De ahí los modelos binomiales negativos PCR seguidos de secuenciación.

Answer 5

Intentaré dar una interpretación mecanicista simplista que me ha resultado útil a la hora de pensar en esto.

Supongamos que tenemos una cobertura uniforme perfecta del genoma antes de la preparación de la biblioteca, y observamos $\mu$ lecturas que cubren un sitio por término medio. Digamos que la secuenciación es un proceso que recoge un fragmento original de ADN, lo somete a un proceso estocástico que pasa por PCR, submuestreo, etc., y obtiene una base del fragmento con una frecuencia $p$ y un fracaso por lo demás. Si la secuenciación continúa hasta $\mu\frac{1-p}{p}$ fallos, puede modelarse con una binomial negativa binomial negativa, $NB(\mu\frac{1-p}{p}, p)$ .

Calculando los momentos de esta distribución, obtenemos el número esperado de aciertos $\mu\frac{1-p}{p}\frac{p}{1-p} = \mu$ según sea necesario. Para la varianza del número de aciertos, obtenemos $\sigma^2 = \mu(1-p)^{-1}$ - la frecuencia con la que falla la preparación de la biblioteca para un fragmento aumenta la varianza en la cobertura observada.

Si bien lo anterior es una descripción un poco artificial del proceso de secuenciación, y se podría hacer un modelo generativo adecuado de los pasos de la PCR, etc, creo que da una idea del origen del parámetro de sobredispersión $(1-p)^{-1}$ directamente de la distribución binomial negativa. Prefiero el modelo de Poisson con la tasa integrada como explicación en general.