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