edgeR入門(mén)

world663 2018-01-08

展開(kāi)全文

edgeR包主要是用于利用來(lái)自不同技術(shù)平臺(tái)的read數(shù)（包括RNA-seq，SAGE或者ChIP-seq等）來(lái)鑒別差異表達(dá)或者差異標(biāo)記（ChIP-seq）。主要是利用了多組實(shí)驗(yàn)的精確統(tǒng)計(jì)模型或者適用于多因素復(fù)雜實(shí)驗(yàn)的廣義線性模型。所以有時(shí)作者也把前者叫做“經(jīng)典edgeR”，后者叫做”廣義線性模型 edgeR“。這里定義的read數(shù)是可以指基因水平、外顯子水平、轉(zhuǎn)錄本水平或者標(biāo)簽水平等，這個(gè)由用戶根據(jù)自己數(shù)據(jù)分析的實(shí)際需要而定。這里作者也列舉了一些差異表達(dá)鑒定方面的文獻(xiàn)：包括edgeR剛發(fā)布時(shí)的文獻(xiàn)–“edgeR: a Bioconductor package for differential expression analysis of digital gene expression data”以及后來(lái)的一些改進(jìn)文章。

1 2	source("http:///biocLite.R") biocLite("edgeR")

快速開(kāi)始

對(duì)一個(gè)軟件不是很熟悉時(shí)，一般都可以看看它的快速簡(jiǎn)介，就能進(jìn)行運(yùn)算了。edgeR的基本計(jì)算步驟就是：
經(jīng)典edgeR操作步驟：

先讀取read數(shù)（比如說(shuō)RNA-seq數(shù)據(jù)，可以利用Bowtie、Tophat及htseq-counts這么個(gè)流程獲?。粋€(gè)行為基因（或者其它特征）列為樣本元素為read數(shù)的矩陣；
構(gòu)建分組變量，主要是對(duì)樣本進(jìn)行分組（分為處理組和對(duì)照組，利用factor函數(shù)）
建立基因表達(dá)列表，利用DEGList函數(shù)，參數(shù)counts為read數(shù)矩陣，group為上一步的分組變量，還有其它參數(shù)，后面有介紹
估計(jì)與計(jì)算各種其它參數(shù)，計(jì)算各樣本內(nèi)的標(biāo)準(zhǔn)化因子（calcNormFactors函數(shù)）（由于RNA-seq測(cè)序過(guò)程中技術(shù)因素可能會(huì)產(chǎn)生樣本特異性效應(yīng)，所以標(biāo)準(zhǔn)化一般是在這種情況下才使用的，也主要是測(cè)序深度即和文庫(kù)大小有關(guān)，可以通過(guò)DEGList函數(shù)的lib.size=colSums(counts)來(lái)設(shè)置。），估計(jì)生物學(xué)變異系數(shù)（BCV，Biological cofficient of variation，在經(jīng)典edgeR中，BCV就是負(fù)二項(xiàng)分布中的離散度的平方根，所以估計(jì)BCV也是估計(jì)離散度，edgeR先認(rèn)為所有的基因有同樣的離散度，estimateCommonDisp算得，在此基礎(chǔ)上計(jì)算基因間范圍內(nèi)的離散度，estimateTagwiseDisp算得）
通過(guò)檢驗(yàn)來(lái)計(jì)算鑒別差異表達(dá)基因

x <- read.delim("fileofcounts.txt",row.names="Symbol") # 讀取reads count文件

group <- factor(c(1,1,2,2)) #分組變量前兩個(gè)為一組，后一個(gè)為一組，每個(gè)有兩個(gè)重復(fù)

y <- DGEList(counts=x,group=group) # 構(gòu)建基因表達(dá)列表

y <- calcNormFactors(y) # 計(jì)算樣本內(nèi)標(biāo)準(zhǔn)化因子

y <- estimateCommonDisp(y) #計(jì)算普通的離散度

y <- estimateTagwiseDisp(y) #計(jì)算基因間范圍內(nèi)的離散度

et <- exactTest(y) # 進(jìn)行精確檢驗(yàn)

topTags(et) # 輸出排名靠前的差異表達(dá)基因信息

廣義線性模型 edgeR操作步驟：

也是要讀取raw read count 數(shù)據(jù)（同上）；
構(gòu)建實(shí)驗(yàn)設(shè)計(jì)變量也就是一個(gè)矩陣；
構(gòu)建基因表達(dá)列表并進(jìn)行標(biāo)準(zhǔn)化處理
利用實(shí)驗(yàn)設(shè)計(jì)變量及基因表達(dá)列表來(lái)估計(jì)估計(jì)生物變異的離散度（dispersion）、基因間離散度估算
擬合廣義線性模型并利用似然比檢驗(yàn)來(lái)鑒別差異表達(dá)基因

rawdata <- read.delim("TableS1.txt", check.names=FALSE, stringsAsFactors=FALSE); #讀取原始文件

y <- DGEList(counts=rawdata[,4:9], genes=rawdata[,1:3]); # 建立基因表達(dá)列表

y$samples$lib.size <- colSums(y$counts); #設(shè)置各個(gè)樣本的庫(kù)大小

y <- calcNormFactors(y); #利用庫(kù)的大小來(lái)進(jìn)行標(biāo)準(zhǔn)化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); #估計(jì)離散度

y <- estimateGLMTrendedDisp(y, design);

y <- estimateGLMTagwiseDisp(y, design);

fit <- glmFit(y, design); #擬合廣義線性模型

lrt <- glmLRT(fit); topTags(lrt); #利用似然比檢驗(yàn)來(lái)檢驗(yàn)差異表達(dá)基因

topTags(lrt); #輸出靠前的差異表達(dá)基因

實(shí)例：

下面通過(guò)文檔中提供的一個(gè)RNA-seq的例子來(lái)具體說(shuō)下edgeR的使用流程（經(jīng)典edgeR），GLM方法的就不講了，更復(fù)雜的見(jiàn)edgeR的文檔。
數(shù)據(jù)： pnas_expression 文本數(shù)據(jù)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

#加載軟件包

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個(gè)和3個(gè)重復(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)建基因表達(dá)列表進(jìn)行標(biāo)準(zhǔn)化因子計(jì)算

y <- DGEList(counts=data[,1:7], group=targets$Treatment, genes=data.frame(Length=data[,8]));

colnames(y) <- targets$Label;

dim(y);

# [1] 37435 7

#過(guò)濾表達(dá)量偏低的基因

#基因在至少3個(gè)樣本中得count per million（cpm）要大于1

keep <- rowSums(cpm(y)>1) >= 3;

y <- y[keep,];

dim(y)

# [1] 16494 7

#重新計(jì)算庫(kù)大小

y$samples$lib.size <- colSums(y$counts);

#3. 進(jìn)行標(biāo)準(zhǔn)化因子計(jì)算默認(rèn)使用TMM方法

y <- calcNormFactors(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 ...

#這里主要是通過(guò)圖形的方式來(lái)展示實(shí)驗(yàn)組與對(duì)照組樣本是否能明顯的分開(kāi)

#以及同一組內(nèi)樣本是否能聚的比較近即重復(fù)樣本是否具有一致性

plotMDS(y);

#4. 估計(jì)離散度

y <- estimateCommonDisp(y, verbose=TRUE)

# Disp = 0.02002 , BCV = 0.1415

y <- estimateTagwiseDisp(y);

plotBCV(y);

#5. 通過(guò)檢驗(yàn)來(lái)鑒別差異表達(dá)基因

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. 定義差異表達(dá)基因與基本統(tǒng)計(jì)

summary(de <- decideTestsDGE(et)); # 默認(rèn)選取FDR = 0.05為閾值

#輸出

# [,1]

# -1 2094 #顯著下調(diào)

# 0 12060 #沒(méi)有顯著差異

# 1 2340 #顯著上調(diào)

#圖形展示檢驗(yàn)結(jié)果

detags <- rownames(y)[as.logical(de)];

plotSmear(et, de.tags=detags);

abline(h=c(-1, 1), col="blue");

關(guān)于樣本沒(méi)有重復(fù)時(shí)該怎么做

當(dāng)然最好還是有重復(fù)，實(shí)在沒(méi)有重復(fù)的話，edgeR的文檔也有提到，那就是選擇一個(gè)認(rèn)為合理的離散度的值。
通常如果是實(shí)驗(yàn)控制的好的人類(lèi)數(shù)據(jù)，那么選擇BCV=0.4，比較好的模式生物選擇BCV=0.1，技術(shù)重復(fù)的話選擇BCV=0.1。
這里作者舉了個(gè)簡(jiǎn)單的模擬數(shù)據(jù)的例子：

bcv <- 0.2

counts <- matrix( rnbinom(40,size=1/bcv^2,mu=10), 20,2)

y <- DGEList(counts=counts, group=1:2)

et <- exactTest(y, dispersion=bcv^2)

本站是提供個(gè)人知識(shí)管理的網(wǎng)絡(luò)存儲(chǔ)空間，所有內(nèi)容均由用戶發(fā)布，不代表本站觀點(diǎn)。請(qǐng)注意甄別內(nèi)容中的聯(lián)系方式、誘導(dǎo)購(gòu)買(mǎi)等信息，謹(jǐn)防詐騙。如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容，請(qǐng)點(diǎn)擊一鍵舉報(bào)。

轉(zhuǎn)藏 分享

QQ空間 QQ好友新浪微博微信

獻(xiàn)花（0） +1

來(lái)自： world663 > 《生信分析》

舉報(bào)/認(rèn)領(lǐng)