bioconductor系列教程之一分析基因芯片中(质量控制)

使用FitPLM生成权重,残差及NUSE图像

使用simpleAffy对基因芯片质量所做出的评估都是基于平均值的,它有一个默认的假设,那就是对于每一块芯片,质量是均匀的,不会因在同一块芯片内部还随着位置的不同而质量发生较大的变化。但事实上往往不是这样。如果我们把芯片划分成许多小格的话,我们会发现格与格之间的质量也是有差异的,这可能是由于芯片印刷的问题,也可能是杂交过程出现的问题,导致了同一块芯片也出现质量不均匀的情况。这种对待芯片的态度,就好象MAS5和RMA之间的差别。

为此,我们导入affyPLM库。这个库会基于RMA算法来对数据进行预处理,然后给出与质量控制相关的直观图像。

> library(affyPLM)
Loading required package: preprocessCore

Attaching package: 'affyPLM'

The following object(s) are masked from 'package:stats':

    resid, residuals, weights

> Pset <- fitPLM(Data) 
> x11()
> par(mfrow=c(2,2))
> image(Data[,2]);image(Pset,type="weights",which=2,main="Weights");
> image(Pset, type="resids", which=2, main="Residuals");
> image(Pset, type="sign.resids", which=2, main="Residuals sign");

什么是所谓的权重(weight)呢?什么是残差呢(residuals)?如何看懂权重残差图呢?首先我们要了解的是affyPLM库在探针水平(probe-level-model)拟合时所使用的回归方法是最小二乘法,而普通最小二乘法假设误差项的方差是不变的。然而,在基因芯片的应用计算过程中,这一假设是不成立 的。所以引入了加权最小二乘法来进行回归。而这个权重呢,就体现着这种方差的变化。残差就是指期待值和观测值 之间的差异。因为所有的探针都是随机分布在基因芯片上的,所以我们应该期待权重和残差的分布也应该是随机的。具体到权重残差图片中来,一般的,在权重图中,绿色代表较低的权重(接近0),白色、灰色代表较高的权重(接近1);在残差图中,红色代表正的高残差,白色代表低残差,蓝色代表负残差。而sign.resids图使用红色表达正的残差,蓝色表达负的残差。因为我们期待权重和残差都是随机分布的,所以我们应该看到的是绿色均匀分布的权重图,或者红蓝均匀分布的残差图。示例的图片中,我们可以看到在一定程度上芯片的右侧有部分权重及残差分布并不均匀。我们也会注意到,图中会出现一些白色的条块,这是正常的现象,因为有些时候,探什会按照GC比排布从布导致白斑的出现。

Figure 2权重及残差图
Figure 2权重及残差图

在同一组实验中,即使是相互比较的对照组与实验组之间,大部分基因的表达量还是应该保持一致的,何况平行实验之间。当我们使用相对对数表达(Relative Log Expression(RLE))的的箱线图来控制不同组之间的实验质量时,我们会期待箱线图应该在垂直中央相类的位置(通常非常接近0)。如果有一个芯片的表现和其它的平行组都很不同,那说明它可能出现了质量问题。

为了使图的颜色好区别,我们使用RColorBrewer库来直接调用R为我们准备好的颜色组。

> library(RColorBrewer)
> colors <- brewer.pal(12,"Set3") 
> x11()
> Mbox(Pset, ylim = c(-1,1), col=colors,main="RLE")
Figure 3相对对数表达箱线图
Figure 3相对对数表达箱线图

相对标准差(Normalized Unscaled Standard Errors(NUSE))是一种比RLE更为敏感 的质量检测手段。如果你在RLE图当中对某组芯片的质量表示怀疑,那当你使用NUSE图时,这种怀疑就很容易被确定下来。NUSE的计算其实也很简单,它是某芯片基因标准差相对于全组标准差的比值。我们期待全组芯片都是质量可靠的话,那么,它们的标准差会十分接近,于是它们的NUSE值就会都在1左右。然而,如果有实验芯片质量有问题的话,它就会严重的偏离1,进而影响其它芯片的NUSE值偏向相反的方向。当然,还有一种非常极端的情况,那就是大部分芯片都有质量问题,但是它们的标准差却比较接近,反而会显得没有质量问题的芯片的NUSE值会明显偏离1,所以我们必须结合RLE及NUSE两个图来得出正确的结论。

> boxplot(Pset, ylim= c(0.95, 1.20), col=colors, main="NUSE")
Figure 4相对标准差箱线图
Figure 4相对标准差箱线图

37 thoughts on “bioconductor系列教程之一分析基因芯片中(质量控制)

  1. 上面说的MCF7数据我怎么下载不了啊,博主能不能发给我邮箱啊,昨天通宵学习您的博客,中间卡壳了。。
    谢谢。。

    1. 数据很大,邮箱发送不方便。其实你随便找一组数据就可以了。从我提供的本站链接下载应该可以的。

  2. > Data.qc = qc(Data)
    > Data.qc
    请问这里面的qc是哪个函数?我输入完之后提示没有qc这个函数~
    谢谢解答~

    1. 你很细心。这里需要注意的是,表达值不一定都是实值,很有可能是log2转化后的值。所以,这里的3’/5’因为是对数转换过的,依据公式log2(A/B) = log2(A) – log2(B),应该就是减一下。

  3. 欧博士你好!
    我非常好奇bioB的强度表达值是否只用AFFX-r2-Ec-bioB-3_at这个探针组来计算,于是我直接打开version 3的CEL文件,用型号对应的CDF文件来给探针定位,先做了2组实验,分别把AFFX-r2-Ec-bioB-M_at和AFFX-r2-Ec-bioB-5_at这两个探针组的11个PM探针(各11个)的强度修改为0(原本有几位数就保留几位数,补0就行),再用R来计算bioB的强度,发现计算结果和修改之前有细微的差别,这是否能说明M端和5端与bioB强度的计算有一定的关系?后来再做一组实验,只把AFFX-r2-Ec-bioB-3_at的11个PM探针强度置0,用R算出了-5.469394这个数,2^(-5.469394)只是约等于0,但不等于0。
    欧博士能为我解答上述问题吗?

    1. 对于bioB的强度表达值,很多都会使用3’端的值。当然,也可以综合其它的探针。但是依据RNA会从5′->3’降解来看,使用3’端的值会尽可能的排除降解这一因素的干扰。因为我没有对CEL文件做过你这样的修改,所以只能是猜测一下。你可以查看一下你所有的探针,在读入R之后的最小值是多少。是否可能系统为了避免log2(0)这种错误而对所有的值都加上了一个基数?

  4. 还有就是,博文中提到:
    探针的排序以组为单位,随机排列。而每一组,都由20对探针组成。这一组探针被称为探针组(probeset)。
    但是根据我的了解,类似AFFX-BioB-5_at这种探针组有20对探针,但是AFFX-r2-Ec-bioB-5_at这种有11对探针,而计算的话应该是使用有r2的那组探针组

    1. 现在,对探针数的概念已经变化了。已经不要求每组都有20对探针了。所以我原文里的部分是不准确的。谢谢指正。

  5. 欧博士你好,我又来问问题了。
    上次说到我直接修改CEL文件来做对照实验,在只修改与BioB相关的探针组的强度的情况下,发现BioC,BioD,CreX的强度表达值也受影响了(影响很小),所以我怀疑是不是先对做原始数据一个统一的处理,然后再用于计算。就是说,先要把探针水平的数据(CEL水平)变成表达水平的数据,再用于计算。这似乎是预处理的内容,难道是先进行预处理,再做质量控制?

    1. 这个我没有实验过。你很有钻研精神。我猜想,当你修改了BioB的值之后,会影响背景噪音的估计,从而影响其它的强度表达值。但这一点不能肯定。另外,不需要预处理再做质量控制。

  6. 你好,
    以前我总认为MA图里的R是红色,G是绿色,怎么想都想不通红色和绿色是怎么来的。后来再看看你上述的例子(MVA plot),六个样品,右上部分有15张图,C(6,2)=15,所以后来想是不是两个样品对应一个图

    1. 你的理解本身没有什么问题,因为MA plot最开始就是你所谓的R是红色,G是绿色, 它们本身代表两种荧光(Cy5(红), Cy3(绿))的颜色,而这两种荧光就分别代表了两组样品的杂交强度。而后来就直接跳出了荧光的颜色了。

  7. 欧博士,您好,想问您一个问题,
    基因探针中的序列除了基因的序列,为什么有时候也有基因上游的序列和基因间的区段?可否将答复发送到我的邮箱xujun_njupt@126.com,万分感谢。

  8. 欧博士,您好,
    基因芯片平台GPL8434(http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?view=data&acc=GPL8434&id=6113&db=GeoDb_blob33),其中有一些探针是这样的:UP_SOXXXX,或IG_SOXXXX_SOXXXX,对这些探针的描述分别是Upstream sequence 和Intergenic sequence。

    我在做基因差异表达分析的时候得到的差异表达基因有时候出现了这些基因,虽然都被我忽略了,但是我很好奇为什么在设计基因芯片平台GPL8434的时候要用到这些序列呢,而其它的一些平台中并没有这些序列。

    祝好,谢谢

    1. 这应该是一个自定义的芯片。我猜它的设计初衷是为了尽可能多的覆盖整个基因组(实际上它的基因组很小)。一方面可能可以当成对照,另一方面有可能会发现新的转录本。但真实原因如何,我就不清楚了。

  9. 您好,看您的博客很受用。由于刚开始学R,很多地方不太清楚,我用R加载.db包时总提示无法打开数据库文件,例如library(org.Hs.eg.db)时,提示
    Error : loadNamespace()里算’org.Hs.eg.db’时.onLoad失败了,详细内容:
    调用: .local(drv, …)
    错误: Could not connect to database:
    unable to open database file
    错误: ‘org.Hs.eg.db’程辑包或名字空间载入失败,
    试了很多方法都没解决,您能指出到底是什么问题,我要怎么做呢!谢谢!

    1. 你试着重新安装一下org.Hs.eg.db.可以使用下面的命令试一下:
      library(BiocInstaller)
      biocLite(“org.hs.eg.db”, type=”source”)

      1. 非常感谢您,这个问题不知道怎么回事,好像在我装了Rstudio之后就没有了。

  10. 今天下载GEO里边的数据用simpleaffy包进行质控分析时,发生错误,麻烦您看看怎么处理?data是我把cel数据下载到本地,用readaffy读入的。
    >library(“affy”)
    >library(“simpleaffy”, lib.loc=”~/R/win-library/3.3″)
    >qc_data<- qc(data)
    Error in setQCEnvironment(cdfn) :
    Could not find array definition file ' hugene10stv1cdf.qcdef '. Simpleaffy does not know the QC parameters for this array type.
    See the package vignette for details about how to specify QC parameters manually.

    1. simpleaffy并不是针对这一芯片设计的,simpleaffy要求有PM和MM,以及spikeIn,但是新的芯片很可能没有MM探针了。

  11. 老师你好,我想请问一下,在有的GEO数据中是没有CEL文件,或者它的RAW文件是txt格式的。如果没有RAW文件或者文件是txt格式的话,是不是不能做芯片的质量分析呢?

    1. 准确地说是不能按照本文的方法来对非Affymetrix的芯片来做质量分析。不同公司的芯片,分析流程也是不一样的。

发表评论

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