生物信息学入门 使用 RNAseq counts数据进行差异表达分析（DEG）——edgeR 算法 数据 代码 结果解读

       差异表达分析通常作为根据基因表达矩阵进行生物信息学分析的第一步，有助于我们观察基因在不同样本中的表达差异，从而确定要研究的基因和表型之间的联系。常用的基因表达数据来自基因芯片或高通量测序。虽然矩阵看起来差不多，但是由于服从不同的分布，因此在进行差异表达的时候需要用不同的方法。对于一般的生命科学领域科研人员来说，了解晦涩的算法并没有太大价值。本文力求精简，从数据——算法——结果三个方面给出最简单的示范。注意：文中代码仅适用于RNAseq的counts数据！使用的是edgeR算法！

1.数据准备

数据准备包括表达矩阵和分组矩阵。

表达矩阵：

分组矩阵：

第一列为样本名称，第二列为组名称，注意每一列都要有列名

2. 使用edgeR包进行差异分析

首先要安装edgeR包和gplots包

source("http://bioconductor.org/biocLite.R")
biocLite("edgeR")
biocLite("gplots")

读取数据

library("edgeR")
library('gplots')
setwd("C:/Users/lenovo/Desktop/sample")
foldChange=1 #fold change=1意思是差异是两倍
padj=0.05#padj=0.05意思是矫正后P值小于0.05
rt=read.csv("fpkm.csv",header=TRUE,row.names=1,check.names = FALSE)  
#读取矩阵文件，这是输入的数据路径，改成自己的文件名#
exp=as.matrix(rt) #转化为矩阵#
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)#15,16行意思是将带引号的数据转换成数值#
data=data[rowMeans(data)>1,] #去除低表达的数据#

读取分组矩阵

group <- read.csv("datTraits.csv",header=TRUE,row.names=1,check.names = FALSE)
group <- group[,1] #定义比较组，按照癌症和正常样品数目修改#
design <- model.matrix(~group) #把group设置成一个model matrix#

计算步骤

y <- DGEList(counts=data,group=group) #group哪些是正常，哪些是癌症样本，让edgeR可以识别#
y <- calcNormFactors(y) #对因子矫正#
y <- estimateCommonDisp(y)#25,26估计变异系数，即估计方差；估计内部差异程度，看组间差异是否比内部差异大，如果大，可选为差异基因#
y <- estimateTagwiseDisp(y)
et <- exactTest(y,pair = c("healthy","T2D"))
topTags(et) #预览结果
summary(de <- decideTestsDGE(et))  #结果的统计信息，基于FDR
ordered_tags <- topTags(et, n=100000)
allDiff=ordered_tags$table
allDiff=allDiff[is.na(allDiff$FDR)==FALSE,]
diff=allDiff
newData=y$pseudo.counts

输出结果

write.csv(diff, "edgerOut.csv")
diffSig = diff[(diff$FDR < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]#筛选有显著差异的#
#write.table(diffSig, file="diffSig.xls",sep="\t",quote=F)#输出有显著差异表达的到diffSig这个文件#
write.csv(diffSig, "diffSig.csv")
diffUp = diff[(diff$FDR < padj & (diff$logFC>foldChange)),]#foldchange>0是上调，foldchange<0是下调#
#write.table(diffUp, file="up.xls",sep="\t",quote=F)#39-42把上调和下调分别输入up和down两个文件#
write.csv(diffUp, "diffUp.csv")
diffDown = diff[(diff$FDR < padj & (diff$logFC<(-foldChange))),]
#write.table(diffDown, file="down.xls",sep="\t",quote=F)
write.csv(diffDown, "diffDown.csv")

差异表达矩阵制作教程：https://blog.csdn.net/tuanzide5233/article/details/83659768

差异表达的热图绘制详见：https://blog.csdn.net/tuanzide5233/article/details/83659501

使用limma包对基因芯片数据进行差异表达分析教程：https://blog.csdn.net/tuanzide5233/article/details/83541443

GEO芯片数据差异表达分析时需要log2处理的原因：https://blog.csdn.net/tuanzide5233/article/details/88542805

GEO芯片数据差异表达分析时是否需要log2以及标准化的问题：https://blog.csdn.net/tuanzide5233/article/details/88542558