【step by step】菜鸟学TCGA（4）-用edgeR做差异表达分析
大家好，工作太忙，太久没有更新了，哎，泪……
有的同学问我要代码，有的发了，后面的还没有发，一个一个发好累啊，大家有建议吗？
感觉某宝的这个课程也不贵，300多，有经济能力的小伙伴可以自己买，学得快些~
这次讲的是用一个叫edgeR的包，在R中使用，做差异表达分析，做了之后可以相应的做火山图和热图。
大家关于代码有什么疑问，可以在R中输入“？+软件包名字”
如侵犯版权，请告诉我，我把帖子撤了： ）
#source("http://bioconductor.org/biocLite.R") #source("https://bioconductor.org/biocLite.R")#如果第一个不行，用第二个#
#biocLite("edgeR")
#install.packages("gplots")
foldChange=1#6,7行表示的是过滤的标准，fold change=1意思是差异是两倍，padj=0.05意思是矫正后P值小于0.05，如果得到的差异基因过多，则增大foldchange，减小padj#
padj=0.05
setwd("C:\\Users\\YunGu\\Desktop\\TCGA\\video\\TCGA07") #设置工作目录，意思是在给自己做的事划定一个地盘，如果输入的数据已经在这个地盘，输出的file想放在这个地盘，后面的文件不用单独在输入路径，如果输入输出的位置与这个路径不一致，那么还需要重新输#
library("edgeR")
rt=read.table("mRNA.symbol.txt",sep="\t",header=T,check.names=F) #读取矩阵文件，这是输入的数据路径，改成自己的文件名#
rt=as.matrix(rt) #转化为矩阵#
rownames(rt)=rt[,1] #第一列作为行名#
exp=rt[,2:ncol(rt)] #第二列到最后一列是表达的数据#
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)#15,16行意思是将带引号的数据转换成数值#
data=avereps(data) #有的基因出现过多行，把出现多行的gene取平均值#
data=data[rowMeans(data)>1,] #去除低表达的数据#
group=c(rep("normal",4),rep("tumor",178)) #定义比较组，按照癌症和正常样品数目修改#
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("normal","tumor"))
topTags(et)
ordered_tags <- topTags(et, n=100000) #显示排名前十万的#
allDiff=ordered_tags$table
allDiff=allDiff[is.na(allDiff$FDR)==FALSE,]
diff=allDiff
newData=y$pseudo.counts
write.table(diff,file="edgerOut.xls",sep="\t",quote=F)#先把所有差异表达，输入这个叫edgerOut的文件#
diffSig = diff[(diff$FDR < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]#筛选有显著差异的#
write.table(diffSig, file="diffSig.xls",sep="\t",quote=F)#输出有显著差异表达的到diffSig这个文件#
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两个文件#
diffDown = diff[(diff$FDR < padj & (diff$logFC<(-foldChange))),]
write.table(diffDown, file="down.xls",sep="\t",quote=F)
normalizeExp=rbind(id=colnames(newData),newData)
write.table(normalizeExp,file="normalizeExp.txt",sep="\t",quote=F,col.names=F) #输出所有基因校正后的表达值（normalizeExp.txt）
diffExp=rbind(id=colnames(newData),newData[rownames(diffSig),])
write.table(diffExp,file="diffmRNAExp.txt",sep="\t",quote=F,col.names=F) #输出差异基因校正后的表达值（diffmRNAExp.txt），以后可以和其他数据，比如差异基因表达联合分析#
heatmapData <- newData[rownames(diffSig),]#此处使用的是校正后基因的表达量#
#volcano
pdf(file="vol.pdf") #注意输出的火山图的路径#
xMax=max(-log10(allDiff$FDR))+1
yMax=12
plot(-log10(allDiff$FDR), allDiff$logFC, xlab="-log10(FDR)",ylab="logFC",
main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.4)
diffSub=allDiff[allDiff$FDR<padj & allDiff$logFC>foldChange,]
points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="red",cex=0.4)
diffSub=allDiff[allDiff$FDR<padj & allDiff$logFC<(-foldChange),]
points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="green",cex=0.4)
abline(h=0,lty=2,lwd=3)
dev.off()
#heatmap
hmExp=log10(heatmapData+0.001)
library('gplots')
hmMat=as.matrix(hmExp)
pdf(file="heatmap.pdf",width=60,height=90)#注意输出的热图的储存路径#
par(oma=c(10,3,3,7))
heatmap.2(hmMat,col='greenred',trace="none")
dev.off()
这是输出路径中得到的结果
