Distribuciones de montaje, bondad de ajuste, valor de p. ¿Es posible hacer esto con Scipy (Python)?

Question

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

Los genes se '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 cómo el ADN en realidad almacena información.

Sospecho que repetations frecuentes de una secuencia motivo corto en particular (AGTGGAC) en algunos genes que 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 ejemplos funcionales verdaderos en genes y aquellos que parecen similares por casualidad. Para evitar este problema, obtengo secuencias de todos los genes y las concateno en una sola cuerda y las mezclo. La longitud de cada uno de los genes originales se almacenó. 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 forma, el conjunto resultante de secuencias aleatorizadas tiene la misma distribución de longitud, así como la composición general de A, T, G, C. Luego busco el motivo en estas secuencias aleatorias. Performed este procedimiento 1000 veces y promedié los resultados.

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

Así incluso después de 1000 veces la aleatorización del verdadero código genético, no hay genes que tengan más de 6 motivos. Pero en el código genético verdadero, 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. Entonces quiero saber la probabilidad de encontrar tal gen por casualidad. Me gustaría implementar esto en Python, pero no sé cómo.

¿Puedo hacer un análisis de este tipo en Python?

Cualquier ayuda sería apreciada.

Answer 1

In SciPy documentation encontrará una lista de todas las funciones de distribución continuas implementadas. Cada uno tiene a fit() method, que devuelve los parámetros de forma correspondientes.

Incluso si usted no sabe qué distribución a utilizar puede probar muchos distrubutions simultáneamente y elegir el que se adapte mejor a sus datos, al igual que en el código de abajo. Tenga en cuenta que si no tiene idea acerca 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:

- Distribution fitting with Scipy

- Fitting empirical distribution to theoretical ones with Scipy (Python)?