2019-12-15小白新手笔记 GEO下载数据及初步处理可视化

今天介绍另外一种方式 具体哪种好用，我还在检测中

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