#TCGA#利用R包edgeR做miRNA表达谱的差异分析

大家好，欢迎来到IT知识分享网。

• 之前我们通过perl得到表达矩阵，现在利用R做差异分析（只截取部分数据）
  
• 得到的火山图
  
• 得到的热图
  

R代码

#设置好工作路径，并卡好p值和foldChange值
setwd("E:/text")
foldChange <- 1
padj <- 0.05

#加载edgeR包
library(edgeR)
#读入表达矩阵，注意读入后fc是数据框
fc <- read.delim("test.txt",sep="\t",header = TRUE,check.names = FALSE)
#剔除表达量为零的miRNA
fc <- fc[rowMeans(fc[,-1])>0,]
#将表达值转为矩阵
counts <- data.matrix(fc[,2:11])
row.names(counts) <- fc$Tags
#创建tumor，normal因子向量（注意数量是否匹配）
group <- factor(c(rep("tumor",5),rep("normal",5)))
y <- DGEList(counts,group=group) #构建列表
y <- calcNormFactors(y) #计算样本内标准化因子
y <- estimateCommonDisp(y) #计算普通的离散度
y <- estimateTagwiseDisp(y) #计算基因或miRNA范围内的离散度
et <- exactTest(y,pair=c("tumor","normal")) #进行精确检验
topTags(et) #输出排名靠前的差异miRNA信息
ordered_tags <- topTags(et,n=100000) #将差异信息存入列表

#剔除FDR值为NA的行
allDiff <- ordered_tags$table
diff <- allDiff[is.na(allDiff$FDR)==FALSE,]
newData <- y$pseudo.counts #将y中的counts信息通过管道存入newData

#将差异miRNA的logF、logCPM、PValue、FDR存入表格
diffx <- rbind(id=colnames(diff),diff)
write.table(diffx,file="edgeOut.xls",sep="\t",quote=F,col.names=F)
#将差异miRNA：p值大于规定值、logFC大于规定值的logF、logCPM、PValue、FDR存入表格
diffSig <- diff[(diff$FDR < padj & (diff$logFC > foldChange | diff$logFC < (-foldChange))),]
write.table(diffSig,file="diffSig.xls",sep="\t",quote=F)
#将上调miRNA的logF、logCPM、PValue、FDR存入表格
diffUp <- diff[(diff$FDR < padj & (diff$logFC > foldChange)),]
write.table(diffUp,file="up.xls",sep="\t",quote=F)
#将下调miRNA的logF、logCPM、PValue、FDR存入表格
diffDown <- diff[(diff$FDR < padj & (diff$logFC < (-foldChange))),]
write.table(diffDown,file="down.xls",sep="\t",quote=F)

#将差异miRNA在tumor和normal样本中的表达值存入表格
diffExp <- rbind(id=colnames(newData),newData[rownames(diffSig),])
write.table(diffExp,file="diffmiRNAExp.txt",sep="\t",quote=F,col.names=F)

#火山图
pdf(file="vol.pdf")
xMax <- max(-log10(diff$FDR))+1 #规定X轴范围
yMax <- 12 #规定Y轴范围
#用实心点将所有点在图里表示
plot(-log10(diff$FDR), diff$logFC, xlab="-log10(FDR)",ylab="logFC",
     main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.4)
#将上调的miRNA标红
diffSub <- diff[(diff$FDR < padj & diff$logFC > foldChange),]
points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="red",cex=0.4)
#将下调的miRNA标绿
diffSub <- diff[(diff$FDR < padj & diff$logFC < (-foldChange)),]
points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="green",cex=0.4)
#加入一条中线
abline(h=0,lty=2,lwd=3)
dev.off() #关闭图形设备

#热图
library(gplots)
heatmapData <- newData[rownames(diffSig),]
hmExp <- log10(heatmapData+0.001)
pdf(file="heatmap.pdf",width=60,height=30)
par(oma=c(10,3,3,7)) #规定图形边界
heatmap.2(hmExp,col='greenred',key=T,keysize=0.8,cexCol=0.8,cexRow=0.5,trace="none")
dev.off()