搜索
查看: 4846|回复: 0

表达芯片数据分析实战三:GPL6480//Agilent//ArrayExpress

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2016-9-5 22:08:32 | 显示全部楼层 |阅读模式
选择这个例子来做表达芯片数据分析实战的原因很多,首先它是agilent芯片了,不再是我们常见的affy系列,其次,这个数据也不是在GEO数据库,而是在ArrayExpress数据库。同样的,进入该项目主页可以看到关于该项目的几乎所有资料:http://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-3017/
但很奇怪的是什么没有给出做这个芯片数据的文章是什么: Understanding molecular basis involved in obesity is a crucial first step in developing therapeutic strategies against excess in body weight gain Objective: The purpose of our study was to assess if obesity possess unique peripheral blood gene expression profiles, in a pilot study.
但是里面也说清楚了是如何分组的,下载数据的同时也可以下载到metadata信息,也是同样的芯片数据分析流程:
[Python] 纯文本查看 复制代码
## GPL6480//Agilent//ArrayExpress
getwd()
# [url]http://www.bioconductor.org/packages/release/bioc/vignettes/ArrayExpress/inst/doc/ArrayExpress.pdf[/url]
# [url]http://www.bioconductor.org/packages/release/bioc/manuals/Biobase/man/Biobase.pdf[/url]
# [url]http://bioinformatics.oxfordjournals.org/content/25/16/2092.full.pdf[/url]
# [url]https://www.bioconductor.org/packages/3.3/bioc/vignettes/gCMAP/inst/doc/diffExprAnalysis.pdf[/url]
library(ArrayExpress)
## getAE download all raw data for this study
## ArrayExpress just donwload the expressionSet 
rawset = ArrayExpress("E-MTAB-3017")    
# [url]http://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-3017/[/url]
## it's not a AffyBatch object , but a NChannelSet  for this study .
#exprSet=exprs(rawset)
#Investigation description       E-MTAB-3017.idf.txt
#Sample and data relationship    E-MTAB-3017.sdrf.txt
#Raw data (1)                    E-MTAB-3017.raw.1.zip
#Array design                    A-AGIL-28.adf.txt
## A-AGIL-28 - Agilent Whole Human Genome Microarray 4x44K 014850 G4112F (85 cols x 532 rows)
#R ExpressionSet                 E-MTAB-3017.eSet.r

##  for this platform the element names: E, Eb 
exprSet_Eb=assayDataElement(rawset,"Eb")
exprSet_E=assayDataElement(rawset,"E")

library(Biobase)
library(GEOquery)
#Download GPL file, put it in the current directory, and load it:
gpl <- getGEO('GPL6480', destdir=".")
colnames(Table(gpl)) ## [1] 41108    17
head(Table(gpl)[,c(1,6,7)])  ## you need to check this , which column do you need 
write.csv(Table(gpl)[,c(1,6,7)],"GPL6400.csv")


pdata=pData(rawset)
group_list=pdata$Factor.Value.disease.

##I don't know how to create a proper dat , so just stop here !!!

## remove the duplicated probeset 
rmDupID <-function(a=matrix(c(1,1:5,2,2:6,2,3:7),ncol=6)){
  exprSet=a[,-1]
  rowMeans=apply(exprSet,1,function(x) mean(as.numeric(x),na.rm=T))
  a=a[order(rowMeans,decreasing=T),]
  exprSet=a[!duplicated(a[,1]),]
  #
  exprSet=exprSet[!is.na(exprSet[,1]),]
  rownames(exprSet)=exprSet[,1]
  exprSet=exprSet[,-1]
  return(exprSet)
}
exprSet=rmDupID(dat)

rn=rownames(exprSet)
exprSet=apply(exprSet,2,as.numeric)
rownames(exprSet)=rn
#exprSet=log(exprSet)
#boxplot(exprSet,las=2)
library(limma) 
design=model.matrix(~factor(group_list))
fit=lmFit(exprSet,design)
fit=eBayes(fit) 
topTable(fit,coef=2,n=20)
write.csv(topTable(fit,coef=2,n=Inf,adjust='BH'),"obesityVSnormal.DEG.csv")



芯片数据分析的套路很简单,就是下载芯片数据文件,用R包处理成表达矩阵,然后根据分组信息来做差异分析,最后选取合适的显著性的差异基因做一些功能性分析,如果作者有做其它数据,就结合其它数据一起分析。
也可以对表达矩阵做其它的统计学分析来找有意义的功能改变信息,比如GSEA




上一篇:表达芯片数据分析实战二:GPL15207-nash疾病
下一篇:表达芯片数据分析实战四:配对样本差异分析
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2019-3-26 19:01 , Processed in 0.033989 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.