搜索
查看: 5427|回复: 4

表达芯片数据分析实战四:配对样本差异分析

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2016-9-5 22:25:24 | 显示全部楼层 |阅读模式
本次介绍的项目实战是非常新的数据,看ID号就知道了,在GEO数据库主页也有详细介绍:http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE70529
同样也需要读懂该文章:
  • Magkos F, Fraterrigo G, Yoshino J, Luecking C et al. Effects of Moderate and Subsequent Progressive Weight Loss on Metabolic Function and Adipose Tissue Biology in Humans with Obesity. Cell Metab 2016 Apr 12;23(4):591-601. PMID: 26916363
这样分析的数据结果自己才有把握,也好联系到生物学意义。
The purpose of this study was to evaluate the effect of progressive weight loss (5, 10, 15% weight loss) on metabolic function such as multi-organ insulin sensitivity and beta-cell function in obese people. We conducted microarray analysis to determine the effect of progressive weight loss on adipose tissue gene expression profile.
这个芯片平台([HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array [transcript (gene) version])并非是affy的最新款芯片,但被研究的也比较多了,在R里面已经有了对于的注释包,hugene10sttranscriptcluster.db


也可以做同样的分析,但是要注意,里面的样本是配对,跟前面3个有很大的区别,在做统计学检验的时候尤其要注意到。   
[Python] 纯文本查看 复制代码
## HuGene-1_0-st // GPL6244 // obesity 
## paired or not 
rm(list=ls())
library(GEOquery)
library(limma) 
GSE70529 <- getGEO('GSE70529', destdir=".",getGPL = F)
exprSet=exprs(GSE70529[[1]])
library("annotate")
GSE70529[[1]]

platformDB='hugene10sttranscriptcluster.db'
library(platformDB, character.only=TRUE)
probeset <- featureNames(GSE70529[[1]])
#EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))
SYMBOL <-  lookUp(probeset, platformDB, "SYMBOL")
a=cbind(SYMBOL,exprSet)
## 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(a)
rn=rownames(exprSet)
exprSet=apply(exprSet,2,as.numeric)
rownames(exprSet)=rn
exprSet[1:4,1:4]
#exprSet=log(exprSet)
#boxplot(exprSet,las=2)

pdata=pData(GSE70529[[1]])
## check the metadat and find the correct group information
# we must tell R that group should be interpreted as a factor !
individuals=factor(unlist(lapply(pdata$title,function(x) strsplit(as.character(x),"-")[[1]][1])))
treatment=unlist(lapply(pdata$title,function(x) strsplit(as.character(x),"-")[[1]][2]))
treatment=factor(sub('2','',treatment))

## if only take the treatment into accont 
design=model.matrix(~treatment)
fit=lmFit(exprSet,design)
fit=eBayes(fit) 
#vennDiagram(fit,include=c("up","down"))
write.csv(topTable(fit,coef='treatmentB',n=Inf,adjust='BH'),"BvsA_NonPaired.DEG.csv")
write.csv(topTable(fit,coef='treatmentC',n=Inf,adjust='BH'),"CvsA_NonPaired.DEG.csv")
write.csv(topTable(fit,coef='treatmentD',n=Inf,adjust='BH'),"DvsA_NonPaired.DEG.csv")
BvsA_NonPaired=topTable(fit,coef='treatmentB',n=Inf,adjust='BH')
CvsA_NonPaired=topTable(fit,coef='treatmentC',n=Inf,adjust='BH')
DvsA_NonPaired=topTable(fit,coef='treatmentD',n=Inf,adjust='BH')


## if paired analysis :
design=model.matrix(~individuals+treatment)
fit=lmFit(exprSet,design)
fit=eBayes(fit)  
write.csv(topTable(fit,coef='treatmentB',n=Inf,adjust='BH'),"BvsA_Paired.DEG.csv")
write.csv(topTable(fit,coef='treatmentC',n=Inf,adjust='BH'),"CvsA_Paired.DEG.csv")
write.csv(topTable(fit,coef='treatmentD',n=Inf,adjust='BH'),"DvsA_Paired.DEG.csv")
BvsA_Paired=topTable(fit,coef='treatmentB',n=Inf,adjust='BH')
CvsA_Paired=topTable(fit,coef='treatmentC',n=Inf,adjust='BH')
DvsA_Paired=topTable(fit,coef='treatmentD',n=Inf,adjust='BH')


t1=rownames(DvsA_Paired[DvsA_Paired$P.Value<0.05,])
t2=rownames(DvsA_NonPaired[DvsA_NonPaired$P.Value<0.05,])
length(setdiff(t1,t2))   ## 594
length(setdiff(t2,t1))   ## 16 
length(intersect(t1,t2))  ##1131


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




上一篇:表达芯片数据分析实战三:GPL6480//Agilent//ArrayExpress
下一篇:mRNA-seq数据分析实战
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

4

主题

24

帖子

184

积分

注册会员

Rank: 2

积分
184
发表于 2017-3-4 16:43:29 | 显示全部楼层
学习了,收藏。
回复 支持 反对

使用道具 举报

0

主题

5

帖子

77

积分

注册会员

Rank: 2

积分
77
QQ
发表于 2017-11-27 19:50:45 | 显示全部楼层
大神,你好。不好意思打扰你了。我看了你写的那篇“表达芯片数据分析实战四:配对样本差异分析”的帖子,非常受用,非常感谢您。但是有一个小问题想请教一下你,在运行到fit=lmFit(exprSet,design)这个步骤时,显示error,原因时我的行数和列数不匹配,row dimension of design doesn't match column dimension of data object.
> dim(exprSet)
[1] 20502    10
> dim(design)
[1] 4 4
> 我找了很久也不知道如何解决这个问题,非常着急。请大神帮忙指导一下,万分感谢。
回复 支持 反对

使用道具 举报

0

主题

5

帖子

77

积分

注册会员

Rank: 2

积分
77
QQ
发表于 2017-11-27 20:20:40 | 显示全部楼层
应该是我矩阵的样本数量与design的样本数量不一致,如何解决,还望大神帮忙解决一下,万分感谢
回复 支持 反对

使用道具 举报

0

主题

5

帖子

77

积分

注册会员

Rank: 2

积分
77
QQ
发表于 2017-11-27 22:02:32 | 显示全部楼层
大神,我在您的代码中加入了这段语句
aa <- as.data.frame(unlist(SYMBOL))
colnames(aa) <- "gene"
结果运行出来的差异基因的结果logFC完全不一样
例如
PBX1        -1.1576875        8.61917        -3.856555869        0.010115344        0.19606529        -2.629316855
PBX1        1.143156        8.61917        6.449239545        6.16E-05        0.159769924        1.810007399
请问这是什么原因?应该要加吗?
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-5-27 04:07 , Processed in 0.147340 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.