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

教程数据下载

本教程所使用到的数据可以从下面的地址下载:

http://bioinformatics.picr.man.ac.uk/experiments/index.html

并建立一个以mcf7.txt为文件名的数据说明,每一列以tab为间隔。

Factor Channel
MCF7_r1.CEL MCF7standard 1
MCF7_r2.CEL MCF7standard 2
MCF7_r3.CEL MCF7standard 3
a1MCF7_r1.CEL MCF7amplified1ng 1
a1MCF7_r2.CEL MCF7amplified1ng 2
a1MCF7_r3.CEL MCF7amplified1ng 3
a10MCF7_r1.CEL MCF7amplified 10ng 1
a10MCF7_r2.CEL MCF7amplified 10ng 2
a10MCF7_r3.CEL MCF7amplified 10ng 3
a100MCF7_r1.CEL MCF7amplified 100ng 1
a100MCF7_r2.CEL MCF7amplified 100ng 2
a100MCF7_r3.CEL MCF7amplified 100ng 3

您也可以从本站提供的地址下载:数据下载

下面是数据读取的过程:

##变更工作目录
> workingdir <- getwd() > workingdir
[1] "C:/Documents and Settings/jianhong ou/My Documents"
> subdir = "/4R/downloads/MCF7"
> workingdir <- paste(workingdir,subdir,sep="") > setwd(workingdir)
> getwd()
[1] "C:/Documents and Settings/jianhong ou/My Documents/4R/downloads/MCF7"
##调入affy及simpleaff库
> library(affy)
Loading required package: Biobase

Welcome to Bioconductor

  Vignettes contain introductory material. To view, type
  'openVignette()'. To cite Bioconductor, see
  'citation("Biobase")' and for packages 'citation(pkgname)'.

> library(simpleaffy)
Loading required package: genefilter
Loading required package: gcrma

##结构化读入数据
> pd <- read.AnnotatedDataFrame("mcf7.txt",header=TRUE, sep="\t",row.names=1) > pData(pd)
                            Factor Channel
MCF7_r1.CEL           MCF7standard       1
MCF7_r2.CEL           MCF7standard       2
MCF7_r3.CEL           MCF7standard       3
a1MCF7_r1.CEL   MCF7amplified100ng       1
a1MCF7_r2.CEL   MCF7amplified100ng       2
a1MCF7_r3.CEL   MCF7amplified100ng       3
a10MCF7_r1.CEL  MCF7amplified 10ng       1
a10MCF7_r2.CEL  MCF7amplified 10ng       2
a10MCF7_r3.CEL  MCF7amplified 10ng       3
a100MCF7_r1.CEL  MCF7amplified 1ng       1
a100MCF7_r2.CEL  MCF7amplified 1ng       2
a100MCF7_r3.CEL  MCF7amplified 1ng       3
> Data <- ReadAffy(filenames=rownames(pData(pd)),phenoData=pd,verbose = TRUE) 1 reading MCF7_r1.CEL ...instantiating an AffyBatch (intensity a 506944x12 matrix)...done. Reading in : MCF7_r1.CEL Reading in : MCF7_r2.CEL Reading in : MCF7_r3.CEL Reading in : a1MCF7_r1.CEL Reading in : a1MCF7_r2.CEL Reading in : a1MCF7_r3.CEL Reading in : a10MCF7_r1.CEL Reading in : a10MCF7_r2.CEL Reading in : a10MCF7_r3.CEL Reading in : a100MCF7_r1.CEL Reading in : a100MCF7_r2.CEL Reading in : a100MCF7_r3.CEL > class(Data)
[1] "AffyBatch"
attr(,"package")
[1] "affy"
##尝试调用原始数据
> PM = pm(Data)
trying URL 'http://bioconductor.org/packages/2.6/data/annotation/bin/windows/contrib/2.11/hgu133acdf_2.6.0.zip'
Content type 'application/zip' length 1740408 bytes (1.7 Mb)
opened URL
downloaded 1.7 Mb

package 'hgu133acdf' successfully unpacked and MD5 sums checked

The downloaded packages are in
        C:\Documents and Settings\jianhong ou\Local Settings\Temp\Rtmpu6Y8xd\downloaded_packages
> MM = mm(Data)
> gnames <- geneNames(Data) > dim(PM)
[1] 247965     12
> PM[1:5,]
       MCF7_r1.CEL MCF7_r2.CEL MCF7_r3.CEL a1MCF7_r1.CEL a1MCF7_r2.CEL a1MCF7_r3.CEL
129340      1215.3      1153.5      1294.5         131.5         135.8         194.5
213420      3075.3      3185.3      3609.3         123.3         156.5         230.5
396671      2875.0      3024.5      3226.3         128.3         147.3         233.0
82246       8715.0      7934.3      8432.5         496.0         727.8         950.5
430968      6943.3      6791.0      7058.3         476.3         502.3         741.0
       a10MCF7_r1.CEL a10MCF7_r2.CEL a10MCF7_r3.CEL a100MCF7_r1.CEL a100MCF7_r2.CEL
129340          386.3          499.5          464.3           596.0           539.5
213420          774.0          745.5          480.0           827.5           691.5
396671         1086.0          817.5          693.5           901.8           983.5
82246          4698.8         4117.5         3833.3          5800.3          4942.0
430968         4809.5         3702.3         2968.3          4888.8          4124.0
       a100MCF7_r3.CEL
129340           433.8
213420           559.0
396671           705.3
82246           3393.5
430968          3367.3
> MM[1:5,]
       MCF7_r1.CEL MCF7_r2.CEL MCF7_r3.CEL a1MCF7_r1.CEL a1MCF7_r2.CEL a1MCF7_r3.CEL
130052       366.0       386.3       404.3         119.3         138.5         190.0
214132      1605.0      1890.3      2016.8         110.3         192.8         228.0
397383      1787.5      2042.5      2255.8         219.3         269.8         376.3
82958       6588.5      6038.8      6461.0         576.0         812.5        1110.8
431680      3977.3      3734.5      4156.8         281.0         400.0         546.3
       a10MCF7_r1.CEL a10MCF7_r2.CEL a10MCF7_r3.CEL a100MCF7_r1.CEL a100MCF7_r2.CEL
130052          311.5          374.3          344.3           443.0           395.0
214132          755.3          735.8          577.8           926.8           636.3
397383         1676.0         1369.8         1118.0          1640.3          1489.0
82958          4678.0         4587.5         3889.5          6310.5          5193.3
431680         3359.3         2705.0         2158.0          3075.3          2552.5
       a100MCF7_r3.CEL
130052           355.8
214132           573.3
397383          1126.3
82958           3673.0
431680          2131.0
> dim(gnames)
NULL
> gnames[1:5]
[1] "1007_s_at" "1053_at"   "117_at"    "121_at"    "1255_g_at"

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的芯片来做质量分析。不同公司的芯片,分析流程也是不一样的。

发表评论

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