bioconductor系列教程之一分析基因芯片上

至于这些方法之间具体的区别,我不做详细的描述了。我们简单地了解一下lowess(aka loess)和quantile。双色芯片经常使用Lowess(locally weighted linear regression)标准化方法。我们通常认为芯片印刷是均匀的,反应动力学是一致的,显色及色彩捕捉也是均匀的。但实验过程中往往不能达到如此理想的效果,所以即使芯片本身尺寸很小,还是会有不均一的背景以及反应动力学。Lowess把芯片分区块处理,以期待达到更好的标准化效果。Quantile 标准化方法很类似最简单的常数标准化,它也是使用最为广泛的标准化方法。常数标准化只是通过寻找中点的方式简单地定义中点为零,把数据分布的中点的两侧。而quantile标准化方法将数据排序并以指定单位间隔划分成不同区块,并将同一区块内的数据平均后用平均值填充所有该区块内的数值。更多的请参考http://www.bioconductor.org/packages/2.7/bioc/vignettes/affy/inst/doc/builtinMethods.pdf

我们来比较一下基因芯片数据处理当中不同分析方法:

方法 背景较正方法 标准化方法 综合表达定量方法
MAS 5.0 mas constant Mas
dChip mas invariantset liwong
RMA Rma quantile medianpolish
ZAM Rma Loess A robust model is fit (can download R package samr)
GL None Loess Rma

我们已经了解了在R中的两个综合应用MAS和RMA,在R中我们可以简单地使用mas5(Data)或者rma(Data)来直接获得表达定量的结果。它们的区别再次简单总结如下:

Figure 3 基因芯片两种基本预处理方法(MAS5,RMA)流程图。TBW: one-step Tukey BiWeight.
Figure 3 基因芯片两种基本预处理方法(MAS5,RMA)流程图。TBW: one-step Tukey BiWeight.

下面需要探讨的问题是质量控制。这是对原始数据处理的更高一级的要求。

我们首先要对基因芯片的实验设计再多讲两句以便能更好地理解后面的介绍。一般的,基因芯片实验都会同时做两组,实验组及对照组。然后用两组实验结果进行比较来寻找在不同实验条件下差异表达的基因。为了方便,我们用R和G来分别表示它们(通常红色(Cy5)和绿色(Cy3)是较为常用的荧光指示)。

虽然通常情况下,我们认为,荧光的强度大小与mRNA或者cDNA人丰度是成正比的,但是我们知道,尤其是在两极,这个比例关系并不一定存在,比如背景的自发荧光会干扰低表达量的基因检测;而高荧光强度时,荧光检测探头可能会被饱合,而无法区分检测上限以上强度的荧光强度。最简单地用于显示原始数据这种荧光强度依赖效应的方法是绘制R-I(ratio-intensity)图,又称MA图。MA图中定义M=log2(R/G); A=1/2*log2(R*G)。我们来看示例:

> library(affydata) ##装载库文件

> data(Dilution) ##载入数据(库文件中附带的示例数据)

>phenoData(Dilution) ##查看数据内容与结构

An object of class "AnnotatedDataFrame"

sampleNames: 20A, 20B, 10A, 10B

varLabels and varMetadata description:

liver: amount of liver RNA hybridized to array in micrograms

sn19: amount of central nervous system RNA hybridized to array in micrograms

scanner: ID number of scanner used

> pData(Dilution) ##不同组实验数据描述

liver sn19 scanner

20A    20    0       1

20B    20    0       2

10A    10    0       1

10B    10    0       2

##使用MAplot为Dilution绘制MA图。
##将生成一个全部数据的一览图。

> MAplot(Dilution, pairs = TRUE, plot.method = "smoothScatter")

这里的参数pairs设置为TRUE的意思是两两都要比较的意思,默认情况下是只比较实验组与对照组,plot.method设置成smoothScatter会加速绘图过程,在默情况下是对数据中的每一个点都运算绘图,设置成smoothScatter后会抽取部分点运算绘图。结果如Figure 4所示。

Figure 4 MAplot
Figure 4 MAplot

51 thoughts on “bioconductor系列教程之一分析基因芯片上

  1. 求助,我一直不知道cdf文件怎么下载,下载完之后放到哪里,我已使用readaffy,就会有Error in getCdfInfo(object) :
    Could not obtain CDF environment, problems encountered:
    Specified environment does not contain Mapping250K_Nsp
    Library – package mapping250knspcdf not installed
    Bioconductor – mapping250knspcdf not available
    这个错误提示 ,但是mapping250knsp.cdf这个我到处下载到处找都不知道怎么使用他,都快急疯了

  2. 你好!我在运行>eset <- expresso(Data, normalize.method="invariantset",
    bg.correct="rma",
    pmcorrect.method="pmonly",summary.method="liwong")的时候会报错:
    Error in if (bg.correct) { : argument is not interpretable as logical
    请问这是什么原因吗?

    1. bg.correct参数只接受logical value。我想你可能是想传递值”rma”给bgcorrect.method参数。注意那个点(.)。

    1. 算法都是一样的,唯一的不同就是如何芯片的设计不同。比如agilent芯片,可以使用limma.当然也可以使用agilp.

  3. 欧博士,您好!我在按照您的教程学习芯片分析,在学习过程中除了重复教程实例外我还参考了其他资料和加入自己的笔记,但是总体框架还是您的博文的内容。我现在想把学习的笔记放入自己的博客,笔记内不可避免地提到了您的博文内容,不知您是否能授权我在自己的博客上加入这部分内容?我的博客非商业性质。

    1. 非常欢迎。只要注明出处即可。我当然希望分享我的学习笔记,希望更多的人看到,并且讨论。也希望广交朋友,互相交流促进。

  4. 博主你好,我又来问问题了
    Data <- ReadAffy()意思不是读取工作目录下的CEL文件
    吗,然后赋给Data吗?然后再执行image(Data[,1]),的确能画出第一张芯片的灰度图像,然后我想查看数据类型,执行data(Data)就出错了,显示:data set ‘Data’ not found,这是什么回事呢?假如要改,怎么改呢
    谢谢博主

    1. 你要查看一个对象的类型,可以使用class或者mode命令。data命令是调用库中的数据。你可以通过?data来了解更多。

  5. 你好,
    我购买了《R语言与Bioconductor》这本书,,看到“5.3 基因芯片预处理”这一节的内容,里面介绍了AffyBatch数据结构,书中说,这个基因表达矩阵列对应不同的样品,行对应不同的探针组(Probe sets),用一个无重复的索引值表示。但是我认为,行对应的是探针(Probes)而不是探针组(Probe sets),因为一个索引对应一个坐标(索引=Y*列数+X+1),一个坐标对应一个探针

  6. 其实我在那个矩阵里还发现了一个规律:就是每一列的数据都是“小大小大小大……”,例如CLL10.CEL对应183,11524,301,11317,115,别的列都是这样。我再看到灰度图像,灰度图像其实就是一个探针对应一个格子,而第一行是芯片的边框,呈现“黑白黑白黑白”这个规律,假如这个矩阵是1164*1164的,那么前1164个格子应该都是“黑白黑白黑白”,对应的强度是“小大小大小大……”,所以我会认为每列数据呈现的是探针的强度,而前1164个数据都是“小大小大小大”

  7. 正是因为看到了图5-6,所以以我才认为那个矩阵的行对应探针而不是探针组啊……
    后来我又找来了本地的2个CEL文件做了实验:
    rawData<-ReadAffy("GSM286756.CEL","GSM286757.CEL")
    rma<-rma(rawData)
    e<-exprs(rawData)//其实第2,第3行拿来干嘛呢,只是不明白为什么不打上去的话以下的代码会出错
    exprs(rawData)[1:10,1:2]

    显示:
    GSM286756.CEL GSM286757.CEL
    1 83 93
    2 6592 7675
    3 120 116
    4 6645 7955
    5 132 125
    6 63 108
    7 6507 7869
    8 87 99
    9 6595 7780
    10 92 96

    然后我打开GSM286756.CEL,前10个探针的强度是这样显示的:
    CellHeader=X Y MEAN STDV NPIXELS
    0 0 83.0 15.8 25
    1 0 6592.0 1136.3 25
    2 0 120.0 25.2 25
    3 0 6645.0 932.4 25
    4 0 132.0 21.3 25
    5 0 63.0 9.5 25
    6 0 6507.0 812.6 25
    7 0 87.0 17.6 25
    8 0 6595.0 1126.1 25
    9 0 92.0 15.6 25

    所以我认为每一行应该对应一个坐标(X,Y),指定某个探针而不是探针组

    1. 具体你可以查看:http://www.affymetrix.com/support/help/faqs/mouse_430/faq_8.jsp。
      至于为什么你会产生这样的误会,就是因为CEL文件其实已经是重新编码过的芯片表达值的记录,而不是原始扫描图像。

  8. 其实CEL文件是由DAT文件得到的(过程大概就是:DAT文件中是用7*7个像素表示一个探针,然后去掉最外那圈,变成5*5个像素,然后就用这25个像素来求一个平均值,作为该探针的原始强度),在这个过程中,并没有把探针汇总到探针组水平,所以CEL文件里MEAN的那列就是一个个探针的原始强度,而那个矩阵显示的数据和CEL文件里MEAN的那列一模一样,所以也应该是探针水平而不是探针组水平。假如要汇总到探针组水平,则需要CDF文件。
    我们这段时间都在研究与芯片有关的各种文件格式,写了很多读写这些文件的代码(所以我才敢做修改原始文件那种事),也作了一些图,整体来说还是有一定成果的。
    就想请问一下是哪位老师负责编写你们这本书中有关质量控制和预处理算法介绍部分的,我们能和他取得联系吗?我们希望可以交流一下,假如能够给我们的工作给予支持,那真的是感激不尽了!!

    1. 不会意思,之前的错误。我看了一下,你说的是正确的。affybatch中保存的是相同cdf下的probe level的数据,不是probeset水平的数据。这些数据经过预处理之后才会得到校正,在综合后得到probeset水平的数据。由此看来知识放一段时间就忘。

  9. 已经给您回邮件信息了,也不知道您能不能收到qq邮箱发的。arrayanalysis能够作出很多有关质控和预处理的图表(如那两类嵌入式探针的强度趋势图等等,书上就找不到了),其实挺希望您的书改版后增加一些内容。

  10. 你好,有没有关于双通道芯片数据读入到质量控制,再到normalization的步骤啊?谢谢:)

  11. 博主
    您好
    我初次接触芯片数据,我现在又大约600个基因,是一个家族的,我想通过分析芯片数据得到这些基因在植物不同组织中的表达量数据,然后用R画一个热图。我现在遇到的问题有三个:
    第一,我看了一些SCi文献中有提及这一方法,但是流程模糊,对于初学者来说,完全不知所云,因此我想在您这里求助一个大致的分析流程。
    第二,我看到网上关于基因表达的有一些教程,但都是关于基因差异性表达分析的,不知道这个和我所要分析的内容是不是有些相通点呢。
    第三,我们如何识别所得到的表达数据那些是我所关注基因的呢?
    感谢您的指导,谢谢。

    1. 因为问题都太广泛了,很难回答您。但是非常感谢您的留言。我想说的就是,任何一个问题,当你自己试着去找解决方案的时候,往往比你拿别人的答案来抄收获更多。我相信,对于affymetrix公司的芯片,本博客大多都介绍过。对于其它的芯片,可以通过扩展阅读来解决。本博客的目的是为了方便阅读英文材料有一定困难的初学者,以及有相关问题希望讨论的研究人员。所有,网络上会有很多相近的材料,都是可以参照的。本博客也有不少的错误,我努力在改正,但个人水平有限,无法全部改好。对于第三个问题,你可以对芯片表达数据做完注释后,通过找并集的方式找到你感兴趣的基因。

  12. 照猫画虎地执行你的命令:
    eset <- expresso(Data, normalize.method="invariantset",
    bg.correct="rma",
    pmcorrect.method="pmonly",summary.method="liwong")

    出来的错误信息是:
    Error in if (bg.correct) { : argument is not interpretable as logical

    看见你的回答如下。但还是不明白。什么是“那个点(.)”?
    bg.correct参数只接受logical value。我想你可能是想传递值”rma”给bgcorrect.method参数。注意那个点(.)。

  13. 您好!
    执行library(affydata)报错, Error in library(affydata) : there is no package called ‘affydata’;请问问题在哪呢?
    非常感谢

  14. 您好!在《R语言与Bioconductor》一书中,芯片质量控制一节中,提到在箱线图、RNA降解曲线、取类分析等作图,都需要affyPLM包实现,可在您的这部分教程中,并没有提到。那么affyPLM核心作用到底是什么?
    非常感谢!

    1. affyPLM的核心是对probe-level进行基于模型的拟合。书中的部分主要来源于http://pgfe.umassmed.edu/ou/archives/2042/6

  15. 您好!
    我在绘制箱式图时boxplot(Pset, col = cols),图中细胞名称如MDB-MB-231,会自动变成MDA.MB.231。这该怎么调整?
    非常感谢!

  16. 执行MAplot 画图出现错误 MAplot(Data, pairs = TRUE, plot.method = “smoothScatter”)
    Error in simpleLoess(y, x, w, span, degree, parametric, drop.square, normalize, :
    NA/NaN/Inf in foreign function call (arg 1),该怎么调整,上次用这个的作图时候没有提示错误?

  17. 非常感谢您的耐心解答!
    当我运行,Gset <- gcrma(data)时,报错:Error: length(prlen) == 1 is not TRUE。(data: formal class AffyBatch)。
    Q:1)是芯片差异?,我分析的是miRNA-1_0; 2)gcrma的意义是converts an AffyBatch into an ExpressionSet , 我还不太明白,在进行聚类及PCA 分析之前,这一步的意义是什么?直接按您的教程《使用Bioconductor下载GEO 上的数据》第一种方法下载的表达数据是ExpressionSet吗?3)miRNA芯片背景校正、标准化和汇总及之后操作与教程上相仿吗?
    非常感谢!!

    1. 你的问题很多。关于问题2,你可以使用class()函数来查看下载下来的是什么对象。关于gcrma的意义,你需要从什么是rma开始。这需要阅读一下文献:Summaries of Affymetrix GeneChip probe level data. Irizarry RA, Bolstad BM, Collin F, Cope LM, Hobbs B, Speed TP. Nucleic Acids Res 2003: 31(4)。在R当中rma的操作分为了三步:一,rma-bg:Performs an RMA style background adjustment as described in Irizarry et al 2003.就是我让你读的这篇,二,quant-norm:Class for doing quantile normalization.三,med-polish:Performs a median polish to estimate target and probe effects. Resulting summaries are in log2 space by default. Used in summary step of RMA as described in Irizarry et al 2003. 我想你理解了这些,就不难理解gcrma了。
      关于问题3,如果你下载的CEL文件的话,是相仿的。

  18. 前辈,您好!我有问题请教
    所有的芯片(cDNA,miRNA,methylation等)质量控制 标准化等预处理后得到的数据, 在挑选差异表达基因时采用的方法(比如 limma包) 是不是都一样?

  19. 老师,您好。我用illuminaio包读取illumina芯片数据的.idat.文件后,得到了ProbeID,AVG_Signal和NBEAD等数据,怎么利用这些数据计算出表征信号强度和背景差别的显著性的信号检测P值呢?

  20. 您好,想请假怎么把Affymetrix AFFY_6.0 碱基用A、B的文件转换为用AGCT表示?
    比如:一个snp对应的plink-bim文件是这样的,0 SNP_A-1781650 0 0 A B,但是想转换为类似这样的,0 SNP_A-1781650 0 0 G T

  21. 请问老师可以写个oligo分析Affymetrix Human Exon 1.0 ST Array的教程么?我学习了您的博客,但是还是发现很多错误不知道如何解决,实在痛苦

  22. 请问大神,这是我的代码:library(affy)
    data <- ReadAffy() 【文件在当前目录下】
    eset <- rma(data)
    MAplot(eset, pairs=T,which=c(1,3,5), plot.method = "smoothScatter",cex=1)
    然后就报错了Error in (function (classes, fdef, mtable) :
    unable to find an inherited method for function ‘MAplot’ for signature ‘"ExpressionSet"’这是为什么?但是我对data做MAplot图就没有问题
    之后我又写代码write.exprs(eset,file="data.txt")
    library(oligo)
    celFiles <- list.celfiles(getwd())
    affyRaw <- read.celfiles(celFiles)
    library(pd.mogene.2.0.st)
    eset <- rma(affyRaw)
    MAplot(eset, pairs=T,which=c(1,3,5), plot.method = "smoothScatter")
    这次就成功了。
    另外,请问Can you figure out the commands to use so RMA outputs expression index for RefSeq (e.g. NM_000546) vs Gene Symbol (e.g. TP53)这句话指的是什么?

    1. 在做MAplot之前你试着library(affyPLM)一下。oligo本身自带了一个MAplot.
      你的第二个问题我也看不懂。这是作业吗?应该是在考你注释问题吧。你可以试试biomaRt。

  23. 您好,请教个问题啊。现在想看三篇文章的基因芯片中A基因的表达量,但是用的基因芯片都是不一样的,所以说他们有没有什么办法能经过标准化后可以互相比较。或者可不可以找些看家基因来作为对照,看A基因对看家基因的相对表达量,以达到不同芯片之间的比较。谢谢。

  24. 老师您好,我在处理一个[HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array 类型的数据文件,将样本病例组和对照组分组放在两个文件夹中,设定其中一个文件夹为工作目录,用以下语句读取文件夹中的CEL文件
    celFiles <- list.celfiles('myCELs', full.names=TRUE)
    rawData <- read.celfiles(celFiles)
    出现如下错误
    All the CEL files must be of the same type.
    Error in read.celfiles(celFiles) :
    checkChipTypes(filenames, verbose, "affymetrix", TRUE) is not TRUE
    求问这是什原因,该如何处理呢?万分感谢您的回答。

    1. 试一下
      pkgname <- read.celfiles(celFiles[1])@annotation ## could be any one you can read. 或者你可以直接试一下:pkgname <- "pd.hugene.1.0.st.v1" rawData <- read.celfiles(celFiles, pkgname=pkgname)

发表评论

电子邮件地址不会被公开。 必填项已用*标注