Фитинги распределения, качество посадки, р-значение. Возможно ли это сделать с помощью Scipy (Python)?

Question

ВВЕДЕНИЕ: Я биоинформатик. В своем анализе, который я выполняю по всем генам человека (около 20 000), я ищу конкретный мотив короткой последовательности, чтобы проверить, сколько раз этот мотив встречается в каждом гене.

Гены «записаны» в линейной последовательности из четырех букв (A, T, G, C). Например: CGTAGGGGGTTTAC ... Это четырехбуквенный алфавит генетического кода, похожий на секретный язык каждой клетки, именно так ДНК на самом деле хранит информацию.

Я подозреваю, что частые повторения определенной короткой мотивной последовательности (AGTGGAC) в некоторых генах имеют решающее значение в специфическом биохимическом процессе в клетке. Поскольку сам мотив очень короткий, с помощью вычислительных инструментов трудно различить истинные функциональные примеры в генах и те, которые выглядят случайно. Чтобы избежать этой проблемы, я получаю последовательности всех генов и объединяю их в одну строку и перемешиваю. Длина каждого из исходных генов была сохранена. Затем для каждой из исходных длин последовательностей была построена случайная последовательность путем многократного выбора A или T или G или C случайным образом из сцепленной последовательности и переноса ее в случайную последовательность. Таким образом, результирующий набор рандомизированных последовательностей имеет такое же распределение длины, как и общий состав A, T, G, C. Затем я ищу мотив в этих рандомизированных последовательностях. Я выполнил эту процедуру 1000 раз и усреднил результаты.

15000 генов, которые не содержат данный мотив 5000 генов, которые содержат 1 мотив 3000 генов, которые содержат 2 мотива 1000 генов, которые содержат 3 мотива ... 1 ген, содержащий 6 мотивов

Таким образом, даже после 1000-кратной рандомизации истинного генетического кода нет никаких генов, которые бы имели более 6 мотивов. Но в истинном генетическом коде есть несколько генов, которые содержат более 20 экземпляров мотива, что говорит о том, что эти повторы могут быть функциональными, и вряд ли их можно найти в таком изобилии по чистой случайности.

ПРОБЛЕМА: Я хотел бы знать вероятность нахождения гена с, скажем, 20 появлений мотива в моем распределении. Поэтому я хочу знать вероятность случайного нахождения такого гена. Я хотел бы реализовать это в Python, но я не знаю, как.

Могу ли я сделать такой анализ в Python?

Любая помощь будет оценена.

Answer 1

В документации SciPy вы найдете список всех реализованных функций непрерывного распределения.Каждый из них имеет a fit() метод , который возвращает соответствующие параметры формы.

Даже если вы не знаете, какой дистрибутив использовать, вы можете попробовать одновременно несколько диструбуций и выбрать тот, которыйлучше подходит для ваших данных, как в коде ниже.Обратите внимание, что если вы не имеете представления о распределении, может быть трудно подобрать образец.

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

Ссылки:

- Распределениеподгонка с Scipy

- Подгонка эмпирического распределения к теоретическим с помощью Scipy (Python)?