基因表达芯片数据分析——Agilent

feilei_1986

【背景知识】了解下数据格式
GEO数据库基础知识

• GEO Platform (GPL) 芯片平台
• GEO Sample (GSM) 样本ID号
• GEO Series (GSE) study的ID号
  
  

• GEO Dataset (GDS) 数据集的ID号
  

【实例操作】

以蔡南海lab发表的Genome research-2014-Wang文章为例：
ATH NAT array数据上传至NCBI GEO, GSE49382
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE49382


注意这个芯片平台Platforms，这个就是 芯片探针与gene的对应关系。
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL17515



【数据下载】利用R包来进行下载——GEOquery

[AppleScript] 纯文本查看 复制代码

source("http://www.bioconductor.org/biocLite.R")
biocLite("GEOquery")
library(GEOquery)

下载GSE返回的对象--直接根据GSE号下载

[AppleScript] 纯文本查看 复制代码

# 下载基因芯片数据，destdir参数指定下载到本地的地址
gse382 <- getGEO('GSE49382', destdir = ".") ##根据GSE号来下载数据，下载_series_matrix.txt.gz
gpl515 <- getGEO('GPL17515', destdir = ".")  ##根据GPL号下载的是芯片设计的信息, soft文件

[AppleScript] 纯文本查看 复制代码

###已经下载好的数据从此开始##########
# 打开已下载的本地数据
gse382 <- getGEO(filename = 'GSE49382_series_matrix.txt.gz')
gpl515 <- getGEO(filename = 'GPL17515.soft')

【数据分析】用limma进行差异表达分析

自己做好三个数据矩阵（表达矩阵，分组矩阵，差异比较矩阵），然后limma的三个步骤（lmFit,eBayes,topTable）就可以啦

[AppleScript] 纯文本查看 复制代码

# 查看列名
colnames(Table(gpl515))
Table(gpl515)[1:10,1:4] # 前10行前4列信息
# 保存自己想要的信息，在此例中c(1,4)，即第一列和第四列，分别是gse382中的行名ID和其对应的gene name
write.csv(Table(gpl515)[,c(1,4)],"GPL515.csv", row.names = F)

# gse382中的行名ID与gene name的对应关系
genename = read.csv("GPL515.csv")

# 构建表达矩阵
exprSet <- as.data.frame(exprs(gse382)) # 得到表达矩阵，行名为ID，需要转换
# 转换ID为gene name
exprSet$ID = rownames(exprSet)
express = merge(x=exprSet, y=genename, by="ID", all.x = T)
express$ID = NULL
# 去除重复的gene ，保留每个基因最大表达量结果
# 参考：[url=http://www.biotrainee.com/thread-113-1-1.html]http://www.biotrainee.com/thread-113-1-1.html[/url]
rowMeans=apply(express,1,function(x) mean(as.numeric(x),na.rm=T))
express = express[order(rowMeans,decreasing = T),]
express = express[!duplicated(express[,28]),]  #express第28列为gene name

rownames(express) = express[,28]
express=express[,-28]
# 至此，表达矩阵（express）已构建好

# 构建分组矩阵
# 参考：[url=http://www.bio-info-trainee.com/1194.html]http://www.bio-info-trainee.com/1194.html[/url]
pdata = pData(gse382) # 每个sample所对应的信息，包括处理条件等
group_list = subset(pdata, select=title) # Sample的分组信息
group_list$condition = rep(c("c0","h0","r0","c1","h1","r1","c6","h6","r6"), each=3)
design = model.matrix(~0+factor(group_list$condition))
colnames(design) = levels(factor(group_list$condition))
rownames(design) = colnames(express)
# 至此，分组矩阵（design）已构建好

# 构建差异比较矩阵
# 参考：[url=http://manuals.bioinformatics.ucr.edu/home/R_BioCondManual#TOC-Analysis-of-Differentially-Expressed-Genes]http://manuals.bioinformatics.uc ... lly-Expressed-Genes[/url]
contrast.matrix = makeContrasts(c0-c1,c0-c6,h0-h1,h0-h6,r0-r1,r0-r6,levels = design)
# 至此，差异表达矩阵已构建好

fit = lmFit(express,design)
fit2 = contrasts.fit(fit, contrast.matrix)
fit2 = eBayes(fit2)

# 得到两两差异表达的结果
# c0 vs. c1
x = topTable(fit2, coef = 1, n=Inf, adjust.method = "BH", sort.by = "P")
sum(x$adj.P.Val<0.05)
re = x[x$adj.P.Val < 0.05 & (x$logFC > 1 | x$logFC < -1),]  #选取adj.p.value<0.05且|logFC|>1的基因
write.csv(re, "c0-c1_DEG_limma.re.csv",quote = F)
# coef可是column number，也可以是column name，这样就可以指定你所感兴趣的两两比较的结果
# 在此例中coef =1 就是c0-c1的差异表达比较结果

file:///C:/Users/fanyy/AppData/Local/Temp/enhtmlclip/Image.png

[AppleScript] 纯文本查看 复制代码

# 查看差异表达结果分组情况
results = decideTests(fit2,p.value = 0.05)

矩阵中1代表显著上调，-1代表显著下调，0代表无显著差异

【参考】
用GEOquery从GEO数据库下载数据
http://www.bio-info-trainee.com/bioconductor_China/software/GEOquery.html

芯片探针注释基因ID或者symbol，并对每个基因挑选最大表达量探针 | 生信菜鸟团
http://www.bio-info-trainee.com/1502.html

没有必要用R包GEOquery | 生信菜鸟团
http://www.bio-info-trainee.com/1571.html

limma分析参考：
R & Bioconductor - Manuals
http://manuals.bioinformatics.ucr.edu/home/R_BioCondManual#TOC-Analysis-of-Differentially-Expressed-Genes

用limma包对芯片数据做差异分析 | 生信菜鸟团
http://www.bio-info-trainee.com/1194.html

feilei_1986

尚目目 发表于 2017-7-8 17:38
其实这4步，第一步执行完了之后你直接看下gse382，应该是个list。一般gse382[[1]]就是你第3步得到的gse38 ...

恩，是的，前两行和后两行都能读取数据，我写后两行的目的是，我分两天做的数据分析，头一天先下载数据，第二天再分析的时候直接从已经下载的数据直接读取。 如果一次性做完，后两行可以不要的。

尚目目

下载 series_matrix 不查看下是否 normalize 了么？我每次习惯性检查下数据是否经过预处理。最后发现质量参差不齐，还是决定尽量直接下RAW自己做预处理算了。

feilei_1986

尚目目 发表于 2017-7-7 12:25
下载 series_matrix 不查看下是否 normalize 了么？我每次习惯性检查下数据是否经过预处理。最后发现质量参 ...

嗯嗯，看过了，是normalize过后的。应该是用RAW data进行处理的，我还不会…… 先从作者已经上传好的表达矩阵开始做，下一步开始学习处理RAW DATA ~  

尚目目

................................

gse382 <- getGEO('GSE49382', destdir = ".") ##根据GSE号来下载数据，下载_series_matrix.txt.gz
gpl515 <- getGEO('GPL17515', destdir = ".")  ##根据GPL号下载的是芯片设计的信息, soft文件
gse382 <- getGEO(filename = 'GSE49382_series_matrix.txt.gz')
gpl515 <- getGEO(filename = 'GPL17515.soft')

其实这4步，第一步执行完了之后你直接看下gse382，应该是个list。一般gse382[[1]]就是你第3步得到的gse382对象了，同样gpl382也储存在第一步gse382对象里，你不用像3、4步这样再读入一次，以免浪费时间。

渊梦无痕

feilei_1986 发表于 2017-7-13 13:50
恩，是的，前两行和后两行都能读取数据，我写后两行的目的是，我分两天做的数据分析，头一天先下载数据， ...

如果你下载的数据　不是标准化过的　　数值有正有负　　　该选用何种方法进行标准化　　ｑｕａｎｔｉｌｅ？百分位数法？

尚目目

feilei_1986 发表于 2017-7-13 13:50
恩，是的，前两行和后两行都能读取数据，我写后两行的目的是，我分两天做的数据分析，头一天先下载数据， ...

我都读入之后存成rda，这个之后再读入速度明显会快一些，读文件感觉太慢了。

咔吧耀

为什么要用最大表达量，一般不是平均值吗

渊梦无痕

尚目目 发表于 2017-7-14 10:45
我都读入之后存成rda，这个之后再读入速度明显会快一些，读文件感觉太慢了。 ...

agilent用什么normalization方法？

渊梦无痕

咔吧耀 发表于 2017-7-15 22:43
为什么要用最大表达量，一般不是平均值吗

2种都可以  有相关文献  可以自己去看看