芯片差异分析就是 对每个基因做t-test
为什么要进行FDR校正
多重检验问题:同时识别两个总体的成千上万个基因中哪些基因的表达水平有显著性差异，错误的概率在增加
设定假设检验水准为α（通常取0.05），若p< α，则拒绝原假设（基因差异表达）；反之，接受原假设

检验一个基因，真阳性概率为0.95，基因为假阳性概率1-0.95；检验两个基因，真阳性概率为0.95^{2，至少一个基因为假阳性概率为1-0.95}2；……；检验n次，真阳性概率为0.95^{n，至少一个基因为假阳性概率为（1-0.95）}n,因此随着检验次数增多，假阳性率会越来越大

image.png

一.数据描述

image.png

二. 实验目的
通过控制FDR来减小假阳性率，提高我们筛选差异基因的准确率
三. matlab代码

%抽取奇数列A(:,1:2:end) 偶数列A(:,2:2:end)

filename='gene_exprs.txt'; %这个文件是通过数据数据预处理的(用R)，是一个20486*12的表达谱矩阵
gene_exprs=importdata(filename);
A=gene_exprs.textdata;
%表达谱三要素
gene_exprs.genesymbol=A(2:end,1); %提取genesymbol
gene_exprs.sampleData=A(1,2:end);%提取样本名
gene_exprs.label=[1;1;1;1;1;1;0;0;0;0;0;0]; %标签，1是被BV6处理的实验组，0是被DMSO处理的对照组
%以下是数据预处理
Imputedata=knnimpute(gene_exprs.data,15);%处理缺失值
maboxplot(Imputedata,gene_exprs.sampleData)
xlabel('Samples')
geneexprs.data=quantilenorm(Imputedata,'MEDIAN',1);
geneexprs.data=quantilenorm(Imputedata,'display',1);

%差异基因筛选BH，Bonferroni法
for i=1:size(gene_exprs.data,1)
[h,p,ci]=ttest2(gene_exprs.data(i,1:6),gene_exprs.data(i,7:12),0.05);
Sigtotal(i)=p;
end
Siggene=gene_exprs.genesymbol(Sigtotal<=0.05);
pval=Sigtotal;
function [sorted_p,sort_index,q]=FDR_base(pval,method)
  if strcmp(method,'BH')
      [sorted_p,sort_index]=sort(pvals);
     for i=1:size(gene_exprs.data,1);
q(i)=pval(i)*20486/sort_index(i);
     end
  end
  Siggene_BH=gene_exprs.genesymbol(q<=0.05);
  if strcmp(method,'Bonferroni')
       for i=1:size(gene_exprs.data,1);
           q(i)=pval(i);
       end
  end
  Siggene_bf=gene_exprs.genesymbol(q<=(0.05/20486));
Save
save  D:\实验\实验五\DEGs_BH Siggene_BH %先写保存路径，再写需要保存的矩阵
save  D:\实验\实验五\DEGs_bf Siggene_bf