今天是生信星球陪你的第991天
公众号里的文章大多数需要编程基础，如果因为代码看不懂，而跟不上正文的节奏，可以来找我学习，相当于给自己一个新手保护期。我的课程都是循环开课，点进去咨询微信↓
生信分析直播课程(9月30日下一期)
生信新手保护学习小组（10月初下一期）
单细胞陪伴学习小组（9.16下一期）
作者：一条没有梦想的咸鱼
24.9.13 投稿
经常在期刊投稿过程中有看到杂志要求提供多重假设检验校正的结果。原始P值、校正后P值，一直没怎么特别搞清楚。最近看到小洁老师在更新Python相关内容，也自学了一些Python统计的东西，与大家分享一下。
1.啥时候需要做P值校正呢？
当同一个数据集，进行的假设检验的次>=2次的时候就需要进行P值校正了。典型的场景就是差异分析，同时对成千上万个基因做了假设检验。还有多组数据两两比较的时候，通常都是做了多次假设检验。
2.为啥要做P值得校正呢？
我们知道通常假设检验P值的阈值一般设为P<0.05。统计学上一般把发生率小于5%这一事件称为小概率事件，认为在一次试验中几乎不会发生。P<0.05的意义就是假阳性（一类错误）的发生的概率小于5%，在一次试验中几乎不会发生。但当我们同时进行多次假设检验的时候假阳性发生的概率就会随着检验次数的增加而逐渐累积，超过5%这时候就需要对P值进行校正了，防止出现过多假阳性的结果。
3.怎样做P值校正呢?
下面介绍两种常用校正方法与Python实现
#载入pandas包
import pandas as pd
#读入deseq2的差异分析结果
data=pd.read_csv(r"C:\Users\result_M.csv",index_col=0)#index_col=0参数类似于R中的row.name=1,将第一列设为数据框行名（index）
data
3.1手动计算-bonferroni校正
Bonferroni校正这是一种严厉而保守的校正方法，对假阳性的控制非常严格
data['bonferroni']=data['pvalue']*len(data) #将原始P值乘以假设检验的次数（表格中基因的数量）就为校正后的P值
data
花花补充：可以计算结果中看到有一些大于1的值，在后面的statmodels包里是把这些值改为了1，因为p值范围是0-1
3.2 手动计算-BH校正
BH校正即Benjamini and Hochberg法，与Bonferroni相比是一种更为温和的校正方法，在数据分析更为常用（deseq2默认使用BH校正）
要先排序再计算，因为公式里用到了排序后的位置
data.sort_values(by='pvalue',ascending=True,inplace=True)#按原始P值大小，对基因进行升序排列
data.reset_index(inplace=True)#升序排列后将rowname恢复到数据列中，此时的rowname就是基因的具体排序
data
#计算BH校正后P值
data['BH_fdr']=data['pvalue']*len(data)/(data.index+1)
#即将bonferroni法校正的每个基因的p值除以它的排序就是BH校正后P值（这里给rowname统一加了1是因为python的索引是默认从0开始的不是1）
data#这里我们看到BH_fdr列我们计算的p值与deseq2计算的校正后P值（padj列）完全一样
3.3 statsmodels 包计算-bonferroni校正
其实使用statsmodels包可以快速计算，手动计算只是为了理解
from statsmodels.stats.multitest import multipletests
#bonferroni校正
bonferroni_corrected=multipletests(pvals=data['pvalue'],method='bonferroni')[1]#这里我们只需要在pvals= 参数输入原始p值，method参数选择相应校正方法就可以完成计算
data['bonferroni']=bonferroni_corrected
data
3.4 statsmodels 包计算-BH校正
#BH校正
BH__corrected=multipletests(pvals=data['pvalue'],method='fdr_bh')[1]#输入原始p值，method选择'fdr_bh'就好了
data['BH_fdr']=BH__corrected
data
PS:常见的软件输出结果都是直接校正P值，也有统计软件是P值不变，直接调整P值的阈值的，其实是一个道理哦，比如原始P值0.01校正为0.05，也可以原始P值不变还是0.01，但将P值的阈值变为0.01而不是通常的0.05。
