我们接下来通过一个案例介绍利用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 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!