Dopasowywanie dystrybucji, dobro dopasowania, wartość p. Czy można to zrobić za pomocą Scipy (Python)?

Question

WPROWADZENIE: Jestem bioinformatykiem. W mojej analizie, którą wykonuję na wszystkich ludzkich genach (około 20 000), szukam konkretnego motywu krótkiej sekwencji, aby sprawdzić ile razy ten motyw występuje w każdym genie.

Geny są "pisane" w liniowej sekwencji w czterech literach (A, T, G, C). Na przykład: CGTAGGGGGTTTAC ... Jest to czteroliterowy alfabet kodu genetycznego, który jest jak tajny język każdej komórki, tak naprawdę DNA przechowuje informacje.

Podejrzewam, że częste powtarzanie określonej sekwencji krótkiego motywu (AGTGGAC) w niektórych genach ma kluczowe znaczenie w konkretnym procesie biochemicznym w komórce. Ponieważ sam motyw jest bardzo krótki, trudno jest narzędziom obliczeniowym odróżnić prawdziwe funkcjonalne przykłady w genach od tych, które wyglądają podobnie przez przypadek. Aby uniknąć tego problemu, otrzymuję sekwencje wszystkich genów i łączę się w jeden ciąg i tasuje. Długość każdego z oryginalnych genów została zapisana. Następnie dla każdej z pierwotnych długości sekwencji skonstruowano losową sekwencję przez wielokrotne wybieranie A lub T lub G lub C losowo ze sprzężonej sekwencji i przeniesienie jej do losowej sekwencji. W ten sposób otrzymany zestaw losowych sekwencji ma taki sam rozkład długości, jak również ogólny skład A, T, G, C. Następnie szukam motywu w tych losowych sekwencjach. Wykonałem tę procedurę 1000 razy i uśredniłem wyniki.

15000 geny nie zawierać określony motyw 5000 geny, które zawierają 1 motif 3000 geny, które zawierają 2 motywy 1000 genów, które zawierają 3 motywy ... jeden gen, który zawiera 6 motywy

Dlatego nawet po 1000 razy randomizacji prawdziwego kodu genetycznego nie ma żadnych genów, które mają więcej niż 6 motywów. Ale w prawdziwym kodzie genetycznym jest kilka genów, które zawierają więcej niż 20 wystąpień tego motywu, co sugeruje, że te powtórzenia mogą być funkcjonalne i jest mało prawdopodobne, aby znalazły je w takiej obfitości przez czysty przypadek.

PROBLEM: Chciałbym poznać prawdopodobieństwo znalezienia genu z powiedzmy 20 wystąpień motywu w mojej dystrybucji. Chcę więc poznać prawdopodobieństwo znalezienia takiego genu przez przypadek. Chciałbym wdrożyć to w Pythonie, ale nie wiem jak.

Czy mogę przeprowadzić taką analizę w Pythonie?

Każda pomoc zostanie doceniona.

Answer 1

In SciPy documentation znajdziesz listę wszystkich zaimplementowanych funkcji dystrybucji ciągłej. Każdy ma a fit() method, który zwraca odpowiednie parametry kształtu.

Nawet jeśli nie wiesz, której dystrybucji użyć, możesz wypróbować wiele rozwiązań jednocześnie i wybrać ten, który pasuje lepiej do danych, jak w poniższym kodzie. Zwróć uwagę, że jeśli nie masz pojęcia o rozkładzie, może być trudno dopasować próbkę.

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() 

Referencje:

- Distribution fitting with Scipy

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