搜索
查看: 12057|回复: 11

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

[复制链接]

2

主题

5

帖子

113

积分

注册会员

Rank: 2

积分
113
发表于 2017-7-7 10:02:18 | 显示全部楼层 |阅读模式
本帖最后由 feilei_1986 于 2017-7-7 10:39 编辑

【背景知识】了解下数据格式
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


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:设置shiny的DT表格的每一列居中显示
下一篇:shortread 包,主要用来读取fastq文件(原始数据)
回复

使用道具 举报

2

主题

5

帖子

113

积分

注册会员

Rank: 2

积分
113
 楼主| 发表于 2017-7-13 13:50:14 | 显示全部楼层
尚目目 发表于 2017-7-8 17:38
其实这4步,第一步执行完了之后你直接看下gse382,应该是个list。一般gse382[[1]]就是你第3步得到的gse38 ...

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

使用道具 举报

2

主题

52

帖子

474

积分

中级会员

Rank: 3Rank: 3

积分
474
发表于 2017-7-7 12:25:21 | 显示全部楼层
下载 series_matrix 不查看下是否 normalize 了么?我每次习惯性检查下数据是否经过预处理。最后发现质量参差不齐,还是决定尽量直接下RAW自己做预处理算了。
回复 支持 反对

使用道具 举报

2

主题

5

帖子

113

积分

注册会员

Rank: 2

积分
113
 楼主| 发表于 2017-7-8 16:51:30 | 显示全部楼层
尚目目 发表于 2017-7-7 12:25
下载 series_matrix 不查看下是否 normalize 了么?我每次习惯性检查下数据是否经过预处理。最后发现质量参 ...

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

使用道具 举报

2

主题

52

帖子

474

积分

中级会员

Rank: 3Rank: 3

积分
474
发表于 2017-7-8 17:38:01 | 显示全部楼层
本帖最后由 尚目目 于 2017-7-8 17:39 编辑
................................
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步这样再读入一次,以免浪费时间。
回复 支持 反对

使用道具 举报

4

主题

48

帖子

778

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
778
发表于 2017-7-13 18:03:27 | 显示全部楼层
feilei_1986 发表于 2017-7-13 13:50
恩,是的,前两行和后两行都能读取数据,我写后两行的目的是,我分两天做的数据分析,头一天先下载数据, ...

如果你下载的数据 不是标准化过的  数值有正有负   该选用何种方法进行标准化  quantile?百分位数法?
回复 支持 反对

使用道具 举报

2

主题

52

帖子

474

积分

中级会员

Rank: 3Rank: 3

积分
474
发表于 2017-7-14 10:45:19 | 显示全部楼层
feilei_1986 发表于 2017-7-13 13:50
恩,是的,前两行和后两行都能读取数据,我写后两行的目的是,我分两天做的数据分析,头一天先下载数据, ...

我都读入之后存成rda,这个之后再读入速度明显会快一些,读文件感觉太慢了。
回复 支持 反对

使用道具 举报

13

主题

28

帖子

218

积分

中级会员

Rank: 3Rank: 3

积分
218
发表于 2017-7-15 22:43:43 | 显示全部楼层
为什么要用最大表达量,一般不是平均值吗
回复 支持 反对

使用道具 举报

4

主题

48

帖子

778

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
778
发表于 2017-7-22 17:34:37 | 显示全部楼层
尚目目 发表于 2017-7-14 10:45
我都读入之后存成rda,这个之后再读入速度明显会快一些,读文件感觉太慢了。 ...

agilent用什么normalization方法?
回复 支持 反对

使用道具 举报

4

主题

48

帖子

778

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
778
发表于 2017-7-22 17:35:58 | 显示全部楼层
咔吧耀 发表于 2017-7-15 22:43
为什么要用最大表达量,一般不是平均值吗

2种都可以  有相关文献  可以自己去看看
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-22 00:19 , Processed in 0.037998 second(s), 29 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.