在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）：
- name: 字符串，实验名称
- lab: 字符串，实验室名称
- contact: 字符串，联系方式
- title: 字符串，一句话描述实验的内容
- abstract: 字符串，实验摘要
- url: 字符串，实验相关的网址
- samples: list类，样品的信息
- hybridizations: list类，杂交的信息
- normControls: list类，对照信息，比如一些持家基因（house keeping genes）
- preprocessing: list类，原始数据的预处理过程
- pubMedIds: 字符串，pubMed索引号
- 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来获得更多的资料。
哈哈 经历了几个月的严格封网 终于又翻出来了 读了您的几篇博文 爽呀！
有两个小问题请教，1. 在表达谱芯片中，处理实验组和对照组样本时，是将对照组的值减去呢？还是除呢？ heatmap中，一般是展现norm intensity还是FC？
2.为啥heatmap做出来的图，红色表示负值，绿色代表正值呀？好像反了。。
感谢你的回访。
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))就好了。
非常感谢您的这一系列的文章，真的是帮助非常非常的大！感激不尽啊！
Error in getClass(Class, where = topenv(parent.frame())) : “AnnotatedDataFrame” is not a defined class.
adf <- new("AnnotatedDataFrame", data=pData, varMetadata=metadata)
编译运行这行代码的时候提示上述错误，这是什么原因啊？
请导入library(Biobase)