Bioconductor基础ExpressionSet介绍 6

12 10月 2012   | 程序员

Tags: ·

在Biobase基础包中，ExpressionSet是非常重要的类，因为Bioconductor设计之初是为了对基因芯片数据进行分析，而ExpressionSet正是Bioconductor为基因表达数据格式所定制的标准。它是所有涉及基因表达量相关数据在Bioconductor中进行操作的基础数据类型，比如affyPLM, affy, oligo, limma, arrayMagic等等。所以当我们学习Bioconductor时，第一个任务就是了解并掌握ExpressionSet的一切。

ExpressionSet的组成：

  • assayData: 一个matrix类型或者environment类型数据。用于保存表达数据值。
    当它是一个matrix时，它的行表示不同的探针组（probe sets）（也是features，总之是一个无重复的索引值）的值，它的列表示不同的样品。如果有行号或者列号的话，那么行号必须与featureData及phenoData中的行号一致，列号就是样品名。当我们使用exprs()方法时，就是调取的这个assayData的matrix。
    当它是一个enviroment时，它必须有两个变量，一个就是与上一段描述一致的matrix，另一个就是epxrs，而这个exprs会响应exprs()方法，返回表达值。
  • 头文件：用于描述实验平台相关的数据，其中包括phenoData, featureData，protocolData以及annotation等等。其中
    phenoData是一个存放样品信息的data.frame或者AnnotatedDataFrame类型的数据。如果有行号的话，其行号必须与assayData的列号一致（也就是样品名）。如果没有行号，则其行数必须与assayData的列数一致。
    featureData是一个存放features的data.frame或者AnnotatedDataFrame类型的数据。它的行数必须与assayData的行数一致。如果有行号的话，那么它的行号必须和assayData的行号一致。
    annotation是用于存放芯片类型的字符串，比如hgu95av2之类。
    protocolData用于存放设备相当的数据。它是AnnotatedDataFrame类型。它的维度必须与assayData的维度一致。
  • experimentData: 一个MIAME类型的数据，它用于保存和实验设计相关的资料，比如实验室名，发表的文章，等等。那么什么是MIAME类呢？MIAME是Minimum Information About a Microarray Experiment的首字母缩写，它包括以下一些属性（slots）：
    1. name: 字符串，实验名称
    2. lab: 字符串，实验室名称
    3. contact: 字符串，联系方式
    4. title: 字符串，一句话描述实验的内容
    5. abstract: 字符串，实验摘要
    6. url: 字符串，实验相关的网址
    7. samples: list类，样品的信息
    8. hybridizations: list类，杂交的信息
    9. normControls: list类，对照信息，比如一些持家基因（house keeping genes）
    10. preprocessing: list类，原始数据的预处理过程
    11. pubMedIds: 字符串，pubMed索引号
    12. others: list类，其它相关的信息

    有了这些，所有实验相关的信息基本全备。

ExpressionSet继承了eSet类，属性基本和eSet保持一致。

那么，对于一个ExpressionSet，哪些属性是必须的？哪些有可能缺失呢？很显然，assayData是必须的，其它的可能会缺失，但是不能都缺失，因为那样的话就无法完成数据分析的工作。

接下来，我们试着开始实例化一个ExpressionSet类。
第一步，得到assayData。

> dataDir <- system.file("extdata", package="Biobase")
> exprsFile <- file.path(dataDir, "exprsData.txt")
> exprsFile
[1] "/Library/Frameworks/R.framework/Versions/2.15/Resources/library/Biobase/extdata/exprsData.txt"
> exprs <- as.matrix(read.table(exprsFile, header=TRUE, sep="\t", row.names=1, as.is=TRUE))
> class(exprs)
[1] "matrix"
> dim(exprs)
[1] 500  26
> colnames(exprs)
 [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q" "R" "S" "T" "U" "V" "W" "X" "Y" "Z"
> head(exprs)
                       A         B        C        D        E       F        G       H       I
AFFX-MurIL2_at  192.7420  85.75330 176.7570 135.5750 64.49390 76.3569 160.5050 65.9631 56.9039
AFFX-MurIL10_at  97.1370 126.19600  77.9216  93.3713 24.39860 85.5088  98.9086 81.6932 97.8015
AFFX-MurIL4_at   45.8192   8.83135  33.0632  28.7072  5.94492 28.2925  30.9694 14.7923 14.2399
AFFX-MurFAS_at   22.5445   3.60093  14.6883  12.3397 36.86630 11.2568  23.0034 16.2134 12.0375
AFFX-BioB-5_at   96.7875  30.43800  46.1271  70.9319 56.17440 42.6756  86.5156 30.7927 19.7183
AFFX-BioB-M_at   89.0730  25.84610  57.2033  69.9766 49.58220 26.1262  75.0083 42.3352 41.1207
                        J        K       L       M       N       O       P        Q        R
AFFX-MurIL2_at  135.60800 63.44320 78.2126 83.0943 89.3372 91.0615 95.9377 179.8450 152.4670
AFFX-MurIL10_at  90.48380 70.57330 94.5418 75.3455 68.5827 87.4050 84.4581  87.6806 108.0320
AFFX-MurIL4_at   34.48740 20.35210 14.1554 20.6251 15.9231 20.1579 27.8139  32.7911  33.5292
AFFX-MurFAS_at    4.54978  8.51782 27.2852 10.1616 20.2488 15.7849 14.3276  15.9488  14.6753
AFFX-BioB-5_at   46.35200 39.13260 41.7698 80.2197 36.4903 36.4021 35.3054  58.6239 114.0620
AFFX-BioB-M_at   91.53070 39.91360 49.8397 63.4794 24.7007 47.4641 47.3578  58.1331 104.1220
                        S       T         U        V       W         X       Y         Z
AFFX-MurIL2_at  180.83400 85.4146 157.98900 146.8000 93.8829 103.85500 64.4340 175.61500
AFFX-MurIL10_at 134.26300 91.4031  -8.68811  85.0212 79.2998  71.65520 64.2369  78.70680
AFFX-MurIL4_at   19.81720 20.4190  26.87200  31.1488 22.3420  19.01350 12.1686  17.37800
AFFX-MurFAS_at   -7.91911 12.8875  11.91860  12.8324 11.1390   7.55564 19.9849   8.96849
AFFX-BioB-5_at   93.44020 22.5168  48.64620  90.2215 42.0053  57.57380 44.8216  61.70440
AFFX-BioB-M_at  115.83100 58.1224  73.42210  64.6066 40.3068  41.82090 46.1087  49.41220

第二步，准备头文件。

> #准备phenoData
> pDataFile <- file.path(dataDir, "pData.txt")
> pDataFile
[1] "/Library/Frameworks/R.framework/Versions/2.15/Resources/library/Biobase/extdata/pData.txt"
> pData <- read.table(pDataFile, row.names=1, header=TRUE, sep="\t")
> dim(pData)
[1] 26  3
> rownames(pData)
 [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q" "R" "S" "T" "U" "V" "W"
[24] "X" "Y" "Z"
> summary(pData)
    gender        type        score       
 Female:11   Case   :15   Min.   :0.1000  
 Male  :15   Control:11   1st Qu.:0.3275  
                          Median :0.4150  
                          Mean   :0.5369  
                          3rd Qu.:0.7650  
                          Max.   :0.9800  
> head(pData)
  gender    type score
A Female Control  0.75
B   Male    Case  0.40
C   Male Control  0.73
D   Male    Case  0.42
E Female    Case  0.93
F   Male Control  0.22
> #检查pData的行号是否和exprs的列号一致
> all(rownames(pData) == colnames(exprs))
[1] TRUE
> #然后将pData转换成一个AnnotatedDataFrame实例
> metadata <- data.frame(labelDescription=c("Patient gender", "Case/control status", "Tumor progress on XYZ scale"), row.names=colnames(pData))
> adf <- new("AnnotatedDataFrame", data=pData, varMetadata=metadata)
> adf
An object of class 'AnnotatedDataFrame'
  rowNames: A B ... Z (26 total)
  varLabels: gender type score
  varMetadata: labelDescription
> dim(adf)
   rowNames columnNames 
         26           3 
> varMetadata(adf)
                  labelDescription
gender              Patient gender
type           Case/control status
score  Tumor progress on XYZ scale
> dimLabels(adf)
[1] "rowNames"    "columnNames"
> featureNames(adf)
 [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q" "R" "S" "T" "U" "V" "W"
[24] "X" "Y" "Z"
> sampleNames(adf)
 [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q" "R" "S" "T" "U" "V" "W"
[24] "X" "Y" "Z"
> head(pData(adf))
  gender    type score
A Female Control  0.75
B   Male    Case  0.40
C   Male Control  0.73
D   Male    Case  0.42
E Female    Case  0.93
F   Male Control  0.22
> #准备annotation
> annotation <- "hgu95av2"

第三步，准备experimentData。

> experimentData = new("MIAME", name="Pierre Fermat",
+ lab="Francis Galton Lab",
+ contact="pfermat@lab.not.exist",
+ title="Smoking-Cancer Experiment",
+ abstract="An example ExpressionSet",
+ url="www.lab.not.exist",
+ other=list(notes="Created from text files"))
> experimentData
Experiment data
  Experimenter name: Pierre Fermat 
  Laboratory: Francis Galton Lab 
  Contact information: pfermat@lab.not.exist 
  Title: Smoking-Cancer Experiment 
  URL: www.lab.not.exist 
  PMIDs:  
 
  Abstract: A 3 word abstract is available. Use 'abstract' method.
  notes:
   notes:     
      Created from text files

最后，把它们整合起来。

> exampleSet = new("ExpressionSet", exprs=exprs,
+ phenoData=adf, experimentData=experimentData,
+ annotation="hgu95av2")
> exampleSet
ExpressionSet (storageMode: lockedEnvironment)
assayData: 500 features, 26 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: A B ... Z (26 total)
  varLabels: gender type score
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu95av2

以上，我们一步一步地实例了一个ExpressionSet类。但是在实际的Bioconductor操作中是很少从第一步开始这样实例一个ExpressionSet类的。最常见的情况是我们使用affy中的ReadAffy()函数，读取.CEL文件进入R成为一个AffyBatch实例，而后使用rma()或者mas5()函数对样品进行处理之后得到一个ExpressionSet的实例。我们之所以做上面的练习，只是为了能够更好地理解什么是ExpressionSet类。

接下来需要了解的是ExpressionSet有哪些方法，我们如何操作ExpressionSet类。

> #可以直接使用$符号来访问AnnotatedDataFrame型的数据
> exampleSet$gender[1:5]
[1] Female Male   Male   Male   Female
Levels: Female Male
> #使用pData函数来获得phenoData中的数据
> head(pData(exampleSet))
  gender    type score
A Female Control  0.75
B   Male    Case  0.40
C   Male Control  0.73
D   Male    Case  0.42
E Female    Case  0.93
F   Male Control  0.22
> #使用featureNames来读取features.我们注意到实际上exampleSet的featureData是空的，但是我们仍然可以得到featureNames。
> featureNames(exampleSet)[1:5]
[1] "AFFX-MurIL2_at"  "AFFX-MurIL10_at" "AFFX-MurIL4_at"  "AFFX-MurFAS_at"  "AFFX-BioB-5_at" 
> #得到样品名称
> sampleNames(exampleSet)[1:5]
[1] "A" "B" "C" "D" "E"
> #使用exprs()函数来提取表达值，也就是assayData中的matrix。
> mat <- exprs(exampleSet)
> dim(mat)
[1] 500  26
> #使用phenoData()函数来获取整个phenoData。注意它和pData()函数的区别。
> adf <- phenoData(exampleSet)
> adf
An object of class 'AnnotatedDataFrame'
  sampleNames: A B ... Z (26 total)
  varLabels: gender type score
  varMetadata: labelDescription

对于ExpressionSet最重要的操作就是如何取出子集了。有时候在进行质量分析之后，我们对其中一些样品的数据不满意，想从已经实例化的ExpressionSet中抽取掉，或者我们希望对样品进行分组，都需要使用到Subset的概念。那么如何抽取子集呢？
我们可以象操作矩阵那样对其进行子集操作。

> vv <- exampleSet[1:5, 1:3]
> vv
ExpressionSet (storageMode: lockedEnvironment)
assayData: 5 features, 3 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: A B C
  varLabels: gender type score
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu95av2 
> featureNames(vv)
[1] "AFFX-MurIL2_at"  "AFFX-MurIL10_at" "AFFX-MurIL4_at"  "AFFX-MurFAS_at"  "AFFX-BioB-5_at"

我们还可以使用它的一些属性来对其进行子集操作。需要注意的是行列的位置。

> males <- exampleSet[, exampleSet$gender == "Male"]; males
ExpressionSet (storageMode: lockedEnvironment)
assayData: 500 features, 15 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: B C ... X (15 total)
  varLabels: gender type score
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu95av2 
> dim(males)
Features  Samples 
     500       15 
> AFFX <- exampleSet[grepl("^AFFX", featureNames(exampleSet)), ]
> AFFX
ExpressionSet (storageMode: lockedEnvironment)
assayData: 67 features, 26 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: A B ... Z (26 total)
  varLabels: gender type score
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu95av2 
> dim(AFFX)
Features  Samples 
      67       26

如果想了解更多，请在R中使用?ExpressionSet以及?eSet来获得更多的资料。

6 thoughts on “Bioconductor基础ExpressionSet介绍”

  1. Reply Joe 12月 19,2012 6:52 上午

    哈哈 经历了几个月的严格封网 终于又翻出来了 读了您的几篇博文 爽呀！
    有两个小问题请教，1. 在表达谱芯片中，处理实验组和对照组样本时，是将对照组的值减去呢？还是除呢？ heatmap中，一般是展现norm intensity还是FC？
    2.为啥heatmap做出来的图，红色表示负值，绿色代表正值呀？好像反了。。

    • Reply admin 12月 19,2012 3:16 下午

      感谢你的回访。
      1.一般是绘制norm intensity，也有可能绘制的是log2转换的intensity， 也有可能是fold change。这主要是看你的目的。如果你想表现不同样品之间的变化的话，你一定会将x轴设定成你的所有的样品，并且将它们分组集中起来，而你的y轴将是一些基因或者探针组。这样的图生成出来以后有一定的宽度，比较好看。但是如果你的篇辐有限，并且注重表现多组样品相对于参照的变化的话，可以使用fc.但是这是很少见的。

      2.如果颜色相反的话，可以重新设置一下颜色值。当初我的教程可能没有注意这个问题，你把它反过来设置一下就好了。比如：pheatmap(data, color = colorRampPalette(c(“green”, “white”, “red”))(50))中如果颜色和你想象的相反的话，你就改写成pheatmap(data, color = colorRampPalette(c(“red”, “white”, “blue”))(50))就好了。

  2. Reply kain 10月 30,2013 9:35 下午

    非常感谢您的这一系列的文章，真的是帮助非常非常的大！感激不尽啊！

  3. Reply bob 10月 29,2014 2:54 上午

    Error in getClass(Class, where = topenv(parent.frame())) : “AnnotatedDataFrame” is not a defined class.

Leave a Reply

  

  

  

%d 博主赞过：