基因芯片数据分析（七）：edgeR差异分析实战案例
我们接下来通过一个案例介绍利用edgeR和DESeq2包进行差异分析，本文先介绍edgeR。
包的安装和加载
BiocManager::install( "edgeR" )
library( "edgeR" )读入数据
这里我们用的数据是一个Excel文件，是原始的counts文件（想运行案例，文末获取文件）
# 读取数据
counts <- read.table("gene_counts.xls", sep = "\t", header = T, row.names = 1)导入的数据我们通过head()函数查看前6行。
行名A1,A2,A3,B1,B2,B3为样本名，列名是基因名。
创建分组
设置实验组别，在基因芯片数据分析（五）：edgeR包的基本原理这篇文章中我们介绍基本原理时，有一步需要选择参考样本，在实际分析中，我们可以自己选择参考样本，一般都是对照组作为参考样本，在edgeR中，1代表control样本，2代表case样本。原始数据中有6个样本，对照和实验组各三个。所以我们可以创建一个分组向量。
# 创建分组
group <- c(rep(1,3),rep(2,3))
数据处理
# 创建DGEList类型变量
y <- DGEList(counts=counts, group=group)这一步相当于创建一列表。注意group中的顺序和counts中行名要对应，也就是对照组和实验组要指定正确。这里A1,A2,A3为control，B1,B2,B3为case。
# 数据过滤
keep <- filterByExpr(y)
y <- y[keep, , keep.lib.sizes=FALSE]filterByExpr函数用于确定哪些基因的计数足够大，可以保留在统计分析中。也就是我们在介绍原理(基因芯片数据分析（五）：edgeR包的基本原理)中提到的去除表达值为0的基因，实际分析中不是0，我个人理解，counts为只为个位数的也认为是不表达的，所以这里默认的min.count = 10, min.total.count = 15，当然也可以通过参数自己设定。
filterByExpr函数返回的是一个逻辑值类型数据，所以我们需要通过索引获取过滤后的数据。即去除FALSE的基因。
# 计算标准化因子
y <- calcNormFactors(y)
# 查看标准化因子
y$samples
# 计算离散度
y <- estimateDisp(y)
# 显著性检验
et <- exactTest(y)
# 获取排名靠前的基因，这里设置n=100000是为了输出所以基因
et <- topTags(et, n=100000)
# 转换为数据框类型
et <- as.data.frame(et)
# 将行名粘贴为数据框的第一列
et <- cbind(rownames(et),et)
# 指定列名
colnames(et) <- c("gene_id", "log2FoldChange", "log2CPM", "PValue", "FDR")
# 保存数据到本地
write.table(et, "case-vs-control-all.gene.xls", sep = "\t", col.names = TRUE,
row.names = FALSE, quote = FALSE, na = "")
差异基因筛选
# 差异基因筛选
etSig <- et[which(et$PValue < 0.05 & abs(et$log2FoldChange) > 1),]
# 加入一列，up_down 体现上下调信息
etSig[which(etSig$log2FoldChange > 0), "up_down"] <- "Up"
etSig[which(etSig$log2FoldChange < 0), "up_down"] <- "Down"
# 保存文件
write.table(etSig, "case-vs-control-diff-pval-0.05-FC-2-edgeR.gene.xls", sep = "\t",
col.names = TRUE, row.names = FALSE, quote = FALSE, na = "")
这样我们就可以获取差异表达的基因了，可以用于后面的作图，比如火山图！
案例文件和代码，后台回复：chip-data-edgeR，获取！
本文分享自 MedBioInfoCloud 微信公众号，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ，欢迎热爱写作的你一起参与！
腾讯云开发者
