Распределение распределения, добротность посадки, значение 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, который возвращает соответствующие параметры формы.

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

enter image description here

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