FDR校正的程序实现及严格程度对比

 

FDR校正的程序实现及严格程度对比

 

前言

 

        做统计分析就离不开P value<0.05,而写过科研文章的人也都知道没有经过FDR校正的P值就像一盘散沙,不用风吹,走两步自个儿就散了。 那么FDR校正这个让人又爱又恨的东西是什么呢?又是如何实现呢?

        原理是这样:设总共有m个候选基因,每个基因对应的p值从小到大排列分别是 p(1),p(2),...,p(m),则若想控制fdr不能超过q(如0.05),则只需找到最大的正整数i,使得 p(i)<= (i*q)/m.然后,挑选对应p(1),p(2),...,p(i)的基因做为差异表达基因,这样就能从统计学上保证fdr不超过q。计算方法参考:http://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html

        没听明白?再说简单点就是对统计分析结果产生的P值做一次筛选,对原始P值<0.05这个显著性的真实性用FDR校正后的新P值做判定。

         下面对FDR校正的程序实现基于MATLAB和R分别进行介绍。

 

(一) MATLAB实现

 

mafdr函数

 

1. FDR = mafdr(PValues);                                                        

%最简单的实现方式,基于Storey procedure ( introduced by Storey, 2002),适用于P值数量>1000的情况,否则原则上会崩溃。我用MATLAB测试过,会出现warning,但不会报错。严格程度较低,如果你的 ttest P值不是特别显著(0.01-0.05),可以用这个试试,或许可以过FDR校正.

 

2.FDR=mafdr(P,'BHFDR', true);                                              

%基于linear step-up (LSU) procedure (introduced by Benjamini and Hochberg, 1995)。最常见的FDR校正方式,严格程度较高,但比
Bonferroni校正低,适用于 ttest P值显著(<0.01)。

注:

BHFDR的计算过程

1. 将一系列的p值按照从大到小排序,然后利用下述公式计算每个p值所对应的FDR值:

公式:FDR = p * (n/i), p是pvalue,n是p值个数,最大的P值的i值为n,第二大则是n-1,依次至最小为1。

2. 将计算出来的FDR值作为新p值,如果某一个p值所对应的FDR值大于前一位p值(更大的p值)所对应的FDR值,则放弃公式计算出来的FDR值,选用与它前一位相同的值,因此会产生连续相同FDR值的现象;反之则保留计算的FDR值。

3. 返回p值对应的FDR值。

 

3.FDR=mafdr(P,'Lambda', [0.01:0.01:0.95]);                          

 %指定调整参数λ,用于估计零假设为真的先验概率,Lambda Value可以是:(1)0-1内任意一值  (2)4个以上序列,或以示例的矩阵形式:[first:incr:last],mafdr函数自动选择最优参数。严格程度比1略低,我用第一个不能过校正但是用这个居然可以过,为FDR苦恼得同志们可以被解救了。

 

4.FDR= mafdr(P, 'Method', MethodValue, [0.01:0.01:0.95]);   

%对选择最优参数lambda的方法进行选择,'bootstrap' (default),'polynomial'。第3种方法的延伸。

 

(二) R语言实现

 

p<-c(p1,p2,...,pn)                                                                    %括号内载入P值序列,","间隔

p.adjust(p,method="method value",length(p))                    %通过修改method value变换FDR方法,包括"holm","horchberg","hommel","bonferroni","BH","BY","fdr"。从严格程度上"bonferroni"最可怕,"BH","fdr"与MATLAB的1,2相同,其余在之间。

 

FDR校正的程序实现先说到这里,如果你有其他的方法请在下方留言评论或发邮件与我交流,我的邮箱是liuyuchen0020@163.com

 

 

### 在 MATLAB 中实现 FDR 校正 #### 使用内置函数 `mafdr` MATLAB 提供了一个名为 `mafdr` 的内置函数来执行 FDR 控制。此函数主要用于微阵列数据分析,但也适用于其他类型的多假设检验场景[^3]。 ```matlab % 假设 pvals 是一个包含多个p值的向量 [pFDR, qValues] = mafdr(pvals); ``` 上述代码返回两个变量: - `pFDR`: 经过 FDR 校正后的 p 值。 - `qValues`: 对应于每个原始 p 值的 Q 值 (即最小的 FDR 阈值,在该阈值下原假设可以被拒绝)。 #### 手动实现 Benjamini-Hochberg 方法 如果希望更灵活地控制过程或理解内部机制,则可以通过手动编写 BH 法来进行 FDR 校正: ```matlab function [fdr_corrected_pvalues] = bh_fdr_correction(p_values, alpha) % Sort the p-values and keep track of original indices. [~, idx] = sort(p_values); sorted_p_values = p_values(idx); n = length(sorted_p_values); % Calculate critical values based on rank position. ranks = 1:n; crit_vals = (ranks / n) * alpha; % Find largest i such that P(i)<=i/m*alpha where m is total number of tests. max_i = find(sorted_p_values <= crit_vals, 1, 'last'); if isempty(max_i) fdr_corrected_pvalues = ones(size(p_values)); % No significant results found. else threshold = sorted_p_values(max_i); % Apply correction to all p-values. corrected_sorted_p_values = min([sorted_p_values ./ ((1:n)' / n), ones(n, 1)], [], 2); % Restore order according to input vector. [~, inv_idx] = sort(idx); fdr_corrected_pvalues = corrected_sorted_p_values(inv_idx); end end ``` 这段自定义代码实现了经典的 Benjamini 和 Hochberg 提出的线性逐步方法用于控制错误发现率[FDR][^2]。
评论 9
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值