搜索
查看: 455|回复: 1

[mRNA-seq] GEO芯片数据分析完整代码

[复制链接]

6

主题

6

帖子

105

积分

注册会员

Rank: 2

积分
105
发表于 2019-4-17 11:44:38 | 显示全部楼层 |阅读模式
本帖最后由 dasunjie6 于 2019-4-17 11:59 编辑

rm(list = ls())

#CEL文件下载
library(GEOquery)
CEL_data=getGEOSuppFiles("GSE36807",baseDir=".\\")
setwd("GSE36807/")
untar("GSE36807_RAW.tar") #解压缩文件
files <-dir(pattern="gz$") #查看满足条件文件
sapply(files, gunzip)
Cellfiles<-dir(pattern="CEL$")

#CEL文件读取
if (!require(affy)) {
  source("http://bioconductor.org/biocLite.R")
  biocLite("affy")
  library("affy")
}
raw_data<-ReadAffy(filenames=Cellfiles) #读取cel文件
plName= cdfName(raw_data) #得到平台的类型
probeName= geneNames(raw_data) #得到所有探针

#CEL文件的boxplot
if (!require(RColorBrewer)) {
  install.packages("RColorBrewer")
  library("RColorBrewer")
  #调色板控制包
}
cols <-brewer.pal(8, "Set1")
par(mfrow= c(1,1))
boxplot(raw_data,col= cols,
        main = 'original cell',las=2)

#CEL文件的Histograms
pdf("log_raw_data_hist.pdf")
par(mfrow= c(3,3))
for (i in 1:length(Cellfiles))
{
  hist(log2(intensity(raw_data[,Cellfiles])),breaks = 1000,xlab = "",main =Cellfiles)
}
dev.off()

#Quality control
install.packages('simpleaffy')
library("simpleaffy")

Data.qc<-qc(raw_data) #获取质量分析报告
plot(Data.qc) #图形化显示报告

library("affyPLM")
PLM_raw_data= fitPLM(raw_data) #对数据集做回归计算
par(mfrow= c(1,1))
Mbox(PLM_raw_data,ylim= c(-1,1),col = cols,main= 'RLE',las= 2)
#绘制RLE(相对对数表达)箱线图
#探针组在某个样品的表达值除以该探针组在所有样品中表达值的中位数后区对数
boxplot(PLM_raw_data,ylim= c(0.92,1.22),col = cols,main= 'NUSE',las= 2)
#绘制NUSE(相对标准差)箱线图,较RLE精确
#探针组在某个样品表达值的标准差除以该探针组在各样品中表达值标准差的中位数

#Data normalization
system.time(raw_data_rma<-rma(raw_data))
system.time(raw_data_mas5 <-mas5(raw_data))
exp_data_rma= exprs(raw_data_rma)
exp_data_mas5 = exprs(raw_data_mas5)
hist(exp_data_mas5,main = "mas5")
hist(log2(exp_data_mas5),main = "log2-mas5")
hist(exp_data_rma,main = "rma")
boxplot(exp_data_rma,lwd= 1,col = cols,main= "rma",las=2)
boxplot(exp_data_rma,lwd= 1,col = cols,main= "mas5",las=2)

par(mfrow= c(3,3))
for (i in 1:8)
{hist(((exp_data_rma[,i])),breaks = 1000,xlab = "",col = 'red',main=colnames(exp_data_rma) )}
par(mfrow= c(3,3))
for (i in 1:8)
{hist(((log2(exp_data_mas5[,i]))),breaks = 1000,xlab = "",col = 'red',main=colnames(exp_data_rma) )}

plot(exp_data_rma[,1],exp_data_rma[,2],xlab= colnames(exp_data_rma)[1],ylab=colnames(exp_data_rma)[2],main = "rma")
plot(exp_data_mas5[,1],exp_data_mas5[,2],xlab= colnames(exp_data_mas5)[1],ylab=colnames(exp_data_mas5)[2],main = "mas5")


#matrix文件下载
library(GEOquery)
getwd()
setwd('../')
gset<-getGEO("GSE36807",destdir="./",getGPL = FALSE) # get GSE from GEO
length(gset)
gset<-gset[[1]]
eset<-exprs(gset) #get the matirxof GSE
gsmName= colnames(eset)
samples = pData(gset)[,1:2]

#Data grouping
datasetHC= eset[,grep("Healthy control",samples[,1])]; nHC= ncol(datasetHC)
datasetCD= eset[,grep("Crohn's Disease",samples[,1])]; nCD= ncol(datasetCD)
datasetUC= eset[,grep("Ulcerative colitis",samples[,1])]; nUC= ncol(datasetUC)
par(mfrow= c(1,4))
nlabels= c("Healthy control","Crohn'sDisease","Ulcerativecolitis")
pie(c(nHC,nCD,nUC),main="rainbow",col=rainbow(3),border=rainbow(3),labels=nlabels)
pie(c(nHC,nCD,nUC),main="cm.colors",col=cm.colors(3),border=cm.colors(3),labels=nlabels)
pie(c(nHC,nCD,nUC),main="topo.colors",col=topo.colors(3),border=topo.colors(3),labels=nlabels)
pie(c(nHC,nCD,nUC),main="heat.colors",col=heat.colors(3),border=heat.colors(3),labels=nlabels)

#Differentially Expressed Genes
#Fold changes
#{10, 20, 15, 25, 15} versus{18, 27, 15, 10, 15}
#{100, 120, 134, 90, 105} versus {40, 27, 32, 38, 41}
#Are they differentially expressed?
#  ave1 = 17 versus ave2 = 17 => no change
#  ave1 = 110 versus ave2 = 36 => 4-fold change

#T-test还是Wilcoxon test
#T-test:观察数据分布,要比较的两组数据是否有正态分布;
#        两组数据的方差是否一致;
#Wilcoxon test:无需任何假设,是非参数的统计检验
#二者抉择:t test有更强的power,但是前提是数据的分布必须符合使用它的假设

#Differentially Expressed Genes
p_values <- c()
p_values_ttest <- c()
FC_values <- c()
for (i in 1:nrow(datasetHC))
{
  p_values <-wilcox.test(datasetHC[i,],datasetCD[i,],paired = F)$p.value
  p_values_ttest <-t.test(datasetHC[i,],datasetCD[i,],paired = F)$p.value
  FC_values <-(mean(2 ^ datasetCD[i,]) / mean(2 ^ datasetHC[i,]))
}
hist(-log10(p_values),breaks = 10,main = 'P value') #hist p_values
hist(FC_values,xlim = c(0,9),breaks = 1000,main = 'FC') #hist FC_values

#常见FDR的方法
#Bonferroni阈值法:将阈值设置为α/n, α为单个假设检验下常用的
#p-value阈值,如0.05;n为所有假设检验的总个数

#火山图

#FDR之后的火山图
fdr.pValue<-p.adjust(p_values,method="fdr",length(p_values))
fdr.pValue_test<-p.adjust(p_values_ttest,method="fdr",length(p_values_ttest))


#Up and down regulated probes
pos = intersect(which(p_values < 0.01),union(which(FC_values < 0.5),which(FC_values >2)))
down_regulated_probe = probset[intersect(which(p_values < 0.01),which(FC_values < 0.5))]
up_regulated_probe = probset[intersect(which(p_values < 0.01),which(FC_values >2))]
pos = intersect(which(p_values_ttest < 0.01),union(which(FC_values < 0.5),which(FC_values >2)))
down_regulated_probe_ttest = probset[intersect(which(p_values_ttest < 0.01),which(FC_values <0.5))]
up_regulated_probe_ttest = probset[intersect(which(p_values_ttest < 0.01),which(FC_values >2))]

#VennDiagram韦恩图
library(VennDiagram)
venn.diagram(list(Up_ttest=up_regulated_probe_ttest,Up_Rank=up_regulated_probe),
             fill=c("red","green"), filename="VennDiagram_up.tiff")
venn.diagram(list(down_ttest=down_regulated_probe_ttest,Down_rank=down_regulated_probe), fill=c("blue",'yellow'),filename="VennDiagram_down.tiff")
#heatmap() 热图
heatmap(Data[c(up_regulated_probe[1:10],down_regulated_probe[1:10]),],margins = c(10,10),scale ="row",col = colorRampPalette(c("green", "black", "red"))(100))
heatmap(Data, scale = "row",Rowv = NULL, Colv = NULL, margins = c(10,10),col = colorRampPalette(c("green", "black", "red"))(100))

#Probe注释
###########get GPL annotation ###############
gpl <-annotation(gset)
platf <-getGEO(gpl, AnnotGPL=TRUE,destdir=".\\")
platf_GPL570_table <-data.frame(attr(dataTable(platf), "table"))
rownames(platf_GPL570_table) = platf_GPL570_table$ID
down_regulated_gene_symbol =platf_GPL570_table[down_regulated_probe,"Gene.symbol"]
up_regulated_gene_symbol = platf_GPL570_table[up_regulated_probe,"Gene.symbol"]





上一篇:如何使用TMHMM serverv2.0软件
下一篇:EBI上的数据存疑?
回复

使用道具 举报

0

主题

1

帖子

47

积分

新手上路

Rank: 1

积分
47
发表于 2019-4-29 15:14:36 | 显示全部楼层
down_regulated_probe = probset[intersect(which(p_values < 0.01),which(FC_values < 0.5))]这行代码probset没有赋值,会报错
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2019-5-20 23:43 , Processed in 0.030695 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.