这一讲其实算不上是自学miRNA-seq分析，本质就是affymetrix的mRNA表达芯片数据分析，而且还是最常用的那种GPL570 HG-U133_Plus_2，但是因为是跟miRNA样本配对检测的，而且后面会利用到这两个数据分析结果来做共表达网络分析等等，所以就贴出对该芯片数据的分析结果。文章里面也提到了 Messenger RNA expression analysis identified 731 probe sets with significant differential expression，作者挑选的差异分析结果的显著基因列表如下：## http://journals.plos.org/plosone/article/asset?unique&id=info:doi/10.1371/journal.pone.0108051.s002
## mRNA expression array - GSE60291 (Affymetrix Human Genome U133 Plus 2.0 Array)
hgu133plus2芯片数据太常见了，可以从GEO里面下载该study的原始测序数据，然后用affy,limma包来分析，也可以直接用GEOquery包来下载作者分析好的表达矩阵，然后直接做差异分析。我这里选择的是后者，而且我跟作者分析方法有一点区别是，我先把探针都注释好了基因，然后对每个基因只挑最大表达量的基因。而作者是直接对探针为单位的的表达矩阵进行差异分析，对分析结果里面的探针进行基因注释。我这里无法给出哪种方法好的绝对评价。代码如下：
rm(list=ls())
library(GEOquery)
library(limma)
GSE60291 <- getGEO('GSE60291', destdir=".",getGPL = F)#下面是表达矩阵
exprSet=exprs(GSE60291[[1]])
library("annotate")
GSE60291[[1]]
## 下面是分组信息
pdata=pData(GSE60291[[1]])
treatment=factor(unlist(lapply(pdata$title,function(x) strsplit(as.character(x),"-")[[1]][1])))
#treatment=relevel(treatment,'control')
## 下面做基因注释
platformDB='hgu133plus2.db'
library(platformDB, character.only=TRUE)
probeset <- featureNames(GSE60291[[1]])
#EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))
SYMBOL <- lookUp(probeset, platformDB, "SYMBOL")
## 下面对每个基因挑选最大表达量探针
a=cbind(SYMBOL,exprSet)
## remove the duplicated probeset
rmDupID <-function(a=matrix(c(1,1:5,2,2:6,2,3:7),ncol=6)){
exprSet=a[,-1]
rowMeans=apply(exprSet,1,function(x) mean(as.numeric(x),na.rm=T))
a=a[order(rowMeans,decreasing=T),]
exprSet=a[!duplicated(a[,1]),]
#
exprSet=exprSet[!is.na(exprSet[,1]),]
rownames(exprSet)=exprSet[,1]
exprSet=exprSet[,-1]
return(exprSet)
}
exprSet=rmDupID(a)
rn=rownames(exprSet)
exprSet=apply(exprSet,2,as.numeric)
rownames(exprSet)=rn
exprSet[1:4,1:4]
#exprSet=log(exprSet) ## based on e
boxplot(exprSet,las=2)
## 下面用limma包来进行芯片数据差异分析
design=model.matrix(~ treatment)
fit=lmFit(exprSet,design)
fit=eBayes(fit)
#vennDiagram(decideTests(fit))
DEG=topTable(fit,coef=2,n=Inf,adjust='BH')
dim(DEG[abs(DEG[,1])>1.2 & DEG[,5]<0.05,]) ## 806 genes
write.csv(DEG,"ET1-normal.DEG.csv")
得到的ET1-normal.DEG.csv 文件就是我们的差异分析结果，可以跟文章提供的差异结果做比较，是几乎一模一样的！
如果根据logFC 1.2 p 矫正P 值0.05来挑选，可以拿到806个基因。