注册 |登录

一米阳光分享 http://blog.sciencenet.cn/u/feilei1986 一个生信菜鸟分享自己的实战经验

博文

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

已有 23646 次阅读 2017-7-7 10:37 |个人分类:数据分析|系统分类:科研笔记| Geo, 基因芯片, Agilent, limma, 表达芯片

【背景知识】了解下数据格式
      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] 纯文本查看 复制代码
?
1
2
3
source("http://www.bioconductor.org/biocLite.R")
biocLite("GEOquery")
library(GEOquery)

下载GSE返回的对象--直接根据GSE号下载
[AppleScript] 纯文本查看 复制代码
?
1
2
3
# 下载基因芯片数据，destdir参数指定下载到本地的地址
gse382<-getGEO('GSE49382',destdir =".")##根据GSE号来下载数据，下载_series_matrix.txt.gz
gpl515<-getGEO('GPL17515',destdir =".")  ##根据GPL号下载的是芯片设计的信息, soft文件
[AppleScript] 纯文本查看 复制代码
?
1
2
3
4
###已经下载好的数据从此开始##########
# 打开已下载的本地数据
gse382<-getGEO(filename ='GSE49382_series_matrix.txt.gz')
gpl515<-getGEO(filename ='GPL17515.soft')

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

自己做好三个数据矩阵（表达矩阵，分组矩阵，差异比较矩阵），然后limma的三个步骤（lmFit,eBayes,topTable）就可以啦
[AppleScript] 纯文本查看 复制代码
?
01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
# 查看列名
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的差异表达比较结果

[AppleScript] 纯文本查看 复制代码
?
1
2
# 查看差异表达结果分组情况
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




https://blog.sciencenet.cn/blog-1160733-1065073.html

上一篇：R中批量读入文件--assign+paste函数
收藏 IP: 124.205.77.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (1 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-12-26 22:51

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部