狮山草堂——猪の家分享 http://blog.sciencenet.cn/u/zhumengjin 行为上保守,思想上激进

博文

【R高级教程】专题二:差异表达基因的分析

已有 42630 次阅读 2011-1-12 11:34 |个人分类:科学研究|系统分类:科研笔记| 表达谱, 高级教程, 基因芯片, 差异表达基因

       

        应学生及个别博友的要求,尽管专业博文点击率和反应均很差,但在去San Diego参加PAG会议之前,还是抽时间给出【R高级教程】的第二专题。专题一给出了聚类分析的示例,本专题主要谈在表达谱芯片分析中如何利用Bioconductor鉴定差异表达基因。

        鉴定差异表达基因是表达谱芯片分析pipeline中必须的分析步骤。差异表达基因分析是根据表型协变量(分类变量)鉴定组间差异表达,它属于监督性分类的一种。在鉴定差异表达基因以前,一般需要对表达值实施非特异性过滤(在机器学习框架下属于非监督性分类),因为适当的非特异性过滤可以提高差异表达基因的检出率、甚至是功效。R分析差异表达基因的library有很多,但目前运用最广泛的Bioconductor包是limma。

        本专题示例依然来自GEO数据库中检索号为GSE11787 的Affymetrix芯片的数据,数据介绍参阅专题一。

>library(limma)
>design <- model.matrix(~ -1+factor(c(1,1,1, 2,2,2)))

        这个是根据芯片试验设计,对表型协变量的水平进行design,比如本例中共有6张芯片,前3张为control对照组,后3张芯片为实验处理组,用1表示对照组,用2表示处理组。其他试验设计同理,比如2*2的因子设计试验,如果每个水平技术重复3次,那么可以表示为:design <- model.matrix(~ -1+factor(c(1,1,1, 2,2,2, 3,3,3, 4,4,4)))。接上面的程序语句继续:

>colnames(design) <- c("control", "LPS")
>fit <- lmFit(eset2, design)
>contrast.matrix <- makeContrasts(control-LPS, levels=design)
>fit <- eBayes(fit)
>fit2 <- contrasts.fit(fit, contrast.matrix)
>fit2 <- eBayes(fit2)
>results<-decideTests(fit2, method="global", adjust.method="BH", p.value=0.01, lfc=1.5)
>summary(results)

>vennCounts(results)
>vennDiagram(results)

        比较遗憾的是,目前limma自带的venn作图函数不能做超过3维的高维venn图,只能画出3个圆圈的venn图,即只能同时对三个coef进行venn作图。上面的venn图只有一个coef,太简单了。下面是一个由本实验室芯片数据得出的三个coef的venn图例:

>heatDiagram(results,fit2$coef)

        红色为control中(与LPS相比)的高表达基因,绿色为control中(与LPS相比)的低表达基因,x轴的数字表示差异表达基因在eset2中所处的位置。

>x<-topTable(fit2, coef=1, number=10000, adjust.method="BH", sort.by="B", resort.by="M")
>write.table(x, file="limma.xls", row.names=F, sep="t")

       将结果写入limma.xls文件中,内容包括AveExpr值(比较组间绝对值的平均差异值)、logFC值(差异倍数)、t值、P值、q值(即adj.P.Val值)和B值。一般logFC值、P值、q值和AveExpr值用来作为差异表达的判断标准,比如差异倍数在2倍以上、绝对差异表达在10以上、P值小于0.01等。在Excel文件中,根据各项判断标准排序,可以很容易地得到差异表达基因列表,这个列表可以用来进行后续的分析,如GO注释、基因网络绘制等。

        专题一中提到实际研究中,一般只用差异表达基因进行聚类分析,在R中,根据差异表达结果过滤表达值很简单(具体的值可以依据芯片数据的实际情况设定,比如P值可以设宽松点0.05、logFC的绝对值也可设为1或2、绝对表达差异也可设低一点,如6或8这样的值):

>y <- x[x$P.Value < 0.01 & (x$logFC > 1.5 | x$logFC < -1.5& x$AveExpr > 10),]
>length(y$ID)
>eset3<-eset2[y$ID,]

      经过上面P值、表达倍数差异和绝对差异的过滤,eset3中就只包含差异表达基因了,这样eset3可用来进行聚类分析了。

【除Bioconductor的logo,所有图片均由程序运行所得】



http://blog.sciencenet.cn/blog-295006-403640.html

上一篇:【研究生培养的几种模式】
下一篇:美国Grand Canyon之旅

6 王秀萍 雷栗 高建国 毛培宏 柳东阳 gaoshannankai

发表评论 评论 (14 个评论)

数据加载中...
扫一扫,分享此博文

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

GMT+8, 2021-10-25 18:35

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部