Это задание поможет вам лучше разобраться в методах множественной проверки гипотез и позволит применить ваши знания на данных из реального биологического исследования.
В этом задании вы:
вспомните, что такое t-критерий Стьюдента и для чего он применяется сможете применить технику множественной проверки гипотез и увидеть собственными глазами, как она работает на реальных данных почувствуете разницу в результатах применения различных методов поправки на множественную проверку Основные библиотеки и используемые методы:
Библиотека scipy и основные статистические функции:http://docs.scipy.org/doc/scipy/reference/stats.html#statistical-functions
Библиотека statmodels для методов коррекции при множественном сравнении:
http://statsmodels.sourceforge.net/devel/stats.html
Статья, в которой рассматриваются примеры использования statsmodels для множественной проверки гипотез:
http://jpktd.blogspot.ru/2013/04/multiple-testing-p-value-corrections-in.html
Описание используемых данных
Данные для этой задачи взяты из исследования, проведенного в Stanford School of Medicine. В исследовании была предпринята попытка выявить набор генов, которые позволили бы более точно диагностировать возникновение рака груди на самых ранних стадиях.
В эксперименте принимали участие 24 человек, у которых не было рака груди (normal), 25 человек, у которых это заболевание было диагностировано на ранней стадии (early neoplasia), и 23 человека с сильно выраженными симптомами (cancer).
Ученые провели секвенирование биологического материала испытуемых, чтобы понять, какие из этих генов наиболее активны в клетках больных людей.
Секвенирование — это определение степени активности генов в анализируемом образце с помощью подсчёта количества соответствующей каждому гену РНК.
В данных для этого задания вы найдете именно эту количественную меру активности каждого из 15748 генов у каждого из 72 человек, принимавших участие в эксперименте.
Вам нужно будет определить те гены, активность которых у людей в разных стадиях заболевания отличается статистически значимо.
Кроме того, вам нужно будет оценить не только статистическую, но и практическую значимость этих результатов, которая часто используется в подобных исследованиях.
Диагноз человека содержится в столбце под названием "Diagnosis".
Практическая значимость изменения
Цель исследований — найти гены, средняя экспрессия которых отличается не только статистически значимо, но и достаточно сильно. В экспрессионных исследованиях для этого часто используется метрика, которая называется fold change (кратность изменения). Определяется она следующим образом:
Fc(C,T)= T/C, if T > C; -C/T, if C > T.
где C,T — средние значения экспрессии гена в control и treatment группах соответственно. По сути, fold change показывает, во сколько раз отличаются средние двух выборок.
Инструкции к решению задачи
Задание состоит из трёх частей. Если не сказано обратное, то уровень значимости нужно принять равным 0.05.
Часть 1: применение t-критерия Стьюдента
В первой части вам нужно будет применить критерий Стьюдента для проверки гипотезы о равенстве средних в двух независимых выборках. Применить критерий для каждого гена нужно будет дважды:
для групп normal (control) и early neoplasia (treatment) для групп early neoplasia (control) и cancer (treatment) В качестве ответа в этой части задания необходимо указать количество статистически значимых отличий, которые вы нашли с помощью t-критерия Стьюдента, то есть число генов, у которых p-value этого теста оказался меньше, чем уровень значимости.
Часть 2: поправка методом Холма
Для этой части задания вам понадобится модуль multitest из statsmodels.
1 import statsmodels.stats.multitest as smm В этой части задания нужно будет применить поправку Холма для получившихся двух наборов достигаемых уровней значимости из предыдущей части. Обратите внимание, что поскольку вы будете делать поправку для каждого из двух наборов p-value отдельно, то проблема, связанная с множественной проверкой останется.
Для того, чтобы ее устранить, достаточно воспользоваться поправкой Бонферрони, то есть использовать уровень значимости 0.05 / 2 вместо 0.05 для дальнейшего уточнения значений p-value c помощью метода Холма.
В качестве ответа к этому заданию требуется ввести количество значимых отличий в каждой группе после того, как произведена коррекция Холма-Бонферрони. Причем это число нужно ввести с учетом практической значимости: посчитайте для каждого значимого изменения fold change и выпишите в ответ число таких значимых изменений, абсолютное значение fold change которых больше, чем 1.5.
Обратите внимание, что
применять поправку на множественную проверку нужно ко всем значениям достигаемых уровней значимости, а не только для тех, которые меньше значения уровня доверия. при использовании поправки на уровне значимости 0.025 меняются значения достигаемого уровня значимости, но не меняется значение уровня доверия (то есть для отбора значимых изменений скорректированные значения уровня значимости нужно сравнивать с порогом 0.025, а не 0.05)! Часть 3: поправка методом Бенджамини-Хохберга
Данная часть задания аналогична второй части за исключением того, что нужно будет использовать метод Бенджамини-Хохберга.
Обратите внимание, что методы коррекции, которые контролируют FDR, допускает больше ошибок первого рода и имеют большую мощность, чем методы, контролирующие FWER. Большая мощность означает, что эти методы будут совершать меньше ошибок второго рода (то есть будут лучше улавливать отклонения от H0, когда они есть, и будут чаще отклонять H0, когда отличий нет).
В качестве ответа к этому заданию требуется ввести количество значимых отличий в каждой группе после того, как произведена коррекция Бенджамини-Хохберга, причем так же, как и во второй части, считать только такие отличия, у которых abs(fold change) > 1.5.
In [23]:
import numpy as np
import pandas as pd
import statsmodels
import scipy as sc
import statsmodels.stats.multitest as smm
In [3]:
data = pd.read_csv('gene_high_throughput_sequencing.csv')
In [4]:
data.head()
Out[4]:
In [6]:
data.shape
Out[6]:
In [8]:
data['Diagnosis'].value_counts()
Out[8]:
In [9]:
normal = data[data['Diagnosis'] == 'normal']
In [11]:
neoplasia = data[data['Diagnosis'] == 'early neoplasia']
In [13]:
cancer = data[data['Diagnosis'] == 'cancer']
In [33]:
def get_fc(c, t):
if t >= c:
return t / c
else:
return -c / t
In [34]:
n1, n2 = 0, 0
p1 = np.zeros(data.shape[1] - 2)
p2 = np.zeros(data.shape[1] - 2)
fc1 = np.zeros(data.shape[1] - 2)
fc2 = np.zeros(data.shape[1] - 2)
for i in xrange(2, data.shape[1]):
gene_normal = normal.iloc[:, i]
gene_neoplasia = neoplasia.iloc[:, i]
gene_cancer = cancer.iloc[:, i]
p1[i-2] = sc.stats.ttest_ind(gene_normal, gene_neoplasia, equal_var = False).pvalue
p2[i-2] = sc.stats.ttest_ind(gene_neoplasia, gene_cancer, equal_var = False).pvalue
fc1[i-2] = get_fc(gene_normal.mean(), gene_neoplasia.mean())
fc2[i-2] = get_fc(gene_neoplasia.mean(), gene_cancer.mean())
n1 = len(p1[p1[:] < 0.05])
n2 = len(p2[p2[:] < 0.05])
print n1, n2
with open('task1.1.txt', 'w') as f:
f.write(str(n1))
with open('task1.2.txt', 'w') as f:
f.write(str(n2))
In [39]:
len(fc1[fc1[:] > 1.5]), len(fc2[fc2[:] > 1.5])
Out[39]:
In [43]:
n1, n2 = 0, 0
reject, p1_corr, a1, a2 = smm.multipletests(p1, alpha = 0.025, method = 'holm')
reject, p2_corr, a1, a2 = smm.multipletests(p2, alpha = 0.025, method = 'holm')
for i in xrange(len(p1_corr)):
if p1_corr[i] < 0.025 and np.abs(fc1[i]) > 1.5:
n1 += 1
if p2_corr[i] < 0.025 and np.abs(fc2[i]) > 1.5:
n2 += 1
print n1, n2
with open('task2.1.txt', 'w') as f:
f.write(str(n1))
with open('task2.2.txt', 'w') as f:
f.write(str(n2))
In [44]:
n1, n2 = 0, 0
reject, p1_corr, a1, a2 = smm.multipletests(p1, alpha = 0.025, method = 'fdr_bh')
reject, p2_corr, a1, a2 = smm.multipletests(p2, alpha = 0.025, method = 'fdr_bh')
for i in xrange(len(p1_corr)):
if p1_corr[i] < 0.025 and np.abs(fc1[i]) > 1.5:
n1 += 1
if p2_corr[i] < 0.025 and np.abs(fc2[i]) > 1.5:
n2 += 1
print n1, n2
with open('task3.1.txt', 'w') as f:
f.write(str(n1))
with open('task3.2.txt', 'w') as f:
f.write(str(n2))
In [ ]: