不同差异分析方法拿到的上下调基因影响什么了？

文献里面提到了一个表达量矩阵数据集的两个可能的误差来源，首先是该队列样品收集就耗费了四年所以有批次效应，其次队列里面的性别失衡，这两个问题都是需要解决的。所以研究者们采用了ANOVA model 很严谨的去判别差异基因，方法学如下所示：

这是一个表达量芯片数据集：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE117261，是很经典的的两分组：58 PAH and 25 control lung tissues，然后我也默认走了标准差异分析，以及读取了作者的文献附件里面的差异分析结果，简单的对比了一下：

load('GSE117261/anno_DEG.Rdata')
table(DEG$g)
head(DEG)
DEG=DEG[!duplicated(DEG$name),]
rownames(DEG)=DEG$name
cg = DEG$name[DEG$g!='stable']
my_deg = DEG
load('GSE117261/step1-output.Rdata') 

paper_deg =data.table::fread('paper.txt',data.table = F)
paper_deg=na.omit(paper_deg)
head(paper_deg)
str(head(paper_deg))
paper_deg=paper_deg[nchar(paper_deg[,1])>1,]
paper_deg=paper_deg[!duplicated(paper_deg[,1]),]
rownames(paper_deg)=paper_deg[,1]
colnames(paper_deg)=c('Symbol','FC','FDR')

ids= intersect(rownames(my_deg),rownames(paper_deg))
plot(my_deg[ids,'logFC'],paper_deg[ids,'FC'])
plot(my_deg[ids,'logFC'],log2(paper_deg[ids,'FC']))

我和作者的logFC是基本上没有差异的，不过作者在文章附件给出来的是没有log的FC，然后我看了看我们不同方法判别差异分析的统计学显著的上下调基因的一致性，如下所示：

在作者的标准里面只需要 false discovery rate [FDR] q value , 0.001  就是统计学显著的差异基因，代码如下所示：

pvalue_t=0.001
paper_deg$g=ifelse(paper_deg$FDR > pvalue_t,'stable', #if 判断：如果这一基因的P.Value>0.01，则为stable基因
            ifelse( paper_deg$FC > 1,'up', #接上句else 否则：接下来开始判断那些P.Value<0.01的基因，再if 判断：如果logFC >1.5,则为up（上调）基因
                    ifelse( paper_deg$FC 1,'down','stable') )#接上句else 否则：接下来开始判断那些logFC <1.5 的基因，再if 判断：如果logFC <1.5，则为down（下调）基因，否则为stable基因
)
table(paper_deg$g)

而我们的表达量芯片默认的差异分析需要同时卡logFC，所以有火山图如下所示：

从火山图可以看到我给出的阈值是很奇怪，因为我采用了 cut_logFC logFC)) )  这个统计学 概念，我们可以细致的检查一下我们的差异分析冲突的基因情况。

df=data.frame(
  gene=ids,group=paste0(my_deg[ids,'g'], paper_deg[ids,'g'] )
)
dfl = split(df$gene,df$group)
symbols_list = dfl[-4]
df= melt(dat[unlist(lapply(symbols_list, head)),])
df$group  = rep(group_list,each=36)
head(df)
p "group", y = "value",
               color = "group", palette = "jco",
               add = "jitter")+facet_wrap(~Var1)
#  Add p-value
p + stat_compare_means() 
  theme_bw() 

比如我们先看看3种下调基因：

最上面的基因是我们两个方法都认为是下调基因，但是中间的是我们的方法判定为下调，但是作者的ANOVA model 反对的，最下面的基因就是反过来的。

然后看看上调基因的冲突情况：

基本上就是一样的情况啦，最上面的基因是我们两个方法都认为是上调基因，但是中间的是我们的方法判定为上调，但是作者的ANOVA model 反对的，最下面的基因就是反过来的。

这个时候，虽然我们搞清楚了，确实是两个方法有冲突，而且冲突的原因就在于对统计学p值的容忍度以及对变化倍数的阈值选择。我们的方法倾向于选择那些变化倍数比较大但是可能统计学不那么显著的，而作者的ANOVA model 恰好就反过来了。

这样的话，我们就产生了6种不同的基因列表，是可以进行生物学功能注释的，代码如下所示：