1. 介绍
2. 从头构建一个 ExpressionSet 数据
3. ExpressionSet 的一些基本特性
- 查看 ExpressionSet 的各个元素
- 取子集

本文翻译自 An Introduction to Bioconductor’s ExpressionSet Class ，有删改。

接触生物信息学和 Bioconductor 的话，ExpressionSet 数据类型是一个很基础的概念。之前我也是一直下了文档但是没仔细看。

1. 介绍

ExpressionSet 是 Biobase 包提供的。Biobase 本身是 Bioconductor 项目的一个部分，为很多基因组学数据提供数据类型支持。ExpressionSet 就是用来专门将多种不同来源的数据整合到一起方便数据处理的一种数据类型，很多 Bioconductor 函数输入输出的都是 ExpressionSet 类型的数据。

简单的说，ExpressionSet 把表达数据（assayData 存储芯片、测序等表达数据），表型信息（phenoData 存储样本信息），注释等元信息（featureData, annotation 存储芯片或者测序技术的元数据），以及操作流程（protocolData 存储样本处理相关的信息，通常由厂家提供）和实验（experimentData 用来描述实验相关信息）几种密切相关数据封装在一起，这样我们处理数据的时候不用关心各个数据的细节把它们当作一个整体来看就好。

下面来看怎么一步一步的构建一个 ExpressionSet 数据。

2. 从头构建一个 ExpressionSet 数据

如果我们在 R 里用 affyPLM, affy, oligo 或者 limma 这些包读入了芯片的原始 CEL 数据的话，比如 affy 包的 ReadAffy(), expresso() 和 justRMA() 很可能读出来的数据已经就是 ExpressionSet 类型了。如果不是的话，我们可以用 convert 包的 as 来转换得到（object 是我们要转换的数据）：

library("convert")
as(object, "ExpressionSet")

我们重点来看怎么从头构建一个 ExpressionSet 数据。

开头说过，芯片或者其他一些高通量基因组学技术通常产生几种相关的数据：assayData，phenoData，featureData，protocolData 和 experimentData。要构建 ExpressionSet 数据，我们就需要按照这几个类型数据一个个准备好数据，然后按照一定的规则组合在一起就行了。

Assay data

首先是 assayData，这个数据通常就是我们最关新的存储表达数据的表格。这个数据通常是一个 F (features) 行 x S (samples) 列的表格数据。每一行代表一个 feature（可能是探针、基因、转录本等等），每一列代表一个样本。我们这里就用 Biobase 自带的数据做例子：

library("Biobase")
dataDirectory <- system.file("extdata", package = "Biobase")
exprsFile <- file.path(dataDirectory, "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"
## [20] "T" "U" "V" "W" "X" "Y" "Z"
head(exprs[, 1:5])
##                        A         B        C        D        E
## AFFX-MurIL2_at  192.7420  85.75330 176.7570 135.5750 64.49390
## AFFX-MurIL10_at  97.1370 126.19600  77.9216  93.3713 24.39860
## AFFX-MurIL4_at   45.8192   8.83135  33.0632  28.7072  5.94492
## AFFX-MurFAS_at   22.5445   3.60093  14.6883  12.3397 36.86630
## AFFX-BioB-5_at   96.7875  30.43800  46.1271  70.9319 56.17440
## AFFX-BioB-M_at   89.0730  25.84610  57.2033  69.9766 49.58220

现在我们就可以用 exprs 构建一个初始的 ExpressionSet 数据了：

minimalSet <- ExpressionSet(assayData = exprs)

但是要利用 ExpressionSet 数据类型丰富的特征，我们还要为这个数据添加表型、样本等数据，我们下面继续。

Phenotypic data

表型数据是描述样本相关信息的，比如受试者的性别、年龄，所在分组等等，这些通常也称为变量。表型数据通常是 S (samples) 行 x V (variables) 列的表格。每一行代表一个样本，比如一个受试者、一个组织样本或者一瓶细胞等等，列是一个关于样本的变量。还是拿 Biobase 自带的数据来看看：

pDataFile <- file.path(dataDirectory, "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"
## [20] "T" "U" "V" "W" "X" "Y" "Z"
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
summary(pData)
##     gender              type               score       
##  Length:26          Length:26          Min.   :0.1000  
##  Class :character   Class :character   1st Qu.:0.3275  
##  Mode  :character   Mode  :character   Median :0.4150  
##                                        Mean   :0.5369  
##                                        3rd Qu.:0.7650  
##                                        Max.   :0.9800

示例数据有 26 行 3 列，这个 26 行刚好和前面的表达数据 exprs 的 26 列是一样的。如果我们看一下行列名的话，会发现行列名连顺序都一模一样：

all(rownames(pData) == colnames(exprs))
## [1] TRUE

assayData 和 pData 之间的这种对应关系其实是必须的，如果行列数或者行列名不对应的话，ExpressionSet 会报错。

pData 的每一列都是一个变量，不同的变量可能有不同的数据类型。比如年龄是数值，性别是字符，肿瘤的分期可能是因子等等。读入数据的时候要注意这些列数据类型是否正确，read.table() 的 colClasses 参数可以用来在读入数据的时候设置每一列的数据类型。

有的时候我们可能觉得 pData 里单独的列名不够明确，比如这里的 type 和 score 都不知道具体指的什么。这时候我们还可以通过一个数据框提供额外的元数据：

metadata <- data.frame(
  labelDescription =
    c("Patient gender",
      "Case/control status",
      "Tumor progress on XYZ scale"),
  row.names = c("gender", "type", "score"))

我们这里提供的元数据只有一列 labelDescription，这也是必须要有的一列，想要其他列还可以自己按需添加。

Biobase 提供了一个 AnnotatedDataFrame 数据类型把 pData 和 metadata 封装起来：

phenoData <- new("AnnotatedDataFrame", data = pData, varMetadata = metadata)
phenoData
## An object of class 'AnnotatedDataFrame'
##   rowNames: A B ... Z (26 total)
##   varLabels: gender type score
##   varMetadata: labelDescription

对于一个 AnnotatedDataFrame 类型数据，有一些很实用的函数用来查看该数据中存储的不同的数据：sampleNames()，pData() 和 varMetadata() 等等。另外，AnnotatedDataFrame 能像一个普通的数据框那样取子集：

head(pData(phenoData))
##   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
phenoData[c("A","Z"),"gender"]
## An object of class 'AnnotatedDataFrame'
##   rowNames: A Z
##   varLabels: gender
##   varMetadata: labelDescription
pData(phenoData[phenoData$score>0.8,])
##   gender    type score
## E Female    Case  0.93
## G   Male    Case  0.96
## X   Male Control  0.98
## Y Female    Case  0.94

Annotations and feature data

feature 的元数据，通常我们叫注释数据，和样本的元数据一样重要，因为我们最终得到各种结果都要通过注释落到具体的基因上去，不然得到的结果是一堆探针这是没用的。通常来说很多不同的实验组会用一种芯片，这个很好理解，因为芯片型号就那么多大家都用，重复平台是很常见的。所以呢，每个数据都把芯片的注释文件打包一份进去效率就太低了。所以，通常芯片的注释文件都是单独的 Bioconductor 包。这些注释文件描述了探针与基因、基因的功能等等之间的对应关系，有时候还有有 GO 和 KEGG 等其他来源的信息。annotate 和 AnnotationDbi 包就是用来处理这些注释元数据包的。

由于注释信息由专门的包来提供，所以我们只需要提供平台信息作为注释信息就够了。我们这个数据是 Affymetrix hgu95av2 芯片的：

annotation <- "hgu95av2"

Experiment description

通常关于实验，我们可以提供研究者和实验室信息、联系方式、研究题目等等相关的信息。MIAME(Minimum Information About a Microarray Experiment) 对象可以用来这种信息，我们通过 new() 来创建一个 MIAME 对象：

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")
)

new() 函数通常接受一个类型参数，然后后面是多个成对的类名及其对应的内容存储在各自的类别里。关于 MIAME 数据的更多介绍请参考文档。

Assembling an ExpressionSet

现在 ExpressionSet 所需要的各个部分都准备好了，我们就可以用 ExpressionSet() 创建一个 ExpressionSet 对象：

exampleSet <- ExpressionSet(
  assayData = exprs,
  phenoData = phenoData,
  experimentData = experimentData,
  annotation = annotation
)

exampleSet 和前面的 minimalSet 都是 ExpressionSet 对象，但是前者包含了远比后者更丰富的信息。

3. ExpressionSet 的一些基本特性

从头构建了一个 ExpressionSet 数据后，我们现在来看看这个数据的一些特征。help("ExpressionSet-class") 有一个简单地介绍随时可以查看。

如果我们单纯打印一下 ExpressionSet 数据的话，内容可能很长，所以默认行为只会打印出简单地概要信息：

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 的各个元素

首先，直接用 $ 就能查看 ExpressionSet 里的 phenoData（AnnotatedDataFrame 对象）：

exampleSet$gender[1:5]
## [1] "Female" "Male"   "Male"   "Male"   "Female"
exampleSet$gender[1:5] == "Female"
## [1]  TRUE FALSE FALSE FALSE  TRUE

featureNames()，sampleNames() 和 varLabels 一看名字就知道是干嘛的：

featureNames(exampleSet)[1:5]
## [1] "AFFX-MurIL2_at"  "AFFX-MurIL10_at" "AFFX-MurIL4_at"  "AFFX-MurFAS_at" 
## [5] "AFFX-BioB-5_at"
sampleNames(exampleSet)[1:5]
## [1] "A" "B" "C" "D" "E"
varLabels(exampleSet)
## [1] "gender" "type"   "score"

最后，exprs() 可以把表达矩阵取出来：

mat <- exprs(exampleSet)
dim(mat)
## [1] 500  26

取子集

对 ExpressionSet 数据取子集可能是最常见的操作了。取子基的操作和针对一个表达矩阵取子集类似，行参数针对 feature 而列参数针对 sample。比如我们对上面的数据取子集得到一个新的只有 5 个探针和 3 个样本的新的 ExpressionSet 数据：

vv <- exampleSet[1:5, 1:3]
dim(vv)
## Features  Samples 
##        5        3
featureNames(vv)
## [1] "AFFX-MurIL2_at"  "AFFX-MurIL10_at" "AFFX-MurIL4_at"  "AFFX-MurFAS_at" 
## [5] "AFFX-BioB-5_at"
sampleNames(vv)
## [1] "A" "B" "C"

比如我们只要来自样本为男性的数据：

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

Bioconductor 中的 ExpressionSet 数据类型

文章目录