2019-09-02R包第一天：差异分析包（一）edgeR
2019-10-01 本文已影响0人
一只烟酒僧
本文转载自：陈洪瑜的简书
英文教程请参考：bioconductor
edgeR主要是利用了适用于多组实验的精确统计模型或者适用于多因素复杂实验的广义线性模型。前者叫做“经典edgeR”， 后者叫做”广义线性模型 edgeR“。
其中定义的read数是可以指基因水平、外显子水平、转录本水平或者标签水平等。
经典edgeR的使用方法
x<-read.delim("fileofcounts.txt",row.names="Symbol")# 读取reads count文件
group<-factor(c(1,1,2,2))#分组变量 前两个为一组， 后一个为一组， 每个有两个重复
y<-DGEList(counts=x,group=group)# 构建基因表达列表
y<-calcNormFactors(y)# 计算样本内标准化因子
y<-estimateCommonDisp(y)#计算普通的离散度
y<-estimateTagwiseDisp(y)#计算基因间范围内的离散度
et<-exactTest(y)# 进行精确检验
topTags(et)# 输出排名靠前的差异表达基因信息
广义线性模型 edgeR （可用于无重复样本）
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); #设置各个样本的库大小
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. 载入数据 读取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");
当没有重复时如何计算：
方法1：根据经验给定一个值(BCV, square-root-dispersion). edgeR给的建议是, 如果你是人类数据, 且实验做的很好(无过多的其他因素影响), 设置为0.4, 如果是遗传上相似的模式物种, 设置为0.1, 如果是技术重复, 那么设置为0.01
y_bcv <- y
bcv <- 0.4
et <- exactTest(y_bcv,dispersion = bcv ^2)
方法2: 根据已知一些不会发生改变的基因推测dispersion. 假设你已经知道了一些基因是不会发生变化，例如管家基因，那么我们就可以通过它们来预测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)
