Это задание поможет вам лучше разобраться в методах множественной проверки гипотез и позволит применить ваши знания на данных из реального биологического исследования.
вспомните, что такое 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 (кратность изменения). Определяется она следующим образом:
$$F_{c}(C,T) = \begin{cases} \frac{T}{C}, T>C \\ -\frac{C}{T}, T<C \end{cases}$$где C,T — средние значения экспрессии гена в control и treatment группах соответственно. По сути, fold change показывает, во сколько раз отличаются средние двух выборок.
Задание состоит из трёх частей. Если не сказано обратное, то уровень значимости нужно принять равным 0.05.
В первой части вам нужно будет применить критерий Стьюдента для проверки гипотезы о равенстве средних в двух независимых выборках. Применить критерий для каждого гена нужно будет дважды:
В качестве ответа в этой части задания необходимо указать количество статистически значимых отличий, которые вы нашли с помощью t-критерия Стьюдента, то есть число генов, у которых p-value этого теста оказался меньше, чем уровень значимости.
In [31]:
import pandas as pd
import scipy.stats
In [32]:
df = pd.read_csv("gene_high_throughput_sequencing.csv")
control_df = df[df.Diagnosis == 'normal']
neoplasia_df = df[df.Diagnosis == 'early neoplasia']
cancer_df = df[df.Diagnosis == 'cancer']
In [33]:
# scipy.stats.ttest_ind(data.Placebo, data.Methylphenidate, equal_var = False)
In [34]:
genes = filter(lambda x: x not in ['Patient_id', 'Diagnosis'], df.columns.tolist())
In [35]:
control_vs_neoplasia = {}
neoplasia_vs_cancer = {}
In [36]:
for gene in genes:
control_vs_neoplasia[gene] = scipy.stats.ttest_ind(control_df[gene], neoplasia_df[gene], equal_var = False).pvalue
neoplasia_vs_cancer[gene] = scipy.stats.ttest_ind(cancer_df[gene], neoplasia_df[gene], equal_var = False).pvalue
In [37]:
print control_df['LOC643837'],neoplasia_df['LOC643837']
In [38]:
scipy.stats.ttest_ind(control_df['LOC643837'], neoplasia_df['LOC643837'], equal_var = False).pvalue
Out[38]:
In [39]:
control_vs_neoplasia_df = pd.DataFrame.from_dict(control_vs_neoplasia, orient = 'index')
control_vs_neoplasia_df.columns = ['control_vs_neoplasia_pvalue']
In [40]:
neoplasia_vs_cancer_df = pd.DataFrame.from_dict(neoplasia_vs_cancer, orient = 'index')
neoplasia_vs_cancer_df.columns = ['neoplasia_vs_cancer_pvalue']
In [42]:
neoplasia_vs_cancer_df
Out[42]:
In [43]:
pvalue_df = control_vs_neoplasia_df.join(neoplasia_vs_cancer_df)
pvalue_df.head()
Out[43]:
In [44]:
pvalue_df[pvalue_df.control_vs_neoplasia_pvalue < 0.05].shape
Out[44]:
In [45]:
pvalue_df[pvalue_df.neoplasia_vs_cancer_pvalue < 0.05].shape
Out[45]:
Для этой части задания вам понадобится модуль multitest из statsmodels.
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)!
In [46]:
import statsmodels.stats.multitest as smm
In [47]:
pvalue_df['control_mean_expression'] = control_df.mean()
pvalue_df['neoplasia_mean_expression'] = neoplasia_df.mean()
pvalue_df['cancer_mean_expression'] = cancer_df.mean()
In [48]:
def abs_fold_change(c, t):
if t > c:
return t/c
else:
return c/t
In [49]:
pvalue_df['control_vs_neoplasia_fold_change'] = map(lambda x, y: abs_fold_change(x, y),
pvalue_df.control_mean_expression,
pvalue_df.neoplasia_mean_expression
)
pvalue_df['neoplasia_vs_cancer_fold_change'] = map(lambda x, y: abs_fold_change(x, y),
pvalue_df.neoplasia_mean_expression,
pvalue_df.cancer_mean_expression
)
In [50]:
pvalue_df['control_vs_neoplasia_rej_hb'] = smm.multipletests(pvalue_df.control_vs_neoplasia_pvalue, alpha=0.025, method='h')[0]
pvalue_df['neoplasia_vs_cancer_rej_hb'] = smm.multipletests(pvalue_df.neoplasia_vs_cancer_pvalue, alpha=0.025, method='h')[0]
In [51]:
pvalue_df[(pvalue_df.control_vs_neoplasia_rej_hb) & (pvalue_df.control_vs_neoplasia_fold_change > 1.5)].shape
Out[51]:
In [52]:
pvalue_df[(pvalue_df.neoplasia_vs_cancer_rej_hb) & (pvalue_df.neoplasia_vs_cancer_fold_change > 1.5)].shape
Out[52]:
Данная часть задания аналогична второй части за исключением того, что нужно будет использовать метод Бенджамини-Хохберга.
Обратите внимание, что методы коррекции, которые контролируют FDR, допускает больше ошибок первого рода и имеют большую мощность, чем методы, контролирующие FWER. Большая мощность означает, что эти методы будут совершать меньше ошибок второго рода (то есть будут лучше улавливать отклонения от H0, когда они есть, и будут чаще отклонять H0, когда отличий нет).
В качестве ответа к этому заданию требуется ввести количество значимых отличий в каждой группе после того, как произведена коррекция Бенджамини-Хохберга, причем так же, как и во второй части, считать только такие отличия, у которых abs(fold change) > 1.5.
In [57]:
pvalue_df['control_vs_neoplasia_rej_bh'] = smm.multipletests(pvalue_df.control_vs_neoplasia_pvalue, alpha=0.025, method='fdr_i')[0]
pvalue_df['neoplasia_vs_cancer_rej_bh'] = smm.multipletests(pvalue_df.neoplasia_vs_cancer_pvalue, alpha=0.025, method='fdr_i')[0]
In [60]:
pvalue_df.control_vs_neoplasia_rej_bh.value_counts()
Out[60]:
In [54]:
pvalue_df[(pvalue_df.control_vs_neoplasia_rej_bh) & (pvalue_df.control_vs_neoplasia_fold_change > 1.5)].shape
Out[54]:
In [55]:
pvalue_df[(pvalue_df.neoplasia_vs_cancer_rej_bh) & (pvalue_df.neoplasia_vs_cancer_fold_change > 1.5)].shape
Out[55]:
In [ ]: