摘要
在本章中我们将展示怎么使用Bioconductor中的元数据。获得一个微阵列实验的探针号并把它们映射到靶基因之后，我们想要使用基因的注释和产品来更好的解释实验结果。通常在初步分析中会使用基因注释，这样可以缩窄需要考虑的数据集，改善多重检验的问题，或者发现特定的生物学假设。
8.1 我们的数据
我们使用ALL的一个子集，建议你可以翻阅第一章中对这个数据和随后的代码的描述，在这段代码里面我们加载了该数据集然后选了携带BCR/ABL或者ALL1/AF4易位的肿瘤样本。
library("ALL")
data("ALL")
types = c("ALL1/AF4", "BCR/ABL")
bcell = grep("^B", as.character(ALL$BT))
ALL_af4bcr = ALL[, intersect(bcell, which(ALL$mol.biol %in% types))]
ALL_af4bcr$mol.biol = factor(ALL_af4bcr$mol.biol)
我们应用非特异过滤的方法来移除那些可能不含有信息的探针集。这一步使用genefilter里面的nsFilter函数。**nsFilter函数进行差异筛选的默认参数是IQR。当组间的样本数量相近时，这是一个合理的选择。这与BCR/ABL、ALL1/AF4子集中的例子不一样，其中ALL1/AF4阳性组的样本量非常小。
table(ALL_af4bcr$mol.biol)
ALL1/AF4 BCR/ABL
10 37
在计算IQR中50%的极端值作为离群值会被丢弃，这样测量主要由更加大的BCR/ABL阳性组主导。我们采用例如标准差的非机器方法来测量变异度，但这样会使得离群值更容易被滤过。实际上我们会察看更加极端的分位数的范围，如下：
qrange <- function(x) diff(quantile(x, c(0.1, 0.9)))
library("genefilter")
filt_af4bcr = nsFilter(ALL_af4bcr, require.entrez=TRUE, require.GOBP=TRUE, var.func=qrange, var.cutoff=0.5)
ALLfilt_af4bcr = filt_af4bcr$eset
现在让我们加载必要的工具和注释数据包
library("Biobase")
library("annotate")
library("hgu95av2.db")
我们的第一步使用rowttests进行两组间比较并选择前100个基因。
rt = rowttests(ALLfilt_af4bcr, "mol.biol")
names(rt)
"statistic" "dm" "p.value"
练习 8.1
- 画出t值和p值的直方图，如下图。
练习 8.2
- 创建一个ExpressionSet对象ALLsub，使其含有最小的p值的400个探针集。
有多种不同的变异个可能的改进来完成这个探针集筛选程序，但这不是本章目的。这里的目的是获取一个合理的探针子集并用它来展示元数据的使用。
练习 8.3
- 映射到同一个EntrezID的探针集在ALL中有多少个？在ALLsub中有多少个？
练习 8.4
- 画出CD44基因的表达图，如下图左部。
练习 8.5
- 画出如下图右部的条形图，它显示了ALLsub中的探针映射的基因在每条染色体上的基因数。
在下面的代码框中，我们展示了如何使用annaffy包来产生一列基因的HTML表格，这里是ALLsub中的400个基因。
在Bioconductor2.1版本中，annaffy包不支持使用SQLite技术的hgu95av2.db注释包，使用的是R环境哈希表的旧式的hgu95av2注释包，这个哈希表是我们加载它的目的。在这个数据内容里面hgu95av2.db和hgu95av2是等价的。
library("annaffy")
library("hgu95av2")
anncols = aaf.handler(chip="hgu95av2")[c(1:3, 8:9, 11:13)]
anntable = aafTableAnn(featureNames(ALLsub), "hgu95av2", anncols)
saveHTML(anntable, "ALLsub.html", title="The Features in ALLsub")
localURL = file.path("file:/", getwd(), "ALLsub.html")
在这个文件中我们可以指向HTML浏览器。
browseURL(localURL)
为了继续后面的练习，我们再次hgu95av2注释包（仍然保持hgu95av2.db加载）。
detach("package:hgu95av2")
8.2 每个基因的多个探针集
注释包hgu95av2.db提供了阵列中出现的基因的信息，包括它们的EntrezGene标识符、Unigene簇标识符、基因名、染色体定位、GO注释和通路关系 (Wheeler et al., 2007; Mulder et al., 2007)。尽管基因一词有多个意义，我们用Entrez数据库中的条目来定义它(Maglott et al., 2007)。但是随之一个问题是有些基因出现在一个芯片的多个探针集中。
对这个问题没有一个简单的答案，因为生物学远比微阵列设计的一个基因一个转录本的模型复杂(ENCODE Project Consortium et al., 2007)。在一些例子中，你可能想怎么处理一个探针集对应一个基因的模型，而另一个探针集是另一种模型。在另外一些例子中，有多个探针集对应一个基因，但是只有一个对应另外一个，但这个不应该影响你的推理方法。
练习 8.6
-
选择一些映射到同样的基因的探针集对，画出他们的表达值图，如下图。
EntrezGene7013(TERF1)
8.3 策略和过度呈现
在本章后面我们将考虑使用GO注释数据并试着发现选出的基因的生物学意义(The Gene Ontology Consortium, 2000)。在这一部分我们先处理一个更一般的问题。假设我们你可以把基因分成k组。比如k可是24且代表人染色体位置。对差异表达基因的选择定义了另一个分组：那些含有低p值和另外一些不含的。这两种策略可以用于定义一个二维列联表。在我们的例子中，含有24行2列，每一个列表元素包含有符合分类标准的基因数量。
问题是两种策略（特定染色体位置和是否低表达）是否关联。对列联表有几种不同的检验方法，比如，皮尔逊χ2检验-由R函数chisq.test提供，Fisher's确定性检验-由函数fisher.test函数提供。
总的来说，我们进行分析时需要注意，单个基因可能有多个探针集。尽管注释包中的注释是根据每一个基因给定的，但是多个探针集映射到同一个EntrezGene ID的信息是冗余的。得出的结论可能是误导性的，比如三个冗余信息都被使用了。对于目前的数据ALLsub，这不是问题，因为我们之前使用过nsFilter进行过滤。
这里我们把Entrez ID作为初始的关键词的策略，这样不同类型的注释就可以相互映射。根据感兴趣的问题，需要不同的方法，比如当考虑与蛋白密切相关的注释时。
练习 8.7
- 创建一个含有gene_id、染色体两列的数据框chr，这此数据框中EntrezGene ID和染色体位相互映射。
练习 8.8
- 创建在ALLfilt_af4bcr数据中差异表达基因关联EntrezGene IDs和染色体号映射的二联表（记得使用此前创建的EGsub的向量）。使用函数fisher.test和chisq.test来检验关联性。（你可能需要查询fisher.test函数的参数simulate.p.value的帮助页面）
一旦我们建立了BCR/ABL、ALL1/AF4中差异表达基因和染色体定位，我们尝试确定更明确的关联，比如考虑残值。
一个不常使用的方法是对每一个染色体进行列联表的独立检验。这经常使用超几何样本模型（基因在不在特定染色体中，和它们是否差异表达）。注意这和对每一个染色体的二维列表使用Fisher's检验是一样的。
备用练习
- 用超几何分布对每一个染色体进行检验。这些结果和上面发现的结果对比有何不同？报告对每一个染色体的统计。如果任何一个p值是显著的，是否多于或者少于比你期望的基因数量？
8.3.1 染色体定位
在其它的设定中，我们对基因定位感兴趣，它们周围有什么基因，或者在过表达前根据定位进行分组。
Bioconductor注释包中包含染色体定位的信息在注释对象中以CHRLOC为前缀。起始位置以整数表示，它们的绝对值是转录起始位点的对染色体5'端的距离。基因正义链是正值反义链是负值。
练习 8.9
- 在ALLsub数据中的正义链含有多少个探针集？
另一个与染色体定位相关的有用的概念是染色体区带。这些信息可由MAP注释对象（例如，hgu95av2MAP）获取。在一些例子中，精确的定位未知，给出的是一个范围。一个对GSEA使用这些的例子可见第13章。
8.4 使用GO
GO是一个描述基因的结构化词汇表，根据基因的分子功能、生物学过程和细胞组成(The Gene Ontology Consortium, 2000)。分子功能描述基因在生化层面的作用，但不涉及该作用何时何地发生。生物过程描述基因涉及的生物学客体。细胞组成描述亚细胞定位和大分子复合物。细胞组分包括核内膜和泛素连接酶复合物。
GO词条可以由两种关系联结：is a和part of。比如，核染色体是一个（is a)染色体，核是细胞的一部分(part of)。本体论是定向非循环(DAG)的结构。一个词条的父类是一个DAG中更一般的GO词条，这些由is a和part of关系联结成的链。例如，生物过程词条己糖生物合成有两个父类，己糖代谢和单糖生物合成。为了精确和简明，所有的GO索引都使用一个GO:前缀和七位数标签，比如GO:0008094。
这里有三个基本的任务，一个是定向层次结构、确定父类词条和子类词条，从全图中取出感兴趣的子图。第二个与GO标签相关的映射处理是自然语言描述，第三个是GO词条和基因或者基因产物之间的映射。在此处有关于处理这些任务的软件工具的包描述。
一个基因集的GO图是关于这些基因的GO注释词条和父类词条的联合结果。
用PARENT和CHILDREN映射可处理寻找不同词条的父类和子类词条。为了寻找“GO:0008094”我们可以使用如下代码：
library("GO.db")
as.list(GOMFCHILDREN["GO:0008094"])
结果如图：
我们用词条offspring来指向一个节点的所有的后代（子类、孙类...）。类似的，我们使用词条ancestor来指向一个节点的父类、祖父类...
as.list(GOMFOFFSPRING["GO:0008094"])
结果如图：
8.4.1 功能分析
annotate和GOstats提供了处理GO的很多必要的功能。其它可以考虑统计分析的包是topGO和goTools。
hyperGTest函数计算导入的GO图中所有GO条目的基因的超几何p值。基本想法是简单。考虑一个辨别基因是否差异表达的实验方法。基因集叫做gene universe。那么对于一个特定的GO词条，全基因被分割为被该词条注释和未被注释的基因。你可以认为这是一个桶，所有的基因是桶里的珠子。被该条目注释的基因是黑色珠子，其它是白色的。你可以问在差异表达的基因中黑色的珠子是否特别的高或者全基因的黑色珠子的总频率相一致。如果黑珠子较多，我们可以认为该GO条目和感兴趣的基因列之间存在相关性。
该检验会多次运行，每一个GO策略运行一次。因为基因注释在GO词条中方式，有人会关注怎么处理发生的多重检验。topGO和GOstats都可以用来处理这些问题。更多的细节见在第14章中的条件性检验图和Falcon and Gentleman (2007)。基因富集分析（GSEA）是另一种处理方法，详见第13章。
library("GOstats")
为了检验我们起先定义的全基因并创建一个参数对象来设定分析选项。
affyUniverse = featureNames(ALLfilt_af4bcr)
uniId = hgu95av2ENTREZID[affyUniverse]
entrezUniverse = unique(as.character(uniId))
params = new("GOHyperGParams", geneIds=EGsub, universeGeneIds=entrezUniverse, annotation="hgu95av2", ontology="BP", pvalueCutoff=0.001, conditional=FALSE, testDirection="over")
mfhyper = hyperGTest(params)
调整一个多重检验的合适方法并不简单。检验方法不是独立的，p值与每一个策略的大小相关（比如注释到该策略的基因数），样本分布也不清晰。尽管如此，很多人还是使用一些多重检验调整。我们推荐一个把未调整的p值和评估效果的大小的程序性方法，比如比值比。我们可以画出未调整p值的条形图，如下图。
hist(pvalues(mfhyper), breaks=50, col="mistyrose")
我们可以看见左边的峰值，接近于0，暗示有些策略富集的基因可能是随机的。
练习 8.10
- GO策略看起来明显是过代表的。你可以使用p值设定为0.001（提示：GOHyperGResult的总结方法会有帮助）。这里有模式吗？
练习 8.11
- 使用GOTERM注释对象来检索某些策略的更广泛的描述。
8.5 其它可以获得的注释
Bioconductor的注释包也包含了Entrez ID号和PubMed IDs号。这些可以联合如pm.getabst和pm.abstGrep函数来用关键词自动下载和搜索PubMed摘要。每一个注释包也提供了到PFAM和Prosite号的映射。
8.6 biomaRt
在前面的部分中，我们通过静态注释包使用了Bioconductor中可用的注释信息。在大量的生物数据库如Ensembl和Uniprot (Flicek et al., 2007; UniProt, 2007)中，有大量有效的附加信息。我们可以在线获取这些数据库，在R中，使用包biomaRt。这里我们用biomaRt来查询差异表达基因的3'UTR。这些信息可以用来进行后续的调节序列分析，如microRNA结合位点。我们需要先创建一个Mart类的实例来存储链接到数据库的信息。所有可用的BioMart网页服务可由函数listMarts来展示，函数head返回前面几个项目。
library("biomaRt")
head(listMarts())
结果如图：
我们举Ensembl为例：
mart = useMart("ensembl")
通常BioMart数据库含有对歌数据集，我们可以函数listDatasets来查看可用的数据集。
head(listDatasets(mart))
我们的数据来自人微阵列，所以想使用hsapiens_gene_ensembl集来进行处理，并需要更新Mart对象。
ensembl = useDataset("hsapiens_gene_ensembl", mart=mart)
对于Ensembl数据库，biomaRt对大部分的文物提供一系列便捷的函数。getGene函数用query IDs向量相关基因的查询名字、描述和染色体定位。getGo可以用来获取GO注释，getSequences返回不同类型的序列信息。getSNP和getHomolog对查询SNP数据和把基因号从一个物种映射到另一个非常有效。
练习 8.12
- 用getSequence来获取差异基因的3'UTR序列。查看该函数的主页来学习函数参数。想想我们的基因集有可以获得的哪类基因IDs。
biomaRt可以用不同的方式返回多种类型的数据。为了理解它一般的查询API如何工作，我们需先学习filter和attribute。一个过滤定义了对查询的限制，属性定义了我们想返回的值，比如这些基因的PFAM的IDs。你可以用listFilters来获取过滤列。
head(listFilters(ensembl, group="GENE:"))
还有用listAttributes返回可以获取的属性。
head(listAttributes(ensembl, group="PROTEIN:"))
对于一些BioMart数据库，特别是Ensembl，可以获得很多属性和过滤，你可以用参数group控制上述函数列出数。biomaRt的一般交涉界面由函数getBM提供。
练习 8.13
对于我们的差异表达基因，发现相关的蛋白质结构域。这样的结构域存储在PFAM、Prosite、InterPro数据库。在这些数据库中找出结构域ID。
8.7 注释包的数据库版本
在Bioconductor2.1版本中，基于环境的注释包逐步被新包取代，这新包中数据存储在SQLite数据库中。所有的数据库变体都含有.db后缀，如此hgu95av2注释包变异本叫做hgu95av2.db。每一个替更都确保数据库变本与基于环境的包操作基本一样。但是，在关系型数据库API的许多操作都比之前的简单。
这些新包的基本交界包含在AnnotationDbi包中，它们的vignette提供了更多细节以备查询。每一个.db注释包都输出一些函数用来直接进入数据库，这些函数的名字是胡乱的，这样在每一个包中他们都有不同的名字。在hgu133a.db中这个函数是hgu133a_dbconn。
library("hgu133a.db")
dbc = hgu133a_dbconn()
几乎所有的基于环境的注释包支持的交互界面也被db版本支持。比如你可以直接用get, mget, $, [[.这些并准的取子集或者提取的函数在环境中也可以使用。
> get("201473_at", hgu133aSYMBOL)
[1] "JUNB" mget(c("201473_at","201476_s_at"), hgu133aSYMBOL)
$`201473_at`
[1] "JUNB"
$`201476_s_at`
[1] "RRM1"
> hgu133aSYMBOL$"201473_at"
[1] "JUNB"
> hgu133aSYMBOL[["201473_at"]]
[1] "JUNB"
在接下来的例子中可见看到新范式的某些优势，我们创建一个表格来计数每一个GO策略中的条目数。基于环境的包使用的代码如下（在db中也能使用）：
goCats = unlist(eapply(GOTERM, Ontology))
gCnums = table(goCats)[c("BP","CC", "MF")]
我们可以使用xtable包来展示结果。
library("xtable")
xtable(as.matrix(gCnums), display=c("d", "d"), caption="Number of GO terms per ontology.", label="ta:GOprops")
| x | |
|---|---|
| BP | 12916 |
| CC | 2007 |
| MF | 7878 |
下面的代码使用的是新的db函数，它更快。
query = "select ontology from go_term"
goCats = dbGetQuery(GO_dbconn(), query)
gCnums2 = table(goCats)[c("BP","CC", "MF")]
identical(gCnums, gCnums2)
[1] TRUE
我们可以搜索含chromosome的GO词目，首先建一个SQL询问语句，然后调用它。
query = paste("select term from go_term where term", "like '%chromosome%'")
chrTerms = dbGetQuery(GO_dbconn(), query)
nrow(chrTerms)
[1] 82
head(chrTerms)
下一步我们展示怎么寻找“transcription factor binding”的GO标识符，并用来获取所有有该注释的Entrez Gene IDs。这里我们把注意力放在HG-U133A芯片覆盖的基因。另外你一可以选择使用org.Hs.eg.db包，它注释整个人类基因组。
query = paste("select go_id from go_term where", "term = 'transcription factor binding'")
tfb = dbGetQuery(GO_dbconn(), query)
tfbps = hgu133aGO2ALLPROBES[[tfb$go_id]]
table(names(tfbps))
IDA IEA IEP IGI IMP IPI ISS NAS NR TAS
106 83 2 2 8 59 30 91 27 366
练习 8.14
- 单词transcription factor中有多少GO条目？
8.7.1 映射标识
在这一部分我们探讨一个高级的话题。基于环境的系统的一个问题是探针ID在大多数映射中作为初始关键词，这使得在没有初始关键词的实体之间进行映射比较困难。基于数据的系统比较直接。我们考虑一个更普遍的问题，把基因标识符映射到其它形式的标识。当出版物使用一个基因的符号或者名字而不是系统性的数据库标识时，这个问题就出现了。
下一节的代码由四个函数组成，三个辅函数和一个从标识符映射到Entrez Gene IDs的主函数findEGs。我们需要知道表的结构来写辅函数，这些辅函数基本是对SQL语句的R封装。hgu95av2_dbschema函数可用来获取关于架构的所有信息。
queryAlias = function(x) {
it = paste("('", paste(x, collapse="', '"), "'", sep="")
paste("select _id, alias_symbol from alias",
"where alias_symbol in", it, ");")
}
queryGeneinfo = function(x) {
it = paste("('", paste(x, collapse="', '"), "'", sep="") paste("select _id, symbol from gene_info where", "symbol in", it, ");")
}
queryGenes = function(x) {
it = paste("('", paste(x, collapse="', '"), "'", sep="")
paste("select * from genes where _id in", it, ");")
}
findEGs = function(dbcon, symbols) {
rs = dbSendQuery(dbcon, queryGeneinfo(symbols))
a1 = fetch(rs, n=-1)
stillLeft = setdiff(symbols, a1[,2])
if( length(stillLeft)>0 ) {
rs = dbSendQuery(dbcon, queryAlias(stillLeft))
a2 = fetch(rs, n=-1)
names(a2) = names(a1)
a1 = rbind(a1, a2)
}
rs = dbSendQuery(dbcon, queryGenes(a1[,1]))
merge(a1, fetch(rs, n=-1))
}
我们先查看这些标识符当前是否可用，然后对找不到的搜索同义词表，看是否有更新过的名称。findEGs函数的前两个询问返回标识（表格中的a1和a2列）和一盒内嵌在SQLite数据库（第一列）中的标识号。最后一个对这些标识号询问语句提取了相应的EntrezGene IDs。
findEGs(dbc, c("ALL1", "AF4", "BCR", "ABL"))
| _id | symbol | gene_id | |
|---|---|---|---|
| 1 | 20 | ABL | 25 |
| 2 | 543 | BCR | 613 |
| 3 | 3781 | AF4 | 4299 |
| 4 | 3946 | ABL | 4547 |
结果中的三列分别是内在ID，标识符和EntrezGene ID（gene_id）。
8.7.2 其它的能力
在许多情形中，你可能想反向映射。意思是如果你有从Affymetrix ID到标识符的映射注释，你可能想有从标识符到Affymetrix ID的映射。这可以用revmap得到。
s1 = revmap(hgu133aSYMBOL)
s1$BCR
[1] "202315_s_at" "214623_at" "217223_s_at"
另一个在AnnotationDbi中有用的工具是toTable函数。它用Bimap类的实例作为输入，返回一个data.frame对象。你可以询问AnnotationDbi包的文档来获取更多细节。在下面的代码中用toTable获取GO信息。
toTable(hgu133aGO["201473_at"])
网友评论