今天介绍另外一种方式 具体哪种好用，我还在检测中
rm(list = ls())
options(stringsAsFactors = F)
library(GEOquery)
setwd("C:\\Users\\zhouwenqing789\\Desktop\\DFEF")
eset1<-getGEO("GSE83521",destdir = ".",getGPL = F) #下载了平台注释文件以及不是矩阵的形式,list。在电脑里面文件打不开。但是可以读取。
class(eset1)
exp1=exprs(eset1[[1]])#基因的表达矩阵
dim(exp1)
exp1[1:4,1:4]
boxplot(exp1,las=2)#做个厢式图看一下数据的分布情况
gpl1=eset1[[1]]@annotation #获取注释信息
pd1=pData(eset1[[1]])#pData和exprs函数都可以处理这个表达对象，从而分别得到样品描述矩阵（每个样本的信息）和样品表达量矩阵
fd1=fData(eset1[[1]])#基因的信息？？？
write.table(eset1,"eset1.txt",col.names = T,row.names = F,quote = F)
save(eset1,exp1,pd1,gpl1,file = "step1.Rdata")
load("step1.Rdata")#载入数据
下载另外一组数据
est2=getGEO("GSE89143",destdir = ".",getGPL = F)
class(est2)
eset2=est2
exp2=exprs(eset2[[1]])
boxplot(exp2,las=2) #发现数据差异太大
exp2=log2(exp2+1) #取log值
boxplot(exp2,las=2)#发现中位数部在一条线上，需要用函数normalizeBetweenArrays拉回到正常水平。
library(limma)
exp2=normalizeBetweenArrays(exp2)
boxplot(exp2,las=2)
index=sort.int(pd2$characteristics_ch1,index.return = T) #将pd2中某列按样本和正常组分类
pd2=pd2[index$ix,]#根据分类后重新调整pd2的行顺序
exp2=exp2[,match(rownames(pd2),colnames(exp2))] #将exp2这个矩阵列明按照pd2的行名的顺序重新排列
#提取平台中基因ID与基因名的矩阵
if(T){gpl<-getGEO("GPL19978",destdir = ".")}
dim(gpl)
colnames(Table(gpl))
head(Table(gpl))[,c(1,2)]
ids=Table(gpl)[,c(1,2)]
ids=ids[-c(1:15),]
#找出两个矩阵相同的基因ID，并合成一个矩阵
x1=exp1[rownames(exp1)%in%ids$ID,] #提取有ID的行
x2=exp2[rownames(exp2)%in%ids$ID,] #提取有ID的行
cg=intersect(row.names(x1),row.names(x2)) 两个表达矩阵是否有相同的行
length(cg)
x_merge=cbind(x1[cg,],x2[cg,]) #在找到相同的行时给予列合并
dim(x_merge)
#去批次化 将上面的数据boxplot后发现两组两本的中位数不在一条线上，故需要去批次化
group_list<-c(rep("tumor",6),rep("normal",6),rep(c("tumor","normal"),each=3))
gse<-c(rep("gse83527",12),rep("gse89143",6))
table(group_list,gse)
dat<-x_merge
library(sva)
library(limma)
batch=c(rep("gse83527",12),rep("gse89143",6))
design=model.matrix(~factor(group_list))
ex_b_limma<-removeBatchEffect(dat,batch = batch,design = design)##注意ex_b_limma的格式
boxplot(ex_b_limma)
#找差异基因
fit=lmFit(ex_b_limma,design)#一定要注意design的格式 ，第二列必须是因子
fit=eBayes(fit)
options(digits = 4)
topTable(fit,coef = 2,adjust="BH")
