edgeR
主要是利用了多組實驗的精確統(tǒng)計模型或者適用于多因素復(fù)雜實驗的廣義線性模型戏售。
前者叫做“經(jīng)典edgeR”侨核， 后者叫做”廣義線性模型 edgeR“。
定義的read數(shù)是可以指基因水平灌灾、外顯子水平搓译、轉(zhuǎn)錄本水平或者標簽水平等
安裝：
library(BiocManager）
BiocManager::install("edgeR")
library(edgeR)
經(jīng)典edgeR
x<-read.delim("fileofcounts.txt",row.names="Symbol")# 讀取reads count文件
group<-factor(c(1,1,2,2))#分組變量 前兩個為一組， 后一個為一組锋喜， 每個有兩個重復(fù)
y<-DGEList(counts=x,group=group)# 構(gòu)建基因表達列表
y<-calcNormFactors(y)# 計算樣本內(nèi)標準化因子
y<-estimateCommonDisp(y)#計算普通的離散度
y<-estimateTagwiseDisp(y)#計算基因間范圍內(nèi)的離散度
et<-exactTest(y)# 進行精確檢驗
topTags(et)# 輸出排名靠前的差異表達基因信息
廣義線性模型 edgeR? 可用于無重復(fù)樣本
rawdata<-read.delim("TableS1.txt",check.names=FALSE,stringsAsFactors=FALSE);??#讀取原始文件
y<-DGEList(counts=rawdata[,4:9],genes=rawdata[,1:3]);??# 建立基因表達列表
y$samples$lib.size<-colSums(y$counts);??#設(shè)置各個樣本的庫大小
y<-calcNormFactors(y);??#利用庫的大小來進行標準化TMM
Patient<-factor(c(8,8,33,33,51,51));??#因素1
Tissue<-factor(c("N","T","N","T","N","T"));#因素2
data.frame(Sample=colnames(y),Patient,Tissue);
design<-model.matrix(~Patient+Tissue);??#建立分組變量
rownames(design)<-colnames(y);
y<-estimateGLMCommonDisp(y,design,verbose=TRUE);??#估計離散度
y<-estimateGLMTrendedDisp(y,design);??
y<-estimateGLMTagwiseDisp(y,design);
fit<-glmFit(y,design);??#擬合廣義線性模型
lrt<-glmLRT(fit);topTags(lrt);??#利用似然比檢驗來檢驗差異表達基因
topTags(lrt);????#輸出靠前的差異表達基因
edgeR示例
#加載軟件包
library("edgeR",verbose=0);
# 1. 載入數(shù)據(jù) 讀取read count數(shù)
data<-??read.delim("pnas_expression.txt",row.names=1,stringsAsFactors=FALSE);
head(data);
#輸出
# lane1 lane2 lane3 lane4 lane5 lane6 lane8??len
# ENSG00000215696???? 0???? 0???? 0???? 0???? 0???? 0???? 0??330
# ENSG00000215700???? 0???? 0???? 0???? 0???? 0???? 0???? 0 2370
# ENSG00000215699???? 0???? 0???? 0???? 0???? 0???? 0???? 0 1842
# ENSG00000215784???? 0???? 0???? 0???? 0???? 0???? 0???? 0 2393
# ENSG00000212914???? 0???? 0???? 0???? 0???? 0???? 0???? 0??384
# ENSG00000212042???? 0???? 0???? 0???? 0???? 0???? 0???? 0?? 92
dim(data);
# [1] 37435???? 8
#2. 構(gòu)建分組變量
#分為 Control組和DHT組??分別為4個和3個重復(fù)
targets<-data.frame(Lane=c(1:6,8),Treatment=c(rep("Control",4),rep("DHT",3)),
??????????????????????Label=c(paste("Con",1:4,sep=""),paste("DHT",1:3,sep="")));
targets
#輸出
#?? Lane Treatement Label
# 1????1????Control??Con1
# 2????2????Control??Con2
# 3????3????Control??Con3
# 4????4????Control??Con4
# 5????5????????DHT??DHT1
# 6????6????????DHT??DHT2
# 7????8????????DHT??DHT3
#3. 創(chuàng)建基因表達列表 進行標準化因子計算
y<-DGEList(counts=data[,1:7],group=targets$Treatment,genes=data.frame(Length=data[,8]));
colnames(y)<-targets$Label;
dim(y);
# [1] 37435???? 7
#過濾表達量偏低的基因
#基因在至少3個樣本中得count per million（cpm）要大于1
keep<-rowSums(cpm(y)>1)>=3;
y<-y[keep,];
dim(y)
# [1] 16494 7
#重新計算庫大小
y$samples$lib.size<-colSums(y$counts);
#3. 進行標準化因子計算 默認使用TMM方法
y<-calcNormFactors(y);
y
#輸出
# An object of class "DGEList"
# $counts
# Con1 Con2 Con3 Con4 DHT1 DHT2 DHT3
# ENSG00000124208??478??619??628??744??483??716??240
# ENSG00000182463?? 27?? 20?? 27?? 26?? 48?? 55?? 24
# ENSG00000124201??180??218??293??275??373??301?? 88
# ENSG00000124207?? 76?? 80?? 85?? 97?? 80?? 81?? 37
# ENSG00000125835??132??200??200??228??280??204?? 52
# 16489 more rows ...
#
# $samples
# group lib.size norm.factors
# Con1???? 1?? 976847????1.0296636
# Con2???? 1??1154746????1.0372521
# Con3???? 1??1439393????1.0362662
# Con4???? 1??1482652????1.0378383
# DHT1???? 1??1820628????0.9537095
# DHT2???? 1??1831553????0.9525624
# DHT3???? 1?? 680798????0.9583181
#
# $genes
# [1] 2131 5453 4060 3264??945
# 16489 more rows ...
#這里主要是通過圖形的方式來展示實驗組與對照組樣本是否能明顯的分開
#以及同一組內(nèi)樣本是否能聚的比較近 即重復(fù)樣本是否具有一致性
plotMDS(y);
#4. 估計離散度
y<-estimateCommonDisp(y,verbose=TRUE)
# Disp = 0.02002 , BCV = 0.1415
y<-estimateTagwiseDisp(y);
plotBCV(y);
#5. 通過檢驗來鑒別差異表達基因
et<-exactTest(y);
top<-topTags(et);
top
#輸出
# Comparison of groups:??DHT-Control
# Length????logFC????logCPM????????PValue?????????? FDR
# ENSG00000151503?? 5605 5.816156??9.716866??0.000000e+00??0.000000e+00
# ENSG00000096060?? 4093 5.004454??9.950606??0.000000e+00??0.000000e+00
# ENSG00000166451?? 1556 4.683425??8.850401 2.297717e-249 1.263285e-245
# ENSG00000127954?? 3919 8.120955??7.216393 1.534440e-229 6.327264e-226
# ENSG00000162772?? 1377 3.316701??9.743797 7.975243e-216 2.630873e-212
# ENSG00000115648?? 2920 2.598440 11.474677 6.984860e-180 1.920138e-176
# ENSG00000116133?? 4286 3.244446??8.791930 1.290432e-174 3.040627e-171
# ENSG00000113594??10078 4.111120??8.055613 3.115276e-161 6.422921e-158
# ENSG00000130066????868 2.609899??9.989778 6.009018e-155 1.101253e-151
# ENSG00000116285?? 3076 4.201846??7.361640 6.299060e-149 1.038967e-145
#6. 定義差異表達基因與基本統(tǒng)計
summary(de<-decideTestsDGE(et));# 默認選取FDR = 0.05為閾值
#輸出
# [,1]
# -1??2094?? #顯著下調(diào)
# 0??12060?? #沒有顯著差異
# 1?? 2340?? #顯著上調(diào)
#圖形展示檢驗結(jié)果
detags<-rownames(y)[as.logical(de)];
plotSmear(et,de.tags=detags);
abline(h=c(-1,1),col="blue");
當沒有重復(fù)時如何計算：
方法1：根據(jù)經(jīng)驗給定一個值(BCV, square-root-dispersion). edgeR給的建議是, 如果你是人類數(shù)據(jù), 且實驗做的很好(無過多的其他因素影響), 設(shè)置為0.4, 如果是遺傳上相似的模式物種, 設(shè)置為0.1, 如果是技術(shù)重復(fù), 那么設(shè)置為0.01
y_bcv <- y
bcv <- 0.4
et <- exactTest(y_bcv,dispersion = bcv ^2)
方法2: 根據(jù)已知一些不會發(fā)生改變的基因推測dispersion. 假設(shè)你已經(jīng)知道了一些基因是不會發(fā)生變化，例如管家基因，那么我們就可以通過它們來預(yù)測dispersion.
nosig <- names(res@.Data[res@.Data ==0,])
gene_name <- sample(nosig,500,replace = FALSE)
y1_gene <- rownames(y1)
housekeeping <- which(y1_gene %in% gene_name)
y1 <- y
y1$samples$group <- 1
y0 <- estimateDisp(y1[housekeeping,],trend = "none",tagwise =FALSE)
y$common.dispersion <- y0$common.dispersion
design <- model.matrix(~group)
fit <- glmFit(y,design)
lrt <- glmLRT(fit)
genes <- decideTestDGE(lrt,p.value=0.05,lfc = 0)
補充：表達式(~)
~0+group :不包括比較矩陣
~group:包括了比較矩陣
design<-model.matrix(~group)
