edgeR使用學習【轉載】

轉自：http://yangl.net/2016/09/27/edger_usage/

1.Quick start

edgeR使用學習【轉載】

2. 利用edgeR分析RNA-seq鑒别差異表達基因：

#加載軟體包
library("edgeR",verbose=0);

# 1. 載入資料 讀取read count數
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. 建構分組變量
#分為 Control組和DHT組 分别為4個和3個重複
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. 建立基因表達清單 進行标準化因子計算
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 ...

#這裡主要是通過圖形的方式來展示實驗組與對照組樣本是否能明顯的分開
#以及同一組内樣本是否能聚的比較近 即重複樣本是否具有一緻性
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. 定義差異表達基因與基本統計
summary(de <- decideTestsDGE(et)); # 預設選取FDR = 0.05為門檻值

#輸出
# [,1] 
# -1 2094 #顯著下調
# 0 12060 #沒有顯著差異
# 1 2340 #顯著上調

#圖形展示檢驗結果
detags <- rownames(y)[as.logical(de)];
plotSmear(et, de.tags=detags);
abline(h=c(-1, 1), col="blue");

//這個是分為 Control組和DHT組，檢驗這兩組的差異表達基因。

//中間又一步是去除表達量過低的基因。

讀取read count數
建構分組變量
建立基因表達清單進行标準化因子計算，過濾表達量偏低的基因，進行标準化因子計算預設使用TMM方法
估計離散度
通過檢驗來鑒别差異表達基因
定義差異表達基因與基本統計

轉載于:https://www.cnblogs.com/BlueBlueSea/p/10295712.html