Распределение распределения, добротность посадки, значение p. Можно ли это сделать с помощью Scipy (Python)?

ВВЕДЕНИЕ: Я биоинформатик. В моем анализе, который я выполняю на всех генах человека (около 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?

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

Ответ 1

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

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

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