edgeR中的那些坑（glmQLFit和glmQLFTest）

library(edgeR)library(statmod)rawdata<-read.delim("Back_skin.counts.txt", sep="\t", check.names=FALSE, stringsAsFactors=FALSE, row.names=1)y<-DGEList(counts=cbind(rawdata[,2],rawdata[,3],rawdata[,4],rawdata[,1],rawdata[,5],rawdata[,6]), genes=row.names(rawdata))### for counts filterkeep <- rowSums(cpm(y) >0.5) >=1table(keep)y <- y[keep, ,keep.lib.sizes=FALSE]#keep.lib.size=FALSE表示重新计算文库大小。y<-calcNormFactors(y)myfactors<-factor(c("WT","WT","WT","FST","FST","FST"))data.frame(Sample=colnames(y), myfactors)design<-model.matrix(~myfactors)rownames(design)<-colnames(y)y<-estimateDisp(y, design, robust=TRUE)fit<-glmFit(y, design)lrt<-glmLRT(fit)deg<-topTags(lrt, n=nrow(lrt$table), adjust.method="BH")$tablewrite.table(deg, file="Deglist_Back_skin_WT_vs_FST.xls", sep="\t", quote=FALSE, row.names=FALSE, col.names=TRUE)

网络上的差异分析代码数不甚数，贴个代码装个逼（还嫩呢）。不过说实话，搞科研这一行的，你去测个序，不会点差异分析怎么能行呢？但是问题来了，cummeRbund 和DESeq, edgeR, limma等等的区别是什么？我应该选哪一种进行差异分析呢？这些R包里面都运用了什么模型、什么检验方法进行差异分析的呢？如果搞不清这些问题，差异结果难以令人信服，讨论时，你也不会自行的说“我的分析是靠谱的”，很有可能被老板怼。

1. 基因差异表达分析 cummeRbund 和DESeq, edgeR, limma的区别是什么？（以前看过一篇别人比较的博文，“DEGseq无重复比较，DESeq有重复比较，edgR有无重复都可用，limma适于芯片，用前考虑batcheffect”，简单借用如下）：

        四款软件的输入文件不一样，这是很自然的。每个软件都需要准备各自的输入文件，格式不同那个肯定的，内容不同那就是取决于软件怎么分析了。结果不一样也很自然。分析方法不一样，结果不可能完全一致。四款软件都是R bioconductor的包。

1）cummeRbund，有参转录组tuxedo组件之一。

tuxedo=tophat+cufflinks+cummeRbund，cummeRbund特别用于这个分析流程。每个都是一套礼服的一部分。

2）DESeq，EMBL开发的。当然还有个DESeq2。

3）edgeR，WEHI开发的。

4）limma，芯片分析用的。虽然手册有用于RNA-seq分析的章节，不过我没用过limma做RNA-seq的差异分析。

其中1，2，3三款软件都是用于测序数据。一般就是比较DESeq和edgeR。Nature protocols有[Count-based differential expression analysis of RNA sequencing data using R and Bioconductor]（放一个能看的Count-based differential expression analysis of RNA sequencing…）由以上知，edgeR和DESeq用哪个都成。我是用edgeR，它的手册写得也足够好。

看了看别人的回答，发现FPKM是需要考虑的一点。

现在凡是用FPKM/RPKM求表达量的软件都别用了，只用拿TPM算的。具体原因见RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome；Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples现在涉及表达量的软件一般都会转算TPM的。Ref:trinityrnaseq/trinityrnaseq转录本定量一节。

2.edgeR分析，glmQLFit还是glmLRT：

大神们在bioconductor和biostars上有激烈的讨论，https://support.bioconductor.org/p/68629/，内容很多，我总结了一些要点：

1）如果在比较差异的过程中考虑到存在一些多变的、不确定的因素，最好采用glmQLFit和glmQLFTest（即perform quasi-likelihood F-tests），这一方法适用于大多数分析（RNAseq、chip-seq，Hi-C）。

2）glmFit和glmLRT（即likelihood ratio tests），适用于无重复的差异比较。（实际上改参数对应的方法只用了卡方检验：chi-squared test）。

3）基于实验设计完善（一般基于很好的重复），有人提出了大多数情况下，glmQLFit and glmQLFTest优于glmLRT进行差异分析。只有少数没有重复存在的实验设计中，采用glmLRT（以先验估计拟合二项分布结合广义线性模型（GLMs）寻找差异基因）。

4）有人认为，表达量数据在ERCC normalization（External RNA Control Consortium，ERCC spike-in）之后（即使是在有重复的差异比较下），LRT（glmLRT）才是做差异分析更好的选择。然而在阅读了这篇文章（https://www.nature.com/articles/nbt.2931）之后，又有其他声音说ERCC normalization (for bulk RNA-seq data, presumably?) can yield odd behaviour prior to DE analyses。【涉及ERCC方法的这篇文章我稍后再来解读】。

3.总之：

做差异优先考虑DESeq2和edgeR，且优先考虑使用glmQLFit，如果遇到miRNA（一般表达数目较少，相比于一般mRNA表达数目而言）分析或者没有重复的比较，则优先使用DEGseq和glmLRT（edgeR）。前面还提到了batcheffect，对于这个问题，在做差异时设计好你的design即可；limma等removebatcheffect的结果适用于一些作图和展示，并不是用于差异分析！！！