normal - Distribución de ajustes, bondad de ajuste, p-valor. ¿Es posible hacer esto con Scipy(Python)?

scipy stats norm 0 1

Question

INTRODUCCIÓN: Soy bioinformática. En mi análisis que realizo en todos los genes humanos (aproximadamente 20 000) busco un motivo particular de secuencia corta para verificar cuántas veces aparece este motivo en cada gen.

Los genes están ''escritos'' en una secuencia lineal en cuatro letras (A, T, G, C). Por ejemplo: CGTAGGGGGTTTAC ... Este es el alfabeto de cuatro letras del código genético que es como el lenguaje secreto de cada célula, es la forma en que el ADN almacena la información.

Sospecho que las repeticiones frecuentes de una secuencia de motivo corto particular (AGTGGAC) en algunos genes son cruciales en un proceso bioquímico específico en la célula. Dado que el motivo en sí es muy corto, es difícil con las herramientas computacionales distinguir entre los verdaderos ejemplos funcionales en los genes y aquellos que se ven similares por casualidad. Para evitar este problema, obtengo secuencias de todos los genes y las concatené en una sola cuerda y las barajé. La longitud de cada uno de los genes originales fue almacenada. Luego, para cada una de las longitudes de secuencia originales, se construyó una secuencia aleatoria seleccionando repetidamente A o T o G o C al azar de la secuencia concatenada y transfiriéndola a la secuencia aleatoria. De esta manera, el conjunto resultante de secuencias aleatorias tiene la misma distribución de longitud, así como la composición general A, T, G, C. Luego busco el motivo en estas secuencias aleatorias. Realicé este procedimiento 1000 veces y promedié los resultados.

15000 genes que no contienen un motivo dado 5000 genes que contienen 1 motivo 3000 genes que contienen 2 motivos 1000 genes que contienen 3 motivos ... 1 gen que contiene 6 motivos

Entonces, incluso después de 1000 veces la aleatorización del código genético verdadero, no hay genes que tengan más de 6 motivos. Pero en el verdadero código genético, hay algunos genes que contienen más de 20 apariciones del motivo, lo que sugiere que esta repetición podría ser funcional y es poco probable que los encuentre en tal abundancia por pura casualidad.

PROBLEMA: Me gustaría saber la probabilidad de encontrar un gen con, digamos, 20 ocurrencias del motivo en mi distribución. Así que quiero saber la probabilidad de encontrar tal gen por casualidad. Me gustaría implementar esto en Python, pero no sé cómo.

¿Puedo hacer tal análisis en Python?

Cualquier ayuda sería apreciada.

Answer 1

En la documentación de SciPy encontrará una lista de todas las funciones de distribución continua implementadas. Cada uno tiene un método de fit() , que devuelve los parámetros de forma correspondientes.

Incluso si no sabe qué distribución utilizar, puede probar muchas distrubuciones simultáneamente y elegir la que mejor se adapte a sus datos, como en el código a continuación. Tenga en cuenta que si no tiene idea de la distribución, puede ser difícil ajustar la muestra.

import matplotlib.pyplot as plt
import scipy
import scipy.stats
size = 20000
x = scipy.arange(size)
# creating the dummy sample (using beta distribution)
y = scipy.int_(scipy.round_(scipy.stats.beta.rvs(6,2,size=size)*47))
# creating the histogram
h = plt.hist(y, bins=range(48))

dist_names = [''alpha'', ''beta'', ''arcsine'',
              ''weibull_min'', ''weibull_max'', ''rayleigh'']

for dist_name in dist_names:
    dist = getattr(scipy.stats, dist_name)
    param = dist.fit(y)
    pdf_fitted = dist.pdf(x, *param[:-2], loc=param[-2], scale=param[-1]) * size
    plt.plot(pdf_fitted, label=dist_name)
    plt.xlim(0,47)
plt.legend(loc=''upper left'')
plt.show()

Referencias:

- Distribución ajustada con Scipy.

- ¿Ajustando la distribución empírica a las teóricas con Scipy (Python)?