vlambda博客
学习文章列表

setwd('~\\下载数据\\')load(file = "mRNA_exprSet.Rda")###metadata <- data.frame(names(mRNA_exprSet)[-1])for (i in 1:length(metadata[,1])) { num <- as.numeric(substring(metadata[i,1],14,15)) if (num %in% seq(1,9)) {metadata[i,2] <- "tumor"} if (num %in% seq(10,29)) {metadata[i,2] <- "normal"}}

mycounts <- mRNA_exprSetmycounts=as.matrix(mycounts)  ##DGEList对象要求矩阵格式rownames(mycounts)=mycounts[,1]exp=mycounts[,2:ncol(mycounts)]dimnames=list(rownames(exp),colnames(exp))data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)data=avereps(data) ##如果有重名的基因，就自动取平均值，方便后面的分析data=data[rowMeans(data)>1,]  ##这个基因在所有的样品上基本都不表达或者是低表达，就删除。

group <- as.vector(metadata$sample) #按照癌症和正常样品数目修改design <- model.matrix(~group) ##变成一个实验设计对象y <- DGEList(counts=data,group=group)  ##变成edgr包可以识别的数据结构，data需要是矩阵格式

# cpm normalizationcountsPerMillion <- cpm(y)countCheck <- countsPerMillion > 1keep <- which(rowSums(countCheck) >= 2)degs.keep <- y[keep,]dim(degs.keep) ##16889 182

y <- calcNormFactors(y) ##对因子进行校正

y <- estimateCommonDisp(y) ##为所有基因都计算同样的dispersiony <- estimateTrendedDisp(y) ##根据不同基因的均值--方差之间的关系来计算dispersion，并拟合一个trended modely <- estimateTagwiseDisp(y) ##计算每个基因的dispersion，并利用经验性贝叶斯方法（empirical bayes）压缩至Trend Didspersion。##个人认为这一项相当于GLM中每个基因的beta值plotBCV(y)

et <- exactTest(y,pair = c("normal","tumor")) #检验两组时间的变异系数## ##类似配对T检验ordered_tags <- topTags(et, n=100000)##P值的校正 采用BH方法

allDiff=ordered_tags$table  allDiff=allDiff[is.na(allDiff$FDR)==FALSE,]diff=allDiff  ##提取所有基因差异分析后的结果，见下表1newData=y$pseudo.counts  ##提取校准后表达的数据，见下表2

###定义两个筛选条件foldChange=1 #差异在2倍以上padj=0.05 ##校正后的P值小于0.05###筛选P值小于0.05且logFC的绝对值大于1的基因diffSig = diff[(diff$FDR < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]write.table(diffSig, file="diffSig.xls",sep="\t",quote=F)## 662 4###显示上调基因###筛选P值小于0.05且logFC的大于1的基因diffUp = diff[(diff$FDR < padj & (diff$logFC>foldChange)),]write.table(diffUp, file="up.xls",sep="\t",quote=F)## 36 4##显示下调基因###筛选P值小于0.05且logFC的小于-1的基因diffDown = diff[(diff$FDR < padj & (diff$logFC<(-foldChange))),]write.table(diffDown, file="down.xls",sep="\t",quote=F)## 626 4###输出所有基因校正后的表达normalizeExp=rbind(id=colnames(newData),newData)write.table(normalizeExp,file="normalizeExp.txt",sep="\t",quote=F,col.names=F)#输出差异基因校正后的表达diffExp=rbind(id=colnames(newData),newData[rownames(diffSig),])write.table(diffExp,file="diffmRNAExp.txt",sep="\t",quote=F,col.names=F) heatmapData <- newData[rownames(diffSig),]

pdf(file="vol.pdf")xMax=max(-log10(allDiff$FDR))+1yMax=12plot(-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()

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()