公共数据库没有错误吗？
2019-11-21 本文已影响0人
mayoneday
我们对GSE17215进行分析
按照原数据3X3分组
image.png
1.注释基因，整理数据
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)#在调用as.data.frame的时，将stringsAsFactors设置为FALSE可以避免character类型自动转化为factor类型
# 注意查看下载文件的大小，检查数据
f='GSE17215_eSet.Rdata'
# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE12452
library(GEOquery)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset <- getGEO('GSE17215', destdir=".",
AnnotGPL = F, ## 注释文件
getGPL = F) ## 平台文件
save(gset,file=f) ## 保存到本地
}
load('GSE17215_eSet.Rdata') ## 载入数据
class(gset) #查看数据类型
length(gset) #
class(gset[[1]])
gset
# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list
a=gset[[1]] #
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
dat<-log2(dat+1)
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
boxplot(dat,las=2)
pd=pData(a) #通过查看说明书知道取对象a里的临床信息用pData
## 挑选一些感兴趣的临床表型。
pd$characteristics_ch1.3
group_list= ifelse(grepl("agent: paclitaxel",pd$characteristics_ch1.3),"paclitaxel",'salinomycin')
table(group_list)
dat[1:4,1:4]
if(F){
library(GEOquery)
#Download GPL file, put it in the current directory, and load it:
gpl <- getGEO('GPL3921', destdir=".")
colnames(Table(gpl))
head(Table(gpl)[,c(1,15)]) ## you need to check this , which column do you need
probe2gene=Table(gpl)[,c(1,15)]
head(probe2gene)
library(stringr)
save(probe2gene,file='probe2gene.Rdata')
}
#
load(file='probe2gene.Rdata')
ids=probe2gene
head(ids)
colnames(ids)=c('probe_id','symbol')
ids=ids[ids$symbol != '',]
ids=ids[ids$probe_id %in% rownames(dat),]
dat[1:4,1:4]
dat=dat[ids$probe_id,]
ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4] #保留每个基因ID第一次出现的信息
dim(dat)
save(dat,group_list,file = 'step1-output.Rdata')
2.检查数据
发现按照3X3分组样本主成分和聚类，有一个样本表现很奇怪
image.png
image.png
3.差异分析
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')
library(limma)
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
head(design)
exprSet=dat
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts("salinomycin-paclitaxel",
levels = design)
contrast.matrix ##这个矩阵声明，我们要把 Tumor 组跟 Normal 进行差异分析比较
deg = function(exprSet,design,contrast.matrix){
##step1
fit <- lmFit(exprSet,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix)
##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
return(nrDEG)
}
deg = deg(exprSet,design,contrast.matrix)
head(deg)
save(deg,file = 'deg.Rdata')
image.png
重新分组
而且adjP值也没有差异，怀疑此样本分组命名有问题
从新查看了分组，把这个单出来的样本放到其他3个一起，分成4x2两组
通过查看样本排序得出分组是：group_list=c(rep("paclitaxel",2),rep('salinomycin',4))
1.检查分组
可以看到样本主成分和热图可以看到样本有呈现相同特点的趋势了
image.png
image.png
2.差异分析
image.png
这时候可以看adjP值是有意义的了
image.png
3.富集分析
kegg_up_down.png
kegg_up_down_gsea.png
kk.diff.dotplot.png
kk.down.dotplot.png
kk.up.dotplot.png
两次不同分组logFC的散点图（比较两次的结果）
rm(list = ls()) ## 魔幻操作，一键清空~
load(file = 'deg.Rdata')
load(file = 'degxiuzheng.Rdata')
degxiu<-deg
a<-data.frame(xiu<-degxiu$logFC,yuan<-degyuan$logFC)
a$xiu<-degxiu$logFC
a$yuan<-degyuan$logFC
# 基函数：colour设置分组
a$group<-degxiu$logFC#随便赋值生成一个组，方便下面代码运行
#根据一致性进行分组
for(i in 1:nrow(a)){
a[i,"group"]<-ifelse(abs(a[i,1]-a[i,2])<1,"yizhi","buzhiyi")
}#绝对值小于1就认为他是一致的
table(a$group)
ggplot(a, aes(x = xiu, y = yuan, colour = group)) +
# 散点图函数
geom_point()
image.png
把logFC相差絕對值小于1的认为是一致画图
可以看到，它成斜角的不一致的数据挺多的
