服务热线
178 0020 3020
R代码1-GEO数据库下载与差异分析
#从GEO数据库上下载数据并分析差异
#https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE70689
#打开上面网址记录GSE70689和GPL平台的GPL7868
#加载GEOquery包，用于下载GEO数据
library("GEOquery")
#用于差异表达分析
library("limma")
#设置当前工作目录,所下载的数据或者生成的数据都在如下目录下
setwd("D:/geo/GSE70689/")
#@_@
#获得基因芯片的表达数据
gse <- getGEO('GSE70689',destdir="./")
#获得基因芯片平台的数据
gpl <- getGEO('GPL7868', destdir = "./")
#^_^
#exprSet获得表达矩阵
exprSet <- exprs(gse[[1]])
dim(exprSet) #获得exprSet的行数和列数，输出如下
#[1] 24611 24
head(exprSet) #查看exprSet的前6行数据，第一列为基因芯片探针名称，第一行为样本信息GSM...
#上面head结果可以查看到一些表达值，应当注意是否已经通过log2变化，如果值一般不超过20，并且不是整数，或者还有负数，基本上认为是已经经过转化的
#一般GEO上面芯片的数据都经过log2转化，如果没有，则可以通过 exprSet = log2(exprSet+1)来转换
#pData获得样本信息,包括样本处理信息
pData <- pData(gse[[1]])
View(pData) #查看下pData数据
write.csv(pData, file="sampleInfo.csv") #将样品信息保存起来，以后用
#基因芯片探针与常用基因名称之间的对应关系
#不同基因芯片平台的数据不一样，一般至少包括探针名称和基因名称，在本例中第一列ID和第二列SYMBOL就是我们需要的
gpl_anno= Table(gpl)
head(gpl_anno) #查看gpl_anno的前6行
#@_@
#将探针ID和基因的名称Symbol提取出来
ID2Symbol = data.frame("ID"=gpl_anno$ID,"SYMBOL"=gpl_anno$SYMBOL)
#^_^
#基因芯片平台经常会有多个探针对应同一个基因的情况，
#因此在计算差异表达之前需要将多个探针对应的同一个基因的数据进行合并
#合并的模式可以选择最大值、平均值、中间值,建议选择最大值
exprSet_temp = as.data.frame(exprSet) #将矩阵转为列表
exprSet_temp$ID = row.names(exprSet_temp) #将行名称转化为一列名为ID的数值
exprSet_temp = merge(exprSet_temp, ID2Symbol, by.x="ID", by.y="ID", all.x=T) #合并
exprSet_temp = exprSet_temp[,-1] #删除第一列ID
exprSet_temp = aggregate(exprSet_temp, by=list(exprSet_temp[,"SYMBOL"]),max) #根据exprSet_temp中“SYMBOL”例进行分组，以最大值作为最终值
exprSet_temp = exprSet_temp[!is.na(exprSet_temp$SYMBOL) & exprSet_temp$SYMBOL!="",]
row.names(exprSet_temp) = exprSet_temp$SYMBOL
exprSet_temp = exprSet_temp[,2:(ncol(exprSet_temp)-1)]
exprSet = exprSet_temp #重新赋值给exprSet
#将最终的表达矩阵保存起来，以备后面使用
write.csv(exprSet, file = "exprSet.csv")
#获得样本的名称GSM...
sample <- pData$geo_accession
#@_@
#设置组别，即每个样本属于哪个组
#如果样本数户比较少，手动添加分组，根据样本顺序设置分组即可，这个分组名可自定义
#如果样本多，需要代码去分组，这里只是展示
#本例中手动分组如下：
group = c("THC_LPS","Cntr","LPS","THC","LPS","THC","THC_LPS","CBD_LPS","Cntr","Cntr","THC","LPS","CBD","Cntr","THC","LPS","CBD","LPS","LPS","THC_LPS","CBD_LPS","THC_LPS","Cntr","Cntr")
#本例中代码分组如下：
temp = pData$title #获得样本信息"THC_LPS.1" "Cntr.1" "LPS.1" ...
temp2 = data.frame(strsplit(temp,"[.]"))[1,] #分割字符串，“.”很特殊，这里需要加中括号
group = as.character(temp2)
summary(as.factor(group)) #查看分组情况，总共6组，以及每组的数据如下
#CBD CBD_LPS Cntr LPS THC THC_LPS
#2 2 6 6 4 4
#^_^
#构建分组矩阵design, 注意下面代码中sample, group, sample_group
sample_group <- data.frame(sample, group)
design <- model.matrix(~0+factor(sample_group$group))
colnames(design) <- levels(factor(sample_group$group))
rownames(design) <- sample_group$sample
#构建比较矩阵，如果想看CBD_LPS与LPS组之间的基因表达差异，如下输入"CBD_LPS-LPS"，即“组名”-“比较的另一组名”
#@_@
contrast.matrix = makeContrasts(CBD_LPS-LPS,levels = design)
#^_^
fit = lmFit(exprSet,design)
fit2 = contrasts.fit(fit,contrast.matrix)
fit2 = eBayes(fit2)
x = topTable(fit2,coef = 1,n=Inf, adjust.method = "BH", sort.by="P") # x中就保存了基因表达差异信息
head(x) #输出如下
# logFC AveExpr t P.Value adj.P.Val B ID
#Mt2 2.968847 10.655913 14.761116 2.793789e-13 5.226900e-09 19.13249 Mt2
#Vegfa 2.195245 9.826903 11.422124 5.323503e-11 4.460197e-07 14.73048 Vegfa
#Mmp23 1.495654 7.421330 11.252894 7.151954e-11 4.460197e-07 14.47264 Mmp23
#Aqp9 3.271151 9.737254 10.075593 6.101003e-10 2.306821e-06 12.57310 Aqp9
#2210008I11Rik 2.231609 10.656282 10.045244 6.461539e-10 2.306821e-06 12.52159 2210008I11Rik
#2810428C21Rik 2.316433 9.578238 9.973933 7.398004e-10 2.306821e-06 12.40007 2810428C21Rik
write.csv(x,file="CBD_LPS-LPS.csv") #保存结果
#筛选表达值差异大于2倍（log2=1) 且 p<0.05
y = x[abs(x$logFC)>1 & x$P.Value<0.05,]
write.csv(y,file="CBD_LPS-LPS_defference_expression_genes.csv")
#把 94行，109行，113行的名称改下，就可以得到CBD组与Cntr组的差异表达基因
contrast.matrix = makeContrasts(CBD-Cntr,levels = design)
fit = lmFit(exprSet,design)
fit2 = contrasts.fit(fit,contrast.matrix)
fit2 = eBayes(fit2)
x = topTable(fit2,coef = 1,n=Inf, adjust.method = "BH", sort.by="P") # x中就保存了基因表达差异信息
head(x)
write.csv(x,file="CBD-Cntr.csv")
y = x[abs(x$logFC)>1 & x$P.Value<0.05,]
write.csv(y,file="CBD-Cntr_defference_expression_genes.csv")
附件