bioconductor系列教程之一分析基因芯片下(数据分析) 72

本节目录:
背景

教程数据下载

SAM算法

CyberT算法

Rank products(RP)算法

总结

提取差异表达的基因

从基因芯片当中提取生物学的信息需要合理的统计学方法。人们已经为优化传统统计学方法在基因芯片方面的应用做出了多年的努力。但是直到现在,最主要的努力依然还是依据实验设计的差别,用统计学方法提取出差异表达的基因,然后再转回使用实验的方法去验证这个结果。

在提取差异表达的基因时,人们总是会有这两种考虑,一是不可漏过一个,二是不能错杀过多(在英语里称为false discovery rate(FDR)错误发现率)。常见的手段是使用多种统计学方法来分析同样一个结果,尽可能多的得到差异表达的基因,而排除那些假的信号。然而学习和使用多种统计分析手段并不一定对于每一个生物学工作者都是非常容易的,这需要付出时间和努力。在这里,我们尽量多介绍几种常用的统计分析手段,并给出实践中人们常常使用的组合,来帮助你更好的分析自己的数据。

现在常用的分析手段主要有:significance analysis of microarrays(SAM),CyberT和Rank products(RP)三种手段。其中CyberT是bioconductor当中最为常用的分析手段,因为它的算法完整地被limma库实现。但有研究指出,使用SAM和RP算法相结合可能是最佳的方案。其实任何一种算法都是有局限性的,我们需要从根本上对算法有所了解,然后才能有针对性地选择合适的算法。

SAM:Tusher VG, Tibshirani R, Chu G. Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci USA 2001; 98:5116-21

CyberT: Baldi P, Long AD. A Bayesian framework for the analysis of microarray expression data: regularized t-test and statistical inferences of gene changes. Bioinformatics 2001; 17:509-19

RP: Breitling R, Armengaud P, Amtmann A, et al. Rank products: a simple, yet powerful, new method to detect differentially regulated genes in replicated microarray experiments. FEBS lett 2004;573:83-92

Pages: 1 2 3 4 5 6

72 thoughts on “bioconductor系列教程之一分析基因芯片下(数据分析)

  1. Reply joyce 12月 12,2010 9:06 下午

    O(∩_∩)O~ 看到了 多谢~

  2. Reply joyce 12月 12,2010 10:04 下午

    O(∩_∩)O~ 刚看到后面 恭喜恭喜 祝小宝宝健康快乐的成长~

  3. Reply gao shan 4月 29,2011 11:04 上午

    请问为什么e10 有2个,是双通道芯片
    还是实验重复2次,另外实验矩阵的0和1能给解释一下么

    • Reply admin 4月 29,2011 8:15 下午

      e10, E10是两组不同的因子。实验设计矩阵是比较难懂。1和0都是权重系数,在处理的时候,会用这个来乘以相应的数据,自然,乘以1就说明使用这个数据,乘以0就说明不使用。

  4. Reply gao shan 4月 29,2011 11:12 上午

    能够和我们几个搞生物信息的一起总结一下,写本书,信息量太大了
    大家需要合作才能干点事

  5. Reply mimisikai 8月 3,2011 1:41 下午

    您好 能不能分享一下 “bioconductor全自动处理基因芯片CEL文件R程序”?
    谢谢!

  6. Reply mimisikai 8月 23,2011 1:44 下午

    hi 为什么 我设计的matrix后 lmfit总是显示Error in lm.fit(design, t(M)) : incompatible dimensions?

    • Reply admin 8月 23,2011 2:29 下午

      可否把详细的matrix设计提供一下?只这么一句,我也只能得到字面上的意思。

      • Reply mimisikai 8月 23,2011 4:59 下午

        goodArrayEset<-exprs(gcrma(goodArray))
        treatments<-factor(c(1:5),labels=c("HJ","SW","JMX","CYM-042","CYM-144"));treatments
        contrasts(treatments)<-cbind(pre=c(1,1,0,1,1),d3=c(1,1,0,0,0),d7=c(1,0,1,0,0),d10=c(0,0,0,0,0));treatments
        design<-model.matrix(~treatments);design
        design<-subset(design,select=-treatmentsd10);design
        colnames(design) <- c("Intercept","Pre","Day3","Day7");design
        fit<-lmFit(goodArrayEset,design)

        • Reply admin 8月 23,2011 7:50 下午

          @mimisikai, 你的数据有几列?是5列吗?每列都一样长吗?如果有可能的话,把你的goodArray的源给我。

          • Reply mimisikai 8月 24,2011 11:36 上午

            @admin,
            > goodData
            AffyBatch object
            size of arrays=1164×1164 features (11 kb)
            cdf=HG-U133_Plus_2 (54675 affyids)
            number of samples=8
            number of genes=54675
            annotation=hgu133plus2
            notes=
            > goodDataAnnotatedTable
            Time.points Subject Group Array
            HJ pre day 0 HJ 1 1
            HJ 3d day 3 HJ 2 2
            HJ 7d day 7 HJ 3 3
            SW pre day 0 SW 1 4
            SW 3d day 3 SW 2 5
            JMX 7d day 7 JMX 3 8
            CYM-042 Low CYM042 4 9
            CYM-144 High CMY144 4 10

          • Reply admin 8月 24,2011 11:48 上午

            @mimisikai, 很明显,你的goodData是8列,而不是你在design中的5列。所以你需要重新设置你的design。

          • Reply mimisikai 8月 24,2011 12:01 下午

            @admin, 但是 我的数据涉及time-course,一个样本有2-3个时间点,该如何设计呢?

          • Reply mimisikai 8月 24,2011 1:14 下午

            @admin, 我把其设计成 col为样本名称 row为时间点,这样的话 就是5X3的matrix。弄成8列如何分析呢?

          • Reply admin 8月 24,2011 2:01 下午

            @mimisikai, 参见http://www.bioconductor.org/packages/2.8/bioc/vignettes/limma/inst/doc/usersguide.pdf第50页。

      • Reply xqyang 10月 23,2011 7:43 下午

        @admin, 我是一个R语言分析芯片的新人,就是简单的用用SAM等程序,希望您以后有更加的新贴子!谢谢。

  7. Reply mimisikai 8月 23,2011 6:35 下午

    您好!我是个bioconductor的新手,有些问题想请教您。
    我的pipeline参考了该连接http://arrayanalysis.org/main.html
    之后我剔除了很糟糕的数据,进行后续分析
    到这儿,我就有些没有思路了,后门应该按照什么样的流程分析呢?
    BTW 我的数据是Affy Human U133 puls2 RNA chip

    • Reply admin 8月 23,2011 7:52 下午

      @mimisikai, 嗯,对你的问题不是很清楚。QA分析完了之后,就应该normalization&summarization,而后就是isoform了啊。

      • Reply mimisikai 8月 23,2011 9:24 下午

        @admin, 那通常 bg correction、normalization、summarization选用什么算法呢?

        • Reply admin 8月 24,2011 8:38 上午

          @mimisikai, 一般来说我们都使用RMA, MAS5/Plier。它们直接包含了bg correction、normalization、summarization等步骤。GCRMA是RMA的一种,它比RMA多了一步分析GC比并权重的过程。有的说挺有用,但也有的说没太大效果。对于RMA,它使用的是RMA bg correction (Irizarry et al 2003), quantile normalization和median来summarization. 它的效果和MAS5比较:

          MAS5 RMA
          log PM No Yes
          Weighs probes in different probe sets differently Yes Yes
          Robust at low levels NA Yes
          Uses iterative method to estimate probe effect and expression index Yes Yes

  8. Reply mimisikai 8月 24,2011 12:42 下午

    BTW
    readTarget fuction如何运用?读取的txt文件是本身就存在的?还是自己建立的呢?

    • Reply admin 8月 24,2011 1:59 下午

      @mimisikai, 那个txt文件要自己写。当然,你也可以通过
      files<-dir(pattern="*.cel",ignore.case=TRUE);files 这样的命令来获取文件名来替换readTargets。之后的ReadAffy(filenames=targets$filename)就可以使用ReadAffy(filenames=files)。 但是使用自己设计好的target文件有个好处,那就是在后面设计design的时候可以很方便的使用列名来设计,比如说你的txt文件有以下几列:filename, factor, pared。你就可以写成: f<-factor(targets$factor); p<-factor(targets$pared); design<-model.matrix(~-1+f+p); colnames(design)<-gsub("^[f|p]","",colnames(design)); colnames(design)<-make.names(colnames(design)); 来实现。

  9. Reply mimisikai 9月 1,2011 12:58 下午

    用maQualityPlots(mraw)
    总是报错Error: duplicate switch defaults: ‘list(dev =…’ and ”

    Why?

    • Reply admin 9月 1,2011 9:32 下午

      无法重复你的错误啊。我没用过这个函数。如果你想让我试试的话,把你的历史记录传上来,并且把sessionInfo给我。

  10. Reply mimisikai 9月 30,2011 5:34 下午

    hi,您好!
    想请教您几个问题,我分开写。
    1. 我看Affy芯片中(Hum U133plus2)有60个 spike control,但通过qcProbe() 和spikeInProbes()只能调出有限的几个参数,如果我想获得所有spike control的数据,怎么调用?例如我要Phe Dap Thr Lys等做几个unlabeled spikes图。
    2. 如果我有5个group,每个group中有5个sample,我想设定一个color参数,即每个group是一种颜色,group内是该色的类似色。如,group1是红色,则group1中的5个sample依次是粉红、淡红。。。之类的。
    想了好久了,也没写明白,谢谢您!

    • Reply admin 9月 30,2011 9:55 下午

      第一个问题:你先去http://www.affymetrix.com/support/downloads/library_files/hgu133plus2_libraryfile.zip下载注释文件。其中两个文件叫做CD_HG-U133_Plus_2/Full/Controls.ini
      CD_HG-U133_Plus_2/Full/HG-U133_Plus_2/LibFiles/HG-U133_Plus_2.psi,你可以通过这两个文件找到所有你需要的probe的索引号,然后再使用这个索引号就可以拿到它的pm,mm值。
      第二个问题:我不是特别明白你的难点在哪里。如果你想使用颜色来分组的话,可以建立一个三列的数组,第一列是样品名,第二列是颜色值,第三列是数据,你在做图的时候保持颜色值和数据的对应关系就可以了。如果你希望有一个亚分组的概念,那基本上也就是这个思路。

      • Reply mimisikai 10月 1,2011 12:55 上午

        @admin, 针对你说的那两个文件 怎么打开呢?

        • Reply mimisikai 10月 1,2011 1:21 上午

          @mimisikai, 对不起 很抱歉问了个傻问题 文件已经打开了。
          BTW 用什么样的命令能够调用出control的pm 或者mm呢
          pm[x,y]吗? 有时候不知道一个function 就解决不了问题。
          谢谢!

          • Reply admin 10月 1,2011 5:48 上午

            @mimisikai, 你可以阅读一下http://bioconductor.org/packages/2.8/bioc/vignettes/affy/inst/doc/affy.pdf当中关于PM和MM的部分。其实对于function,你如果大约知道名字的话,可以使用?和??来获得帮助,也可以使用google。

          • Reply mimisikai 10月 2,2011 10:50 上午

            @admin, 根据您说的 hum U133 plus2.psi文件中的 左边是ID? location? 右边是intensity?
            我用pm【xx,】调用,出来的matrix是一些序列号,但又不像是ID,因为有的都6-7位数字了,然后是样本intensity。
            若我用probe调用,出来的有probe名称 但是不知道具体的每个probe的强度

          • Reply admin 10月 3,2011 9:51 上午

            psi文件中,第一列是索引号,第二列是探针名称,第三列是探针的个数。你可以先用read.delim读入psi文件,psi<-read.delim("CD_HG-U133_Plus_2/Full/HG-U133_Plus_2/LibFiles/HG-U133_Plus_2.psi",header=F)然后把其中第二列当中spike-in的部分读入control,比如controls<-as.character(psi[2:63,2]),然后再使用pmindex把它们的index值都读出来,比如pi<-pmindex(data,controls),最后你可以使用intensity来读取它们的值,比如intensity(data)[pi[[1]], ]。如果你想读取mm的值,可以使用:mmindex。还有一种相同的操作,比如:indexProbes(data, which="pm", genenames="AFFX-BioB-5_at")等等。你都可以自如地读取原始数据。

      • Reply mimisikai 10月 1,2011 1:26 上午

        @admin, 关于颜色问题,我是想在所有的function前设立col的参数。
        并且我想这个col比较通用,就是说我今后把它整理成package或者一个大的function,无论我有多少样本,只要有几个组,就有几种大颜色对应于组,然后在组内,每个样本又使用该系列的颜色,如group1是红色系,则该组内的所有sample都是红色,但是彼此间有差别,但保持是一个色系。

        • Reply admin 10月 1,2011 6:13 上午

          @mimisikai, 我明白你说的了。这样的一个包,你可以这样想,第一,面向对象编程,第二你要定义好一些大色系的亚色系。第三,需要一个setCol及getCol之类的函数。第四重载或者复写使用到col的函数。比如你定义好颜色,下面的亚色系有十几种。你有多个组,你setCol给它们颜色的过程就是枚取每个组内成员赋取亚颜色。对于颜色定义,你可以使用list,比如col<-list("red"=colors()[grep("red",colors())],"blue"=colors()[grep("blue",colors())],...),调用时就可以使用col$red来调用红色内的亚色系。有组group1<-c("s1","s2"); group2<-c("s3","s4"); group<-list("g1"=group1,"g2"=group2)。你建立一个对象,它有两个成员,一个是group,一个是gpcol,它有四个方法,一个是setGroup, 一个是getGroup, 一个是setCol,一个是getCol。然后就使用setGroup来初始化对象中的group,而后使用setCol对它依次设置颜色。以后在调用其它关于col的函数时,分组数据由getGroup来完成,相应的,颜色由getCol来完成。

  11. Reply mimisikai 9月 30,2011 5:40 下午

    还有个问题,对数据分析时:
    background adjustment中有RMA、MAS5.0、等方法;
    Normalization中有constant、invariantset、quantiles、loess、contrast等方法;
    使用时咋么选择呢?还是都试一遍,然后比较呢?
    之后我又看见summarization中有expresso、RMA、GCRMA等等,这个是不是将前面的bg correction和norm都综合在一起了?就是说,比如只要RMA这一个方程搞定所有的了?
    最后一个问题, 当做完norm后,准备设计design matrix前, filtering probe set是不是必要的步骤?一遍筛选严格程度是多少呢?
    非常感谢!

    • Reply admin 9月 30,2011 10:08 下午

      对,你的理解是正确的,使用RMA或者GCRMA之后,bg correction 和normalization都会自动完成。你仔细阅读相关的文献就会知道,其实RMA就是使用的rma bg correction, quant-norm, & med-polish。而GCRMA就是在此基础上参考了GC比。而之后的filtering probeset的步骤要看你的芯片。如果是3’IVT expression芯片的话,不是必须的。但如果是Whole-Transcript Expression,对于外显子水平是必须的,而对于基因水平,不是必须的。一般来讲,基因水平的表达差异比较都是在之后使用foldchange和p.value, adj.p.val来筛选的。而对于probe水平的筛选,通常来讲,RMA已经做出了筛选,在算法上,RMA首先是做了聚类分析,而后才选择了一部分probe来总结出表达强度。所以你不用过于担心之前的筛选步骤。

  12. Reply mimisikai 10月 4,2011 3:16 下午

    又来请教您了。
    我的数据 来自human U133 plus2。0
    我发现我的数据 在做nsFilter后 原先50000多个原始值,现在只被过滤的剩下25000多个了。我记得上次您说过,针对3IVT的芯片,做不做nsfilter都行。但我这个咋的差距那么大呢?

    • Reply admin 10月 5,2011 8:20 上午

      @mimisikai, nsFilter不单单会过滤掉信号强度弱的,还会去掉不同样品间变化很小的,所以,只剩下一半是很正常的。你不可能期待着你的样品超过一多半都是摆动很大的吧?

  13. Reply mimisikai 10月 6,2011 3:16 下午

    hi 又有两个问题想请教您
    topTable 中,sort.by 及method有很多个选择,通常应该作何选择呢?
    相似的问题,decideTests中method四个方法 separate、global、hierarchical、nestedF,我挨个试了,发现结果有较大差别,应该选择哪个方法呢?

    • Reply admin 10月 6,2011 4:21 下午

      @mimisikai, topTable是非常常用的一种函数,简单地来讲,就是显示表格排序最靠前的部分,因为有排序,所以就会要求选择sort.by。选择sort.by就是在问你对什么最感举兴趣,你如果对foldChange最感兴趣,你就选logFC,如果你认为p.value最重要就选p,如果你认为t检验的值最重要,就选t。或者说你觉着limma本身给的评分最重要,就选B。按我的经验,这个时候一般都使用logFC或者B值。而且一般地,我都会将number设置成nrow(fit),其中fit就是你传入的线性拟合的结果,这样呢,就可以把所有的值全部都拿到,以方便以后改变主意,想换换其它的排序方法。至于decideTests中的method方法,简单地来讲,separate就是等同于使用topTable,并且以B值排序,global就是使用t检验,hierarchical就是使用聚类分析,而nestedF和t检验很类似,它使用的是F检验。至于F检验和t检验的差别就在于F检验验的首先是样本之间的方差是否相同,如果相同就接着使用t检验,如果不同,就是使用t’检验。这个nested之外的意思就是会比直接使用t-test多花些准备的工夫,但很有可能在结果的准确性上更好。

      • Reply mimisikai 10月 7,2011 1:25 下午

        @admin, 太感谢您了!感谢您的耐心讲解!
        不知道您是否处理过Agilent的数据,我发现它的数据结构上,在结果中不显示有些未标识的基因,导致实际结果的probe数和设计probe数量不一致,使得后续用limma或者marray处理时,总是报错。

        • Reply admin 10月 8,2011 7:25 上午

          处理过一些,但从未遇到你说的这个情况。我觉得很可能是你需要自己设定annotation和columns。具体的,需要你把进一步的信息给我,比如说文件头。

  14. Reply mimisikai 10月 9,2011 3:57 下午

    hi 您好!
    不知道您是否可以帮我review一下 我写的关于Affy U133plus2 分析的pipeline?

  15. Reply 糗世界 10月 9,2013 10:09 上午

    谢谢您与我们分享Bioconductor的使用经验,我获益匪浅,特请教如下问题,望博主回答。
    1.本篇讲Bioconductor使用R的版本是?在安装某些软件包时提示某些软件包建造的版本和使用R的版本不一致,如何解决?
    2.可否告诉我下分析外显子程序的密码?谢谢。

    • Reply admin 10月 9,2013 10:32 上午

      版本可能会带来一些影响。我当初写的时候没有留下sessionInfo是一下遗憾。

      没有带有密码的分析外显子的博客。唯一一篇带密码的博文是分析普通芯片的,密码是0123456987

  16. Reply taotao 8月 28,2014 4:26 上午

    博主您好,构建实验设计矩阵这里不是很懂。比如说我有8个样本要构建前四个与后四个样本的对比该如何设计?可否告知您的其他联系方式,比如说QQ之类的,谢谢博主。

  17. Reply Ares 9月 28,2014 10:38 下午

    欧老师您好, 您的博客在讲解CyberT算法时(http://pgfe.umassmed.edu/ou/archives/2083/4)最后为什么不把“怎么根据P.Value对结果报表进行筛选,得到全部差异表达基因”的内容也写上去呢?

  18. Reply Ares 10月 10,2014 10:13 下午

    欧老师,我实在是没在到你说的博文啊,你发个链接可以么。万分感谢

  19. Reply lianghui 11月 1,2014 3:46 上午

    老师您好,我从网上下载了Affymetrix Human Genome U133 Plus 2.0 Array平台的一系列数据,想将其中的差异lncRNA找出来,并进行后续的差异表达分析及GO等分析,具体我该如何操作,有没有详细的教程,非常感谢

  20. Reply yjw 12月 4,2014 6:55 上午

    老师您好,多实验因素的差异分析该怎么做,是应该两两比对吗,我下载了一个78个样本的数据,里面既有时间点,又有多种处理,我该怎么分析?

    • Reply admin 12月 4,2014 9:41 上午

      我的建议是先做简单的,不要考虑那么多,先只考虑两两比对,等你的基础知识掌握了,再学习design的设计,学习的时候,以limma的说明文档为圣经。

  21. Reply yjw 12月 4,2014 7:55 上午

    具体是是老年和年轻,两组,又有损伤,未损伤,和未处理3个水平,训练持续5天,21天两个时间点

  22. Reply yjw 12月 4,2014 9:22 上午

    老师您好,我又仔细想了一下,具体是动物有老年损伤,老年未损伤,老年未分类,以及年轻,共4类, 因素是训练21天,训练5天,以及未训练3种,但是,老年未损伤和老年损伤,又是根据年轻组的训练结果进行区分的,是不是4X3类型的?抱歉耽误您时间,有点复杂,实在想不出应该怎样构建表达矩阵。。

  23. Reply 杨峰 12月 19,2014 3:16 上午

    老师您好!我手头上有个数据跟您的例题很像,毒物处理 24h 72h看应答,然后我照您的流程处理有三个问题想请教您!谢谢!
    1、我的每个组有三个重复,是不是要把 treatments <- factor(c(rep(1:4,each=2)),labels=c("e10","E10","e48","E48"))中的each=2改成each=3就行了
    2、我想把差异表达中p.value1这个条件怎么加上
    3、找到那些共同上调和下调的基因怎么把他们输出出来
    真心感谢您!自己捣鼓了N久后来看了您的这个网站才看到点希望!

    • Reply admin 12月 19,2014 4:06 下午

      1.是的。
      2. 不是很明白你的意思。如果你是想使用p.value + logFC两条来筛选的话,你只需要使用
      res <- topTable(fit2,coef="constract.name.you.want",n=nrow(fit2))得到全部, 然后使用你的条件去筛就好了。 3. 得到up or down regulated基因之后,你可以使用merge,或者intersect来取得相同的部分。你还可以参与更多我的博客。

  24. Reply 杨峰 12月 20,2014 1:55 上午

    谢谢您的耐心解答,我正计划把您的博客全都看一遍,真心感谢您!

  25. Reply 杨峰 12月 21,2014 7:06 上午

    老师您好!不好意思又打扰您!有三个问题想请教您!
    1、decideTests和classifyTests这两个函数有什么区别吗?我想实A:24h对照和处理比较得到差异表达基因,B:72h对照处理得到的差异表达基因
    A与B中相同的基因有哪些不同的有哪些
    2、我按照这个dif<-topTable(fit2,coef="E10",n=nrow(fit2))
    dif<-dif[dif[,"P.Value"]<0.05,]找到的差异基因数目跟您上面用 summary(results)得到的上调和下调基因总数不相同
    3、eset<-rma(ab)与eset.e<-expr(eset)这两个的不同是后面那个变成eset的矩阵了吗?
    谢谢您!祝您工作顺利!万事如意!

    • Reply admin 12月 23,2014 10:07 上午

      问题1,在classifyTests里有这样一段:Note that these functions do not adjust for multiple testing across genes. The adjustment for multiple testing is across the contrasts rather than the more usual control across genes. The functions described here are called by decideTests. Most users should use decideTests rather than using these functions directly.意思正好可以回答你的问题,就是decideTests调用了classifyTestsT/F/P等,建议你最好使用decideTests,而不使用classifyTests.
      问题2, 是会不同的,因为算法根本就不一样。你试试results.sep <- decideTests(fit2, method="separate"); summary(results.sep) 问题3, 对。expr就是只把eset中的一部分提取出来了。你具体的可以看一下什么是ExpressionSet。或者查看一下视频:http://pgfe.umassmed.edu/ou/archives/3413

  26. Reply 尔霖 7月 2,2017 10:38 下午

    博主你好,生物信息新手,但偏偏我要用的芯片都是Agilent平台的,一头雾水不知从何下手啊,试了下limma包和marray包,但没找到能看懂的教程。能否指点一下,一步怎么办,或者有R脚本就更好了,谢谢您先

Leave a Reply

  

  

  

%d 博主赞过: