edgeR

edgeR:Empirical Analysis of Digital Gene Expression Data in R

一个差异性分析的R包，用于RNA-seq或DNA甲基化等相关技术分析。

其原理利用广义线性模型对每个基因或者甲基化位点建模，然后直接比较线性模型的参数。

输入要求：必须是支持该位点的原始read count，而不是经过normalization计算的结果。

对于RNA-seq可以是htseq-count的结果。

对于甲基化分析可以是bismark的结果。

 edgeR是以DEGList的格式储存数据，DEGList是一个以list为基础的数据格式，list的所有方法其都可以使用。

1.构建DEGList：用DEGList（）构建。

>y<-DGEList(counts=x)  #x是read counts的matrix或data.flame。
>group<-c(1,1,2,2)    #关于sample属于哪一个group。
>y<-DGEList(counts=x,group=group) 

 DEGList中主要包括一个 $counts，一个 $samples  ,还有一个可选的$genes （注释）。

$sample:

lib.size:默认为count的每一列总和，代表了该样品的总深度。

2。过滤：

生物学上看，一个基因要被表达成蛋白或是其他的生物功能，则它的表达量应该达到一个最低的水平。
所以在进一步分析前，应该过滤掉一些low counts的基因。这里用cpm（count per million)来表示基因的counts水平。

cpm的计算举例：比如在$count中一个点是10 ，该位点对应的sample的lib.size=50000.

该点经过cpm计算后得到的值（X）：10/50000=X/1M

>keep<-rowSums(cpm(y)>1)>=2  #>=2表示每个group中的samples数最少是2.
>y<-y[keep, ,keep.lib.sizes=FALSE]

3.TMM标准化：

在treated和untreated样品中，常常会有少量的基因在treated样品中高表达，但在untreated样品则正常。在treated样品中，高表达的基因的reads会占据一大部分的library size，

而导致剩余基因被错误的判断为下调。

>y<-calcNormFactors(y)
>y$samples

4.设计矩阵（design matrix）

补充：表达式(~)

~0+group :不包括比较矩阵

~group:包括了比较矩阵

design<-model.matrix(~group)

5.离散度的检测(dispersion)：

补充：

什么是离散度：

 

y<-estimateGLMCommonDisp(y,design,verbose=TRUE)

y<-estimateGLMTrendedDisp(y, design)

y<-estimateGLMTagwiseDisp(y, design)

 

 

6.差异性基因分析

to perform quasi-likelihood F-tests:

fit <- glmQLFit(y,design)

qlf <- glmQLFTest(fit,coef=2)

topTags(qlf)

 

to perform likelihood ratio tests:

fit<-glmFit(y, design)

lrt<-glmLRT(fit)

topTags(lrt)#前10个差异表达基因，内容分别为：genes:logFC:logCPM:F:PValue:FDR

#FC:fold change,，一般取两个样本的比值的均值，由于两个样本之间差异过大，为了缩小数值之间的差异，做对数转换。Log2FC一般选择大于1的。

P VALUE:统计学检验变量，代表差异显著性，一般认为p值小于0.05代表具有显著性差异。

FDR:false discovery rate,

 

 

7.筛选出符合的基因

dt<-decideTestsDGE(lrt)

#decideTestsDGE(object,p.value=0.05,lfc=0)

isDE<-as.logical(dt)

DEnames<-rownames(y)[isDE]

 参考：

http://blog.sina.com.cn/s/blog_5d188bc40102vwci.html

http://blog.sciencenet.cn/blog-508298-776802.html