发表于 | 分类于
| 字数统计: 1,455 | 阅读时长 ≈ 6

前言

我们已经查看了dfHclust 的一些代码。它在很大程度上主要是用于data.frame数据类实例的通用分析：

  • 行可以表示要聚类的对象（就是样本名）；
  • 列可以表示用于比较和分组对象的特征（例如表达值）；
  • rownames 可以为对象添加标签；
  • colnames 可以给特征值添加名称（例如基因名）。

使用名称这种思路可以很容易在app中方便地查看各种基因的表达值。

当我们对基因组级数据进行分析时就有点不太一样了。我们现在将重点放在表达矩阵方面，它有以下特点：

  • 数据通量特别大；
  • 假设我们感兴趣的对象是单个生物学样本，它们已经进行了排列(arrayed)，那么对它们添加棱可能就比较复杂。原文如下：

assuming the objects of interest are the individual biological samples that have been arrayed, the labeling of objects can be complex

  • 假设我们使用ExpressionSet来管理数据，那么rownames colnames 都将倾向于使用无意义的词汇表来进行标记，例如探针数据标识符（它们都是一些数字）和样本标识符（都是一些表示样本名的英文字符，但无法直接看出是什么样本）。

通常来说，这些都是我们最后遇到的麻烦（除了使用tissuesGeneExpression型数据），这里我们有一个比较简单的解决䢍。我们先来看数据量的问题，我们可以通过使用基于括号的子集(subsetting)来方便地管理ExpressionSet数据。

一个简单的app

现在让我们编写一个函数，这个函数可以将ExpressionSet转换为data.frame，然后使用dfHclust()来处理data.frame。这交进我们的app的一个功能，它用于实验数据中的交互式聚类，如下所示：

1
2
3
4
5
6
7
8
9
10
11
library(Biobase)
library(hgu133a.db)
## 
library(ph525x)
esHclust = function(es) {
  emat = t(exprs(es))
  rownames(emat) = sampleNames(es)
  dd = data.frame(emat)
  dfHclust(dd)
}

准备组织特异性表达数据

现在我们转换tissuesGeneExpression包中的matrix/data.frame数据，将它们转换为已经注释了的ExpressionSet数据类型，如下所示：

1
2
3
4
5
6
7
8
9
10
11
library(tissuesGeneExpression)
library(AnnotationDbi)
data(tissuesGeneExpression)
tgeES = ExpressionSet(e)
annotation(tgeES) = "hgu133a.db"
pData(tgeES) = tab
featureNames(tgeES) =
  make.names(mapIds(hgu133a.db, keys=featureNames(tgeES),
    keytype="PROBEID", column="SYMBOL"), unique=TRUE)
## 'select()' returned 1:many mapping between keys and columns
sampleNames(tgeES) = make.names(tgeES$Tissue, unique=TRUE)

以下代码是我们随机挑出一些样本的基因集来放到app中运行

1
2
set.seed(1234)
esHclust( tgeES[1:50, sample(1:ncol(tgeES), size=40) ] )

探索组织分化的聚类分析

这个数据集中的哪些基因对于分组来说最为重要？能找到这些基因的分析工具能否帮助我们找到，并理解这些基因呢？这两个问题都不是特别清楚，我们需要理解构成组织基因表达背后的实验设计，并且完全理解它们。然而，以下代码可以帮助我们确定那些在统计学意义不太可能在所有组织中保持恒定表达的基因。我们使用调和F检验(moderated F test)来计算这些基因，它的零假设是，在所有的组织中，基因的表达是恒定的，代码如下所示：

1
2
3
4
library(limma)
mm = model.matrix(~Tissue, data=pData(tgeES))
f1 = lmFit(tgeES, mm)
ef1 = eBayes(f1)

现在我们根据简单F检验确定了50个最有差异的基因，现在用热图画出来：

1
heatmap(ef1$coef[order(ef1$F,decreasing=TRUE)[1:50],-1], cexCol=.8)

plot of chunk lkheat

以上是相对于参考分类（小脑，cerebellum）的平均表达值进行的可视化。我们也可以将原始的数据进行可视化，如下所示：

1
2
sig50 = rownames(ef1$coef[order(ef1$F,decreasing=TRUE)[1:50],])
heatmap(exprs(tgeES[sig50,]))

plot of chunk lkh2

这都是非常非正式的，只使用热图演示的默认值。我们会继续这样做。在检查平均热图时，我认为以下五个基因是组织分化的潜在特征：

这是非正式地展示数据的方法，我们只使用了默认值。 我认为以下五个基因代表了组织分化的潜在特征，如下所示：

原文为：

This is all very informal, using only default values for heatmap presentations. We’ll continue in this vein. In my inspection of the heatmap of means, I considered the following five genes to be a potential signature for tissue differentiation:

1
sig5 = c("IL26", "ZNF674", "UBC.1", "C7orf25.1", "RPS13")

In the exercises we’ll see whether this is at all satisfactory.

机器学习方法评估区分能力

使用机器学习并不是直接解决可视化的问题，但是对于思考我们可以用可视化做哪些事情，它简单且有用。 MLInterfaces 包可以更容易地使用各种基于R的统计学习工具来处理ExpressionSet 实例。我们使用随机森林来检验一下上面定义的sig50 的有效程度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
library(MLInterfaces)
library(randomForest)
set.seed(1234)
rf1 = MLearn(Tissue~., tgeES[sig50,], randomForestI, xvalSpec("NOTEST"))
RObject(rf1)
## 
## Call:
##  randomForest(formula = formula, data = trdata) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 7
## 
##         OOB estimate of  error rate: 9.52%
## Confusion matrix:
##             cerebellum colon endometrium hippocampus kidney liver placenta
## cerebellum          37     0           0           0      1     0        0
## colon                0    31           1           0      1     1        0
## endometrium          0     4           9           0      2     0        0
## hippocampus          0     0           0          30      0     1        0
## kidney               0     1           0           0     36     2        0
## liver                0     0           0           0      0    26        0
## placenta             0     2           0           0      1     1        2
##             class.error
## cerebellum   0.02631579
## colon        0.08823529
## endometrium  0.40000000
## hippocampus  0.03225806
## kidney       0.07692308
## liver        0.00000000
## placenta     0.66666667

这里我们使用了默认值，我们可以看到50个基因标签能够对这几个组织进行了合理的区分。通过使用MLearn对单行代码进行了简单的修改，我们就可以评估这些基因标签改变后以及计算过程的效果。

参考资料

  1. A shiny app for clustering using genes in ExpressionSet

发表于 | 分类于
| 字数统计: 1,662 | 阅读时长 ≈ 8

NGS的可视化

我们将用以下的包来说明一下如何对NGS数据进行可视化。

1
2
3
4
5
6
#biocLite("pasillaBamSubset")
#biocLite("TxDb.Dmelanogaster.UCSC.dm3.ensGene")
library(pasillaBamSubset)
library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
fl1 <- untreated1_chr4()
fl2 <- untreated3_chr4()

我们将会采用四种方式来查看NGS的数据，分别是IGV，plot,Gviz和ggbio。

IGV

将这些文献从R library的目录中复制到当前的工作目标。首先将工作目录设置为源文件位置。我们需要使用Rsamtools包来对BAM文件加上索引，然后使用IGV来查看，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
file.copy(from=fl1,to=basename(fl1))
## [1] FALSE
file.copy(from=fl2,to=basename(fl2))
## [1] FALSE
library(Rsamtools)
## Loading required package: Biostrings
## Loading required package: XVector
indexBam(basename(fl1))
##       untreated1_chr4.bam 
## "untreated1_chr4.bam.bai"
indexBam(basename(fl2))
##       untreated3_chr4.bam 
## "untreated3_chr4.bam.bai"

IGV可以在这个网站下载: https://www.broadinstitute.org/igv/home

下载的时候需要提供你的电子邮箱，现在我们使用IGV来查看一个基因，例如igs。

如果无法下载IGV，那么可以使用UCSC的 Genome Browser：http://genome.ucsc.edu/cgi-bin/hgTracks

UCSC Genome Browser是一个很好的资源，有许多涉及基因注释，多个物种的保护，以及ENCODE表观遗传信号。但是，UCSC Genome Browser要求你将基因组文件上传到他们的服务器，或将数据放在公共可用服务器上。如果你处理的数据涉及保密，那就另当别论了。

Simple plot

现在我们来看一使用plot()函数来绘制基因。

1
library(GenomicRanges)

注意：如果你使用的是Bioconductor 14，它与R 3.1配合使用，需要加载这个包。如想是Bioconductor 13（它与R 3.0.x配合使用），则不需要加载这个包。

1
2
library(GenomicAlignments)
## Loading required package: SummarizedExperiment

我们从文件fl1中读取比对结果(alignments)。然后我们使用coverage函数来计算碱基对的覆盖。然后我们提取与我们感兴趣的基因重叠的覆盖区子集，并将此覆盖率从RleList转换为numeric向量。我们前面提到，Rle对象是一种压缩对象，例如重复数字会被压缩为2个数字，其中一个是重复的数字，另外一个是表示这个数字重复的次数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
x <- readGAlignments(fl1)
xcov <- coverage(x)
z <- GRanges("chr4",IRanges(456500,466000))
# Bioconductor 2.14
xcov[z]
## RleList of length 1
## $$chr4
## integer-Rle of length 9501 with 1775 runs
##   Lengths: 1252   10   52    4    7    2 ...   10    7   12  392   75 1041
##   Values :    0    2    3    4    5    6 ...    3    2    1    0    1    0
# Bioconductor 2.13
xcov$chr4[ranges(z)]
## integer-Rle of length 9501 with 1775 runs
##   Lengths: 1252   10   52    4    7    2 ...   10    7   12  392   75 1041
##   Values :    0    2    3    4    5    6 ...    3    2    1    0    1    0
xnum <- as.numeric(xcov$chr4[ranges(z)])
plot(xnum)

plot of chunk unnamed-chunk-5

现在我们读取另外一个文件fl2，如下所示：

1
2
3
4
5
6
7
8
9
10
y <- readGAlignmentPairs(fl2)
## Warning in .make_GAlignmentPairs_from_GAlignments(gal, strandMode = strandMode, : 63 pairs (0 proper, 63 not proper) were dropped because the seqname
##   or strand of the alignments in the pair were not concordant.
##   Note that a GAlignmentPairs object can only hold concordant pairs at the
##   moment, that is, pairs where the 2 alignments are on the opposite strands
##   of the same reference sequence.
ycov <- coverage(y)
ynum <- as.numeric(ycov$chr4[ranges(z)])
plot(xnum, type="l", col="blue", lwd=2)
lines(ynum, col="red", lwd=2)

plot of chunk unnamed-chunk-6

我们可以单独放大一个外显子，如下所示：

1
2
plot(xnum, type="l", col="blue", lwd=2, xlim=c(6200,6600))
lines(ynum, col="red", lwd=2)

plot of chunk unnamed-chunk-7

使用转录本数据库提取感兴趣的基因

假设我们想可视化基因igs。我们可以从转录本数据库TxDb.Dmelanogaster.UCSC.dm3.ensGene 中提取它，但我们首先要知道这个基因的Ensembl基因名，现在我们使用前面学到的函数来实现这一功能，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
# biocLite("biomaRt")
library(biomaRt)
m <- useMart("ensembl", dataset = "dmelanogaster_gene_ensembl")
lf <- listFilters(m)
lf[grep("name", lf$description, ignore.case=TRUE),]
##                             name
## 1                chromosome_name
## 17         with_flybasename_gene
## 19   with_flybasename_transcript
## 23  with_flybasename_translation
## 52            external_gene_name
## 66              flybasename_gene
## 67        flybasename_transcript
## 68       flybasename_translation
## 93                 wikigene_name
## 107               go_parent_name
## 208               so_parent_name
##                                      description
## 1                                Chromosome name
## 17                   with FlyBaseName gene ID(s)
## 19             with FlyBaseName transcript ID(s)
## 23                with FlyBaseName protein ID(s)
## 52          Associated Gene Name(s) [e.g. BRCA2]
## 66           FlyBaseName Gene ID(s) [e.g. cul-2]
## 67  FlyBaseName Transcript ID(s) [e.g. cul-2-RB]
## 68     FlyBaseName Protein ID(s) [e.g. cul-2-PB]
## 93                 WikiGene Name(s) [e.g. Ir21a]
## 107                             Parent term name
## 208                             Parent term name
map <- getBM(mart = m,
  attributes = c("ensembl_gene_id", "flybasename_gene"),
  filters = "flybasename_gene", 
  values = "lgs")
map
##   ensembl_gene_id flybasename_gene
## 1     FBgn0039907              lgs

现在我们提取每个基因的外显子，然后提取igs基因的外显子，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
library(GenomicFeatures)
grl <- exonsBy(TxDb.Dmelanogaster.UCSC.dm3.ensGene, by="gene")
gene <- grl[[map$ensembl_gene_id[1]]]
gene
## GRanges object with 6 ranges and 2 metadata columns:
##       seqnames           ranges strand |   exon_id   exon_name
##          <Rle>        <IRanges>  <Rle> | <integer> <character>
##   [1]     chr4 [457583, 459544]      - |     63350        <NA>
##   [2]     chr4 [459601, 459791]      - |     63351        <NA>
##   [3]     chr4 [460074, 462077]      - |     63352        <NA>
##   [4]     chr4 [462806, 463015]      - |     63353        <NA>
##   [5]     chr4 [463490, 463780]      - |     63354        <NA>
##   [6]     chr4 [463839, 464533]      - |     63355        <NA>
##   -------
##   seqinfo: 15 sequences (1 circular) from dm3 genome

最后，我们绘制出这些外显子的范围，如下所示：

1
2
3
4
5
rg <- range(gene)
plot(c(start(rg), end(rg)), c(0,0), type="n", xlab=seqnames(gene)[1], ylab="")
arrows(start(gene),rep(0,length(gene)),
       end(gene),rep(0,length(gene)),
       lwd=3, length=.1)

plot of chunk unnamed-chunk-10

不过实际上这个基因位于负链，因此我们需要改一下箭头的方向，如下所示：

1
2
3
4
5
plot(c(start(rg), end(rg)), c(0,0), type="n", xlab=seqnames(gene)[1], ylab="")
arrows(start(gene),rep(0,length(gene)),
       end(gene),rep(0,length(gene)),
       lwd=3, length=.1, 
       code=ifelse(as.character(strand(gene)[1]) == "+", 2, 1))

plot of chunk unnamed-chunk-11

Gviz

我们将简要展示两个在Bioconductor中可视化基因组数据的软件包。请注意，这两个包中的任何一个都可以用于绘制多种数据的图形。我们将在这里展示如何像以前一样制作覆盖图(converage plot)：

1
2
3
4
5
6
#biocLite("Gviz")
library(Gviz)
## Loading required package: grid
gtrack <- GenomeAxisTrack()
atrack <- AnnotationTrack(gene, name = "Gene Model")
plotTracks(list(gtrack, atrack))

plot of chunk unnamed-chunk-12

提取覆盖。Gviz中的数据最好是是GRanges对象，因此我们需要将RleList转换为GRanges对象，如下所示：

1
2
3
4
5
xgr <- as(xcov, "GRanges")
ygr <- as(ycov, "GRanges")
dtrack1 <- DataTrack(xgr[xgr %over% z], name = "sample 1")
dtrack2 <- DataTrack(ygr[ygr %over% z], name = "sample 2")
plotTracks(list(gtrack, atrack, dtrack1, dtrack2))

plot of chunk unnamed-chunk-13

1
plotTracks(list(gtrack, atrack, dtrack1, dtrack2), type="polygon")

plot of chunk unnamed-chunk-13

ggbio

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#biocLite("ggbio")
library(ggbio)
## Loading required package: ggplot2
## Warning: replacing previous import by 'BiocGenerics::Position' when loading
## 'ggbio'
## Need specific help about ggbio? try mailing 
##  the maintainer or visit http://tengfei.github.com/ggbio/
## 
## Attaching package: 'ggbio'
## The following objects are masked from 'package:ggplot2':
## 
##     geom_bar, geom_rect, geom_segment, ggsave, stat_bin,
##     stat_identity, xlim
autoplot(gene)

plot of chunk unnamed-chunk-14

1
2
3
4
5
autoplot(fl1, which=z)
## reading in as Bamfile
## Parsing raw coverage...
## Read GAlignments from BamFile...
## extracting information...

plot of chunk unnamed-chunk-14

1
2
3
4
5
autoplot(fl2, which=z)
## reading in as Bamfile
## Parsing raw coverage...
## Read GAlignments from BamFile...
## extracting information...

plot of chunk unnamed-chunk-14

参考资料

  1. Visualizing NGS data
  2. IGV: https://www.broadinstitute.org/igv/home
  3. Gviz: http://www.bioconductor.org/packages/release/bioc/html/Gviz.html
  4. ggbio:http://www.bioconductor.org/packages/release/bioc/html/ggbio.html
  5. UCSC Genome Browser: zooms and scrolls over chromosomes, showing the work of annotators worldwide: http://genome.ucsc.edu/
  6. Ensembl genome browser: genome databases for vertebrates and other eukaryotic species:http://ensembl.org
  7. Roadmap Epigenome browser: public resource of human epigenomic data: http://www.epigenomebrowser.org http://genomebrowser.wustl.edu/ http://epigenomegateway.wustl.edu/
  8. “Sashimi plots” for RNA-Seq: http://genes.mit.edu/burgelab/miso/docs/sashimi.html
  9. Circos: designed for visualizing genomic data in a cirlce: http://circos.ca/
  10. SeqMonk: a tool to visualise and analyse high throughput mapped sequence data: http://www.bioinformatics.babraham.ac.uk/projects/seqmonk/

发表于 | 分类于
| 字数统计: 890 | 阅读时长 ≈ 4

使用ggbio包绘制染色体图

我们可以通过各种方法来绘制基因组级规模的数据。Bioiconductor中有两个这样的包，即Gviz 和 ggbio。这两个包旨在绘制基因组级数据结构的图形示意图。

ggbio 包中的 autoplot 方法用途广泛。对于GRanges实例来说，存在于每个范围中的数据都可以绘制成染色体上条带。核布图(karyogram)布局提供了全基因组的视图，不过这种图在处理染色体外序列水平方面有着重要作用。

原文如下：

The karyogram layout gives a genome-wide view, but it can be important to control the handling of extra-chromosomal sequence levels.

以下是肝细胞系的布局：

1
2
3
4
5
6
7
8
9
10
11
12
library(ERBS)
data(HepG2)
library(GenomeInfoDb)  # trim all but autosomal chroms
HepG2 = keepStandardChromosomes(HepG2)
data(GM12878)
GM12878 = keepStandardChromosomes(GM12878)
library(ggbio)
autoplot(HepG2, layout="karyogram", main="ESRRA binding on HepG2")
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

plot of chunk lkd

以下是B细胞系的布局：

1
2
3
4
5
autoplot(GM12878, layout="karyogram", main="ESRRA binding on GM12878")
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

使用Gviz包绘制ESRRA结合位点与基因注释

我们通常会将观察到的数据放在基因组的参考信息中来绘制这些数据的图形。现在我们演示一下如何ESRRA结合的图形。

第一步我们先加载相关的数据、注释包与Gviz包，如下所示：

1
2
3
4
5
library(ERBS)
library(Gviz)
library(Homo.sapiens)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene

ESRRA附近的基因

我们如何确定含有ESRRA以该基因附近的人类基因组片段呢？有几种方法，我们先从ENTREZ ID开始，如下所示：

1
2
3
4
5
6
library(Homo.sapiens)
eid = select(Homo.sapiens, keys="ESRRA", keytype="SYMBOL", columns="ENTREZID")
## 'select()' returned 1:1 mapping between keys and columns
eid
##   SYMBOL ENTREZID
## 1  ESRRA     2101

现在我们获取ESRRA基因的地址，并收拾它相邻基因的地址，并绑定这些基因的符号，如下所示：

1
2
3
4
5
6
allg = genes(txdb)
esrraAddr = genes(txdb, filter=list(gene_id=2101)) # redundant...
esrraNeigh = subsetByOverlaps(allg, esrraAddr+500000)
esrraNeigh$symbol = mapIds(Homo.sapiens, keys=esrraNeigh$gene_id, keytype="ENTREZID",
  column="SYMBOL")
## 'select()' returned 1:1 mapping between keys and columns

使用Gviz 包快速绘制图形，如下所示：

1
plotTracks(GeneRegionTrack(esrraNeigh, showId=TRUE))

plot of chunk lknei

此区域内的ESRRA结合峰

我们从某个样本（GM12878 EBV-transformed B-cell）获取ESRRA结合数据，以及我们目标基因附近的子集，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
data(GM12878)
sc = subsetByOverlaps(GM12878, range(esrraNeigh))
sc
## GRanges object with 6 ranges and 7 metadata columns:
##       seqnames               ranges strand |      name     score       col
##          <Rle>            <IRanges>  <Rle> | <numeric> <integer> <logical>
##   [1]    chr11 [64071338, 64073242]      * |         3         0      <NA>
##   [2]    chr11 [63913452, 63914077]      * |        61         0      <NA>
##   [3]    chr11 [63741506, 63742754]      * |       210         0      <NA>
##   [4]    chr11 [63992515, 63994725]      * |       408         0      <NA>
##   [5]    chr11 [64098744, 64100083]      * |      1778         0      <NA>
##   [6]    chr11 [63639263, 63640021]      * |      2248         0      <NA>
##       signalValue    pValue    qValue      peak
##         <numeric> <numeric> <numeric> <integer>
##   [1]       69.87       310        32      1040
##   [2]       52.48    123.72  1.785156       259
##   [3]       18.24    59.652  2.022276       622
##   [4]        9.04    41.001  1.910095       990
##   [5]        9.03    19.758  2.075721       653
##   [6]        8.65    17.481  2.031517       436
##   -------
##   seqinfo: 93 sequences (1 circular) from hg19 genome

计算表义文字(ideogram)从而给出染色体体的信息

计算过程有点慢。

1
idxTrack = IdeogramTrack(genome="hg19", chr="chr11")

数据合并

我们从顶部开始，用表意文字来确定染色体和染色体上的区域，并且通过观察到的数据和结构信息放大这个区域，如下所示：

1
2
3
4
plotTracks(list(idxTrack, GenomeAxisTrack(), 
   DataTrack(sc[,7], name="ESRRA peak values"), 
   GeneRegionTrack(esrraNeigh, showId=TRUE,
         name="genes near ESRRA"), GenomeAxisTrack()))

plot of chunk dofull

参考资料

  1. Sketching the binding landscape over chromosomes with ggbio’s karyogram layout
  2. Gviz for plotting data with genomic features

发表于 | 分类于
| 字数统计: 4,240 | 阅读时长 ≈ 18

前言

在这一节中我们将学习基因集分析(gene set analysis)背后的统计学原理。在这一部分里，我们要记住三个重要的步骤，这三个步骤都非常重要，但它们又相互独立，分别如下所示：

  • 参与统计的基因集
  • 汇总每个基因集的统计学指标
  • 评估不确定性的方法

准备数据与基因集数据库

1
2
source("http://www.bioconductor.org/biocLite.R")
biocLite("GSEABase")

加载以下包：

1
2
3
4
5
6
library(rafalib)
library(GSEABase)
library(GSE5859Subset)
library(sva)
library(limma)
library(matrixStats)

我们使用下面包中的数据来说明一下，这个包可以从Github中下载，如下所示：

1
2
library(devtools)
install_github("genomicsclass/GSE5859Subset")

我们构建的这个数据含有批次效应，因此我们需要使用SVA对其进行校正，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
library(GSE5859Subset)
data(GSE5859Subset)
library(sva)
library(limma)
library(matrixStats)
library(rafalib)
X = sampleInfo$group
mod<-model.matrix(~X)
svafit <- sva(geneExpression,mod)
## Number of significant surrogate variables is:  5 
## Iteration (out of 5 ):1  2  3  4  5
svaX<-model.matrix(~X+svafit$sv)
lmfit <- lmFit(geneExpression,svaX)
tt<-lmfit$coef[,2]*sqrt(lmfit$df.residual)/(2*lmfit$sigma)
pval<-2*(1-pt(abs(tt),lmfit$df.residual[1]))
qval <- p.adjust(pval,"BH")

我们需要注意一下，在这一节中，我们需要每个基因的评分和差异表达基因基因的候选列表。我们如何获取这个候选基因列表与所表示的基因丰度的单位无关。为了说明这一点，我们会假装我们事先不知道在不同组之间基因表达的关系。

我们构建一个FDR为25%的候选基因列表。这个列表中有23个基因。我们从中能了解到什么生物学信息呢？既然我们知道基因存在着相互作用，那么研究一组基因，而不是单个基因有着重要意义。通常来说，这就是基因集分析的原始想法。

网上有各种提供基因集的资源。这些资源中的基因都是经过缜密筛选的，它们都有着一些共同的特征。这些基因分组的方式是通过通路（例如KEGG）或GO来划分的。在下面的网站上就能找到许多这样的资源：

http://www.broadinstitute.org/gsea/msigdb/index.jsp

在以下的案例中，我闪将使用一个根据染色体进行分组的基因集来分析数据。第一步就是下载数据库，不过你需要在网站注册。

当你下载好文件后，就可以通过以下的命令读入数据库（文件扩展名为gmt，其实就是GSEA分析中的基因集文件），如下所示：

1
gsets <- getGmt("c1.all.v6.1.entrez.gmt") ##must provide correct path to file

这个对象中含有326个基因集。每个基因集使用ENTREZ ID标明了基因，如下所示：

1
2
3
4
5
6
7
8
9
10
length(gsets)
## [1] 326
head(names(gsets))
## [1] "chr5q23"  "chr16q24" "chr8q24"  "chr13q11" "chr7p21"  "chr10q23"
gsets[["chryq11"]]
## setName: chryq11 
## geneIds: 728512, 392597, ..., 360027 (total: 204)
## geneIdType: Null
## collectionType: Null 
## details: use 'details(object)'

我们可以使用以下方式获取ENTREZ ID，如下所示：

1
2
head(geneIds(gsets[["chryq11"]]))
## [1] "728512" "392597" "9086"   "359796" "83864"  "425057"

现在，为了将其应用于我们的数据，我们必须将这些ID映射到Affymetrix ID上。现在我们编写这样一个函数，这个函数会使用Expression Set类，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
mapGMT2Affy <- function(object,gsets){
  ann<-annotation(object)
  dbname<-paste(ann,"db",sep=".")
  require(dbname,character.only=TRUE)
  gns<-featureNames(object)
  ##This call may generate warnings
  map<-select(get(dbname), keys=gns,columns=c("ENTREZID", "PROBEID"))
  map<-split(map[,1],map[,2])
  indexes<-sapply(gsets,function(ids){
    gns2<-unlist(map[geneIds(ids)])
    match(gns2,gns)
    })
  names(indexes)<-names(gsets)
  return(indexes)
}
##create an Expression Set
rownames(sampleInfo)<- colnames(geneExpression)
e=ExpressionSet(assay=geneExpression,
                phenoData=AnnotatedDataFrame(sampleInfo),
                annotation="hgfocus")
##can safely ignore the warning
gsids <- mapGMT2Affy(e,gsets) 
## 'select()' returned 1:many mapping between keys and columns

基于相关检验的方法

现在我们就可以计算差异表达基因是否在某个基因集中富集。那么我们如何计算呢？最简单的方法就是独立卡方检验，这种统计方法我们以前学过。现在我们对这些差异基因是否在Y染色体上的基因集中富集进行计算，如下所示：

1
2
3
4
5
6
7
8
9
tab <- table(ingenset=1:nrow(e) %in% gsids[["chryq11"]],signif=qval<0.05)
tab
##         signif
## ingenset FALSE TRUE
##    FALSE  8769    9
##    TRUE     11    4
chisq.test(tab)$p.val
## Warning in chisq.test(tab): Chi-squared approximation may be incorrect
## [1] 5.264033e-121

因为我们比较的是男性与女性，因此我们使用Y染色体上基因集来进行富集，这样比较便于富集差异基因，更能说明这种算法的原理。但是，由于只有13个基因有显著性（在Y染色体上的只有4个），这种方法对于寻找可能在某些方面具有生物学意义的基因集来说，实际意义不大。例如，X染色体上的基因集是我们感兴趣的基因集么？这里需要注意的是，一些基因在X染色体上并不活性，因此与男性相比，女性体内这样的基因数目更多。

这里我们面临的问题是，FDR值为0.25，这是任意设置的一个截止值。为什么不设为0.05或0.01呢？我们如何设置这个截止值？

如下所示：

1
2
3
4
5
library(rafalib)
mypar(1,1)
qs<-seq(0,1,len=length(tt)+1)-1/(2*length(tt));qs<-qs[-1]
qqplot(qt(qs,lmfit$df.resid),tt,ylim=c(-10,10),xlab="Theoretical quantiles",ylab="Observed") ##note we leave some of the obvious ones out
abline(0,1)

plot of chunk unnamed-chunk-10

基因集相关的统计学

目前大多数进行基因集分析的方法都没有将基因分析为差异表达部分和非差异表达部分。相反，我们会使用一些类似于t检验的方法来进行计算。为了了解这种方法比基因列表的二分法更为强大，现在我们提出一个问题：一些基因是否在X染色体上并没有失活？我们不像前面提到的那样创建一个基因表，而是查看一下X染色体中第一组基因t检验中t值的分布，如下所示：

1
2
3
4
5
6
7
ind <- gsids[["chrxp11"]]
mypar(1,1)
plot(density(tt[-ind]),xlim=c(-7,7),main="",xlab="t-stat",sub="",lwd=4)
lines(density(tt[ind],bw=.7),col=2,lty=2,lwd=4)
rug(tt[ind],col=2)
legend(-6.5, .3, legend=c("present on xp11", "not on xp11"),
     col=c(2,1), lty=c(2,1), lwd=4)

plot of chunk unnamed-chunk-11

从上图我们可以看到基因集的t统计量的分布，它的分布似乎不同于所有基因的分布。那么我们如何地这个结果进行定量呢？

有关基因集分析的第一篇文献是于2001年由[Virtaneva]发表的，他在前面部分中讨论了Wilcoxon Mann-Whitney检验。那篇文献的基本思路就是统计实验组和对照组之间的差异基因在某个特定的基因集中的数目是否比在其它基因集中的数目更大。Wilcoxon比较的是基因集上基因的秩是否在中位数以上。

在这里我们计算了每个基因集的Wilcoxon检验。然后我们对Wilcoxon统计量的结果进行标准(standardize)，也就是将它们转换为均值为0，标准差为1的数据。注意，中心极限理论告诉我们，如果基因集足够大，并且每个基因的效应大小彼此独立，那么这些统计量将遵循正态分布，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
es <- lmfit$coef[,2]
wilcox <- t(sapply(gsids,function(i){
  if(length(i)>2){
  tmp<- wilcox.test(es[i],es[-i])  
  n1<-length(i);n2<-length(es)-n1
  z <- (tmp$stat -n1*n2/2) / sqrt(n1*n2*(n1+n2+1)/12)
  return(c(z,tmp$p.value))
  } else return(rep(NA,2))  
  }))
mypar(1,1)
cols <- rep(1,nrow(wilcox))
cols[grep("chrx",rownames(wilcox))]<-2
cols[grep("chry",rownames(wilcox))]<-3
qqnorm(wilcox[,1],col=cols,cex=ifelse(cols==1,1,2),pch=16)
qqline(wilcox[,1])
legend("topleft",c("Autosome","chrX","chrY"),pch=16,col=1:3,box.lwd=0)

plot of chunk unnamed-chunk-12

从上图我们可以看出来，排名前10位的基因并没像我们预期的那样在X染色体和Y染色体中高表达。使用Wilcoxon（稳健,robust)为基础的背后逻辑是，我们并不想看到一个或两个高表达的基因对基因集分析产生重要影响（Wilcoxon检验是功效比较低，但比较稳健，对异常值不敏感）。虽然这种统计方法能够降低低秩基因的假阳性，但是这牺牲了灵敏度（功效,power）。下面我们来考虑一个低稳健性(less robust)，但是通常情况下更加强大的方法。

如果我们对基个基因集中向左分布（通常来说是表达量少的基因）的基因或者说是向右分布（通常来说是就是表达高的基因）的基因戌兴趣，那么我们可以使用的最简单的统计就是比较两组之间的均值，具体可以参考这篇文献《Gene set enrichment analysis made simple》，其中G是基因在基因集中的索引，我们定义以下的统计量：

零假设是基因集中的基因表达没有差异，这些t检验的均值为0。如果t检验统计量之间互相是独立的（后面我们会学习避免此假设的方法），那么$\sqrt{N} \bar{t}$ 服从以下分布：

以下是这些统计量的qq图：

1
2
3
4
avgt <- sapply(gsids,function(i) sqrt(length(i))*mean(tt[i]))
qqnorm(avgt,col=cols,cex=ifelse(cols==1,1,2),pch=16)
qqline(avgt)
legend("topleft",c("Autosome","chrX","chrY"),pch=16,col=1:3,box.lwd=0)

plot of chunk unnamed-chunk-13

1
2
3
4
5
avgt[order(-abs(avgt))[1:10]]
##    chryp11    chryq11    chrxq13    chr6p22   chr17p13    chrxp11 
## -17.749865 -15.287203   6.592806  -4.558008  -4.545851   4.420686 
##   chr20p13    chr4q35    chr7p22    chrxp22 
##  -4.241892  -4.033233  -3.642068   3.608926

有一些算法就是基于以上原理来转统计量的均值的，例如 Tian et al 2005 , SAFE, ROAST, CAMERA.

我们可以构建其它的得分来检验其它替代假设。例如我们可以检验changes in variances 与这个分布的常规变化 [1], [2]。

基因集的假设检验

上面我们对$\sqrt{N} \bar{t}$ 的正态近似表明，t检验的统计量是独立的，并且分布相同，例如，现在我们看第二个公式，如下所示：

这个公众就是假设t检验统计量是独立的。但是，即使在零假设下，我们也许会考虑到基因集之间是相互关联的。为了了解这一点，我们可以取每个基因集中所有成对相关的平均值，我们就会看到，它们并不都是0，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
N <- sapply(gsids,length)
ind1<-which(X==0)
ind2<-which(X==1)
corrs <- t(sapply(gsids,function(ind){
  if(length(ind)>=2){
    cc1<-cor(t(geneExpression[ind,ind1]))
    cc2<-cor(t(geneExpression[ind,ind2]))
  return(c(median(cc1[lower.tri(cc1)]),
    median(cc2[lower.tri(cc2)])))
  } else return(c(NA,NA))
}))
mypar(1,1)
plot(corrs[N>10,],xlim=c(-1,1),ylim=c(-1,1),xlab="Correlation within females",ylab="Correlation within males")
abline(h=0,v=0,lty=2)

plot of chunk unnamed-chunk-14

注意，在这个案例中，在一个基因集中，基因的检验统计量具有相关性的 $\rho$， 平均t检验统计量的方差如下所示：

在上面公式里，$\rho$ 是正值，方差与 $\{1 + (N-1) \rho \} $有关。虽然常规情况下，在一个基因集内，$\rho$ 并不一定都相等，但是我们仍然可以根据一个基因集内的平均 $\rho$ 来计算一个校正因子(correction factor)。这个公式可以参考这篇文献《A statistical framework for testing functional categories in microarray data》。公式中使用的校正方法还包括Wald统计量。

这些校正因子可以对上述提到的统计量进行简单的汇总。事实上，现在已经有论文提出了基于简单的均值转换方法的算法，这些算法还考虑到了相关性，它们基于渐进近似的ROAST和CAMERA，用于开发统计摘要(summary statistics)与统计推断。在计算机演示案例中，我们将会演示其中一款软件。简单起见，我们这里基于 $\sqrt{N} \bar{t}$ 计算近似，并使用校正因子。这里我们比较一下初始 $\sqrt{N} \bar{t}$ 以及校正后的 $\sqrt{N} \bar{t}$ ，如下所示：

在下图里，我们注意到这个校正会降低分数(socres)（也就是将分数值降低到接近0的程度）。尤其注意的是，第三个最高的值现在非常接近0，它不再显著（参见箭头），如下所示：

1
2
3
4
5
6
7
8
9
10
11
avgcorrs <- rowMeans(corrs)
cf <- (1+(N-1)*avgcorrs)
cf[cf<1] <- 1 ## we ignore negative correlations
correctedavgt <- avgt/sqrt(cf)
parampvaliid <- 2*pnorm(-abs(avgt))
parampval<- 2*pnorm(-abs(correctedavgt))
plot(avgt,correctedavgt,bg=cols,pch=21,xlab="Original",ylab="With correction factor",xlim=c(-7,20),ylim=c(-7,20),cex=1.5)
abline(0,1)
abline(h=0,v=0,lty=2)
thirdhighest <- order(-avgt)[3]
arrows(avgt[thirdhighest]+3,correctedavgt[thirdhighest],x1=avgt[thirdhighest]+0.5,lwd=2)

plot of chunk unnamed-chunk-15

需要注意，标注出来的基因集具有相对高的相关性和大小(size)，它就使得校正因子变得很大，如下所示：

1
2
3
4
5
6
7
8
avgcorrs[thirdhighest]
##   chrxp22 
## 0.0287052
length(gsids[[thirdhighest]])
## [1] 76
cf[thirdhighest]
## chrxp22 
## 3.15289

在这一节中，我们已经使用了平均t检验统计量作为一个简单的案例来说明基因集在零假设下的统计摘要。我们证明了，当基因集中的基因并不独立时，需要一个校正因子。这里要注意的是，我们使用平均t检验统计量只是做了一个简单的说明，如果需要更加严格处理方法来实现均值偏移方法(testing-for-mean-shifts)，则需要参考 SAFE, ROAST 和 CAMERA。

即使在使用了膨胀因子(inflation factor)校正了平均t检验统计量后，零分布也有可能不成立。例如，对于平均t检验统计量是正态的假设而言，原始的样本大小和基因集大小有可能太小。在一次部分里，我们会演示如何在不使用参数假设的情况下，使用置换检验(Permutations)来估计零分布。

置换检验(Permutations)

当我们感兴趣的结果允许以及样本足够的情况下，我们可以使用基因集特定的转换检验来替换参数检验。参考文献中列出了几种这样的方法[1], [2]。对于大多数的富集分布而言，我们可以使用置换来创建一个零分布。例如，在这个案例中，我们对简单的t检验摘要执行置换检验，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
library(matrixStats)
## matrixStats v0.50.1 (2015-12-14) successfully loaded. See ?matrixStats for help.
## 
## Attaching package: 'matrixStats'
## The following object is masked from 'package:GenomicFiles':
## 
##     rowRanges
## The following object is masked from 'package:SummarizedExperiment':
## 
##     rowRanges
## The following objects are masked from 'package:genefilter':
## 
##     rowSds, rowVars
## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians
set.seed(1)
B <- 400 ##takes a few minutes
null <- sapply(1:B,function(b){
 nullX<- sample(X)
 nullsvaX<-model.matrix(~nullX+svafit$sv) ##note that we are not recomupting the surrogate values. 
 nulllmfit <- lmFit(geneExpression,nullsvaX)
 nulltt<-nulllmfit$coef[,2]*sqrt(nulllmfit$df.residual)/(2*nulllmfit$sigma)
 nullavgt <- sapply(gsids,function(i) sqrt(length(i))*mean(nulltt[i]))
 return(nullavgt)
})
permavgt <- avgt/rowSds(null)
permpval<- rowMeans(abs(avgt) < abs(null))

在这个案例中，使用置换得到的结果非常类似于使用了校正因子得到的参数检验结果，我们从比较的结果中就能看出这些信息，如下所示：

1
2
3
plot(correctedavgt,permavgt,bg=cols,pch=21,xlab="Parametric z-score (with correction)",ylab="Permutation z-score",cex=1.5,ylim=c(-5,15),xlim=c(-5,15))
abline(0,1)
abline(h=0,v=0,lty=2)

plot of chunk unnamed-chunk-18

我们看到这两种方法得到的q值也是比较类似：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
tab <- data.frame(avgt=p.adjust(signif(2*pnorm(1-abs(avgt)),2),method="BH"),
                  correctedavgt=p.adjust(signif(2*pnorm(1-abs(correctedavgt)),2),method="BH"),
                  permutations=p.adjust(permpval,method="BH"))
##include only gene sets with 10 or more genes in comparison
tab<-tab[N>=10,]
tab <- cbind(signif(tab,2),apply(tab,2,rank))
tab<-tab[order(tab[,1]),]
tab <- tab[tab[,1]< 0.25,]
tab
##             avgt correctedavgt permutations avgt correctedavgt
## chryp11  1.9e-60       1.5e-37         0.00  1.0             1
## chryq11  4.2e-44       6.3e-02         0.00  2.0             3
## chrxq13  2.4e-06       5.3e-03         0.00  3.0             2
## chr17p13 2.5e-02       1.0e+00         0.00  4.5           124
## chr6p22  2.5e-02       1.0e+00         0.77  4.5           124
## chrxp11  3.4e-02       1.0e+00         0.70  6.0           124
## chr20p13 5.6e-02       9.7e-02         0.00  7.0             4
## chr4q35  9.7e-02       1.0e+00         0.39  8.0           124
##          permutations
## chryp11             4
## chryq11             4
## chrxq13             4
## chr17p13            4
## chr6p22           114
## chrxp11            78
## chr20p13            4
## chr4q35            18

不过置换检验方法的一个局限就是，为了能够得到非常小的p值，我们需要进行多次置换。例如为了得到一个在10e-5到10e-6之间的一个p值，我们霜肆进行100万次置换。对于那些使用了保守的多重比较校正方法（例如Bonferonni校正）的人来说，很有必要进行更多次地的模拟。置换检验的另外一个局限就是，在更复杂的设计中，例如具有许多轭和连续协变量的实验中，如何对置换进行解读并不容易。

参考资料

  1. Gene set analysis
  2. Virtaneva K, “Expression profiling reveals fundamental biological differences in acute myeloid leukemia with isolated trisomy 8 and normal cytogenetics.” PNAS 2001 http://www.ncbi.nlm.nih.gov/pubmed/?term=11158605
  3. Jelle J. Goeman and Peter Bühlmann, “Analyzing gene expression data in terms of gene sets: methodological issues” Bioinformatics 2006. http://bioinformatics.oxfordjournals.org/content/23/8/980

发表于 | 分类于
| 字数统计: 1,713 | 阅读时长 ≈ 8

简单t检验

在这一节中，我们将展示一下使用简单t检验和使用limma层次模型进行差异分析的区别。limma的相关参考文献已经列到参考资料中。

在这一节中，我们还民法了执行limma分析的基本步骤。需要注意的是，limma的功能非常强大，它有数百页的文档，我们不能详细地介绍这个包以及函数，有兴趣的同学可以直接查看它的文档。

spike-in数据集

我们先加载一批由Affymetrix生成的均一化(normalized)后的spike-in数据，如下所示：

1
2
3
4
# biocLite("SpikeInSubset")
library(SpikeInSubset)
data(rma95)
fac <- factor(rep(1:2,each=3))

简单来说，这里一共有16个mRNA的物种，它们都已经使用了固定的浓度，并混合起来使用了hgu95芯片进行了检验。这个子集中对于每个mRNA物种，ExpressionSet中的第一个三元组(trio)和第二个三元组有着固定的不同值。通过pData就能看到，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
pData(rma95)
##                  37777_at 684_at 1597_at 38734_at 39058_at 36311_at
## 1521a99hpp_av06      0.00   0.25     0.5        1        2        4
## 1532a99hpp_av04      0.00   0.25     0.5        1        2        4
## 2353a99hpp_av08      0.00   0.25     0.5        1        2        4
## 1521b99hpp_av06      0.25   0.50     1.0        2        4        8
## 1532b99hpp_av04      0.25   0.50     1.0        2        4        8
## 2353b99hpp_av08r     0.25   0.50     1.0        2        4        8
##                  36889_at 1024_at 36202_at 36085_at 40322_at 407_at
## 1521a99hpp_av06         8      16       32       64      128   0.00
## 1532a99hpp_av04         8      16       32       64      128   0.00
## 2353a99hpp_av08         8      16       32       64      128   0.00
## 1521b99hpp_av06        16      32       64      128      256   0.25
## 1532b99hpp_av04        16      32       64      128      256   0.25
## 2353b99hpp_av08r       16      32       64      128      256   0.25
##                  1091_at 1708_at 33818_at 546_at
## 1521a99hpp_av06      512    1024      256     32
## 1532a99hpp_av04      512    1024      256     32
## 2353a99hpp_av08      512    1024      256     32
## 1521b99hpp_av06     1024       0      512     64
## 1532b99hpp_av04     1024       0      512     64
## 2353b99hpp_av08r    1024       0      512     64

如果要了解这批数据这样设计的原理，现在我们来看一下下面的四张spiked mRNA，原文如下：

To get a feel for the response of the array quantifications to this design, consider the following plots for four of the spiked mRNAs.

1
2
3
4
5
6
7
8
par(mfrow=c(2,2))
for (i in 1:4) {
  spg = names(pData(rma95))
  plot(1:6, exprs(rma95)[spg[i+6],], main=spg[i+6], ylab="RMA",
    xlab="nominal", axes=FALSE)
  axis(2)
  axis(1, at=1:6, labels=pData(rma95)[[spg[i+6]]])
}

plot of chunk lkpl

由于RMA(robust multi-array average,鲁棒多阵列平均 )是以log2转换为定量标准的，因此在归一化的检测值中观察到的差异是比较合理的。但是，在每个设计好的浓度内部，存在着相当大的变异。

rowttests

我们现在使用 genefilter 包中的 rowttests 函数来做一个简单的t检验，如下所示： 

1
2
library(genefilter)
rtt <- rowttests(exprs(rma95),fac)

我们将根据p值，均值，以及特征值是否是spike-in值来设置颜色，如下所示：

1
2
3
mask <- with(rtt, abs(dm) < .2 & p.value < .01)
spike <- rownames(rma95) %in% colnames(pData(rma95))
cols <- ifelse(mask,"red",ifelse(spike,"dodgerblue","black"))

我们现在使用上面定义的颜色来绘制结果。我们将dm乘以-1，因为国我们感兴趣的是，第二组与第一组的差异（这是由lm与limma包默认情况下使用的差异）。spike-in基因是蓝色的，它们大多数有着较小的p值与较大的均值差异。红点表示较小p值的基因，同时也有着较小的均值差异。我们将会看到这些点使用limma后的变化，如下所示：

1
2
3
4
5
with(rtt, plot(-dm, -log10(p.value), cex=.8, pch=16,
     xlim=c(-1,1), ylim=c(0,5),
     xlab="difference in means",
     col=cols))
abline(h=2,v=c(-.2,.2), lty=2)

Note that the red genes have mostly low estimates of standard deviation.

1
2
3
4
rtt$s <- apply(exprs(rma95), 1, function(row) sqrt(.5 * (var(row[1:3]) + var(row[4:6]))))
with(rtt, plot(s, -log10(p.value), cex=.8, pch=16,
              log="x",xlab="estimate of standard deviation",
              col=cols))

plot of chunk unnamed-chunk-5

limma步骤

以下是执行基本limma分析的三个步骤。我们指定coef=2，因为我感兴趣的是组与组之间的差异，而不是截距(intercept)，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
library(limma)
options(digits=3)
fit <- lmFit(rma95, design=model.matrix(~ fac))  # step 1 least squares estimates
colnames(coef(fit))     
## [1] "(Intercept)" "fac2"
fit <- eBayes(fit)           # step 2 moderate the t statistics
tt <- topTable(fit, coef=2)  # step 3 report
tt
##           logFC AveExpr      t  P.Value adj.P.Val    B
## 1708_at  -7.061    7.95 -73.53 7.82e-17  9.87e-13 8.65
## 36202_at  0.853    9.37   9.98 4.94e-07  3.12e-03 4.59
## 36311_at  0.832    8.56   8.36 3.02e-06  1.27e-02 3.57
## 33264_at  0.712    4.92   7.43 9.67e-06  2.71e-02 2.84
## 32660_at  0.655    8.68   7.36 1.07e-05  2.71e-02 2.77
## 38734_at  0.747    6.26   7.19 1.35e-05  2.83e-02 2.62
## 1024_at   0.843    9.70   6.73 2.50e-05  4.40e-02 2.20
## 36085_at  0.645   12.19   6.65 2.79e-05  4.40e-02 2.12
## 33818_at  0.532   12.29   6.45 3.70e-05  5.19e-02 1.92
## 39058_at  0.609    7.53   6.28 4.77e-05  5.69e-02 1.74

topTable 会返回你定义的值排序的前几个基因。你还可以使用topTable通过指定none来要求返回所有的值。需要注意的是，有一列会给出校正后的每个基因的p值。默认情况下，通过p.adjust函数中的Benjamini-Hochberg方法用于校正p值，如下所示：

1
2
3
4
# ?topTable
dim(topTable(fit, coef=2, number=Inf, sort.by="none"))
## [1] 12626     6
# ?p.adjust

在这里，我们比较一下使用limma分析的结果与以前得到的火山图结果。注意，红点现在都位于-log10(p.value)=2之下。此外，表示实际差异的蓝点的p值比以前更高，如下所示：

1
2
3
4
5
limmares <- data.frame(dm=coef(fit)[,"fac2"], p.value=fit$p.value[,"fac2"])
with(limmares, plot(dm, -log10(p.value),cex=.8, pch=16,
     col=cols,xlab="difference in means",
     xlim=c(-1,1), ylim=c(0,5)))
abline(h=2,v=c(-.2,.2), lty=2)

最后我们来构建一张图，用于展示limma是如何将方差估计缩小到一个共同值(common value)的，这能消除由于方差估计过低而导致的假阳性。

在这里，我们为不同的方差估计值设计40个小区间(bin)，每个区间中落入一个基因。我们会删除那些不含有任何基因的区间，原文如下：

Here we pick, for each of 40 bins of different variance estimates, a single gene which falls in that bin. We remove bins which do not have any such genes.

1
2
3
4
n <- 40
qs <- seq(from=0,to=.2,length=n)
idx <- sapply(seq_len(n),function(i) which(as.integer(cut(rtt$s^2,qs)) == i)[1])
idx <- idx[!is.na(idx)]

现在绘制一条直线，这个直线展示了那些基因最初的估计方差到运行了limma之后的估计方差，如下所示：

1
2
3
4
5
6
par(mar=c(5,5,2,2))
plot(1,1,xlim=c(0,.21),ylim=c(0,1),type="n",
     xlab="variance estimates",ylab="",yaxt="n")
axis(2,at=c(.1,.9),c("before","after"),las=2)
segments((rtt$s^2)[idx],rep(.1,n),
         fit$s2.post[idx],rep(.9,n))

plot of chunk unnamed-chunk-10

关于limma校正的一些层级模型资料，可以参考以前的笔记《DALS019-统计模型2-贝叶斯分布与层次分析》。

参考资料

  1. Using limma for microarray analysis
  2. Smyth GK, “Linear models and empirical bayes methods for assessing differential expression in microarray experiments”. Stat Appl Genet Mol Biol. 2004 http://www.ncbi.nlm.nih.gov/pubmed/16646809

发表于 | 分类于
| 字数统计: 1,000 | 阅读时长 ≈ 4

前言

在上一节中，我们集中研究了一对基因来说明两个方面的变异。其中一个基因似乎在不同的小鼠之间存在着很高的变异，但是这种变异隐藏在了合并数据背后。当我们使用合并的数据对不同的品系小鼠进行统计时，获得的p值在表面上看来是非常低的（p<10e-6），但是当我们使用含有几个代表个体的数据进行比较时，我们发现p值又变得比较大（p=0.089）。在这一节中我们会认识到，找到最令人信服的问题就是要研究生物学变异，只有通过良好的实验设计才能直接检测到生物学变异。准确解释生物学变异的来源(origin)和大小(size)需要合适的统计分析。

在这一节上，我们将会介绍基因组实验中的统计推断，在此，我们需要了解以下几个方面的注意事项：

  • 如果需要进行多次检验，那么对于数以千计的特征值（通常是一个基因）来说，每个特征值至少要经历一次检验；
  • 每个特征值都会表现出自身的技术变异与生物学变异；
  • 有可能存在未检测到或未报道的生物学变异（例如以天为单位时间时，不同时间点的检测）（注：原文是there may be unmeasured or unreported sources of biological variation (such as time of day)）；
  • 许多特征值有着内在的联系，这些检验并不是独立的。

我们将会使用到前面章节中提到的一些概念，包括t检验和多重比较；在后面部分里，我们将会根据分层模型(hierarchical models)计算标准偏差估计。

我们先加载合并的数据，如下所示：

1
2
3
4
5
library(Biobase)
library(maPooling)
data(maPooling)
pd=pData(maPooling)
individuals=which(rowSums(pd)==1)

现在我们提取单个小鼠及品系，如下所示：

1
2
3
4
individuals=which(rowSums(pd)==1)
individuals=individuals[-grep("tr",names(individuals))]
y=exprs(maPooling)[,individuals]
g=factor(as.numeric(grepl("b",names(individuals))))

T检验

我们可以使用genefilter包中的rowttest函数来进行t检验，如下所示：

1
2
library(genefilter)
tt=rowttests(y,g)

现在我们来找一些哪些基因有统计学意义。由于某些历史原因，在科学研究中，我们判断是否有统计学意义的截止值是p值为0.01或0.05，现在我们找到那些p值小于0.01的基因数目，如下所示：

1
2
3
4
5
6
NsigAt01 = sum(tt$p.value<0.01)
NsigAt01
## [1] 1578
NsigAt05 = sum(tt$p.value<0.05)
NsigAt05
## [1] 2908

多重检验

我们在前面的章节中提到了多重检验（参考笔记DALS017-高维数据推断2-统计原理），现在我们来进行一个快速总结。

我们是否要报道前面提到的名义上的所有的具有显著性的基因？现在让我们来研究一下，如果将第一组分为两组，强制零假设为真，会发生什么情况。

1
2
3
4
5
6
7
8
9
set.seed(0)
shuffledIndex <- factor(sample(c(0,1),sum(g==0),replace=TRUE ))
nulltt <- rowttests(y[,g==0],shuffledIndex)
NfalselySigAt01 = sum(nulltt$p.value<0.01)
NfalselySigAt01 
## [1] 79
NfalselySigAt05 = sum(nulltt$p.value<0.05)
NfalselySigAt05
## [1] 840

如果我们使用0.05作为截止值，则将出现840个假阳性。我们在前面部分中已经提到了几种校正方法，其中就包括qvalue包中的qvalue函数。经过这样的调整，我们就能降低假阳性，如下所示：

1
2
3
4
5
6
library(qvalue)
qvals = qvalue(tt$p.value)$qvalue
sum(qvals<0.05)
## [1] 1179
sum(qvals<0.01)
## [1] 538

并且在零亿是高的情况下，也不会产生误报，如下所示：

1
2
3
4
5
6
library(qvalue)
nullqvals = qvalue(nulltt$p.value)$qvalue
sum(nullqvals<0.05)
## [1] 0
sum(nullqvals<0.01)
## [1] 0

这以相当常规的方式解决了在固定名义显著性水平上执行多次假设检验时造成的假阳性问题。

参考资料

  1. Inference: t-tests, multiple comparisons

发表于 | 分类于
| 字数统计: 2,653 | 阅读时长 ≈ 11

前言

在后面的章节中，我们将介绍基因组学实验中的统计推断。我们会使用前面几节中介绍的一些概念，其中包括t检验和多重比较。在这一节中，我们先介绍一个在基因组学数据分析中一个极为重要的概念：生物学变异与技术变异(biological and technical variability)以及它们之间的区别。

通常而言，我们在一个种群里所观察到的生物学单位，例如某个个体的变异称为生物学变异(biological variability)。我们将那些对同一个生物学单位进行测量时所观察到的变异称为技术性变异(technical variability，例如同一个样本测3次，这3次就是技术重复，3次数据的变异称为技术变异)。由于在基因组学研究中经常会用到新技术，因此技术重复(technical replicates)在很多时候都会用于评估实验数据。理想情况下，同一样本的检测结果应该是相同的，但在实际检测中同一样本的多次检测会出现一定的偏差，因此使用技术重复我们就能我们就能够评估新技术的技术变异。我还经常使用生物学重复(biological replicates)与技术重复(technical replicates)来分别代指生物学变异和技术变异。

我们在进行统计推断来解释有差异时，需要注意不要混淆生物学变异和技术变异，因为它们所代表的意义是不同的。例如，当我们分析技术重复时，样本仅仅只有一个，而不是代表某个总体，例如人类或小鼠。这一节中，我们通过设计一个实验来说明一下技术重复和生物重复。

合并实验数据

我们将要研究的数据集来源于基因表达微阵列。在这个实验中，我们从两个品系的小鼠中分别随机选择了12只小鼠，并提取了它们的RNA。我们将得到的24个样本使用微阵列的手段进行检测，引外，我们还构成了两个库，这两个库中的样本则是分别来源于这两个品系的12只小鼠的混合样本。

现在我们安装以下包，这个包中含有这些数据，如下所示：

1
2
library(devtools)
install_github("genomicsclass/maPooling")

使用pData函数我们就能查看实验设计。每行表示一个样本，每列表示1只小鼠。例如i,j单元格，中如果是1，则表示RNA数据来源于j只小鼠的样本i。通过行名我们可以知道品系信息（不推荐这种方法来进行分组，你可以向phenoData中添加一个变量，用于指定小鼠品系信息），数据如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
library(Biobase)
library(maPooling)
data(maPooling)
head(pData(maPooling))
##          a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a14 b2 b3 b5 b6 b8 b9 b10 b11
## a10       0  0  0  0  0  0  0  0   1   0   0   0  0  0  0  0  0  0   0   0
## a10a11    0  0  0  0  0  0  0  0   1   1   0   0  0  0  0  0  0  0   0   0
## a10a11a4  0  0  1  0  0  0  0  0   1   1   0   0  0  0  0  0  0  0   0   0
## a11       0  0  0  0  0  0  0  0   0   1   0   0  0  0  0  0  0  0   0   0
## a12       0  0  0  0  0  0  0  0   0   0   1   0  0  0  0  0  0  0   0   0
## a12a14    0  0  0  0  0  0  0  0   0   0   1   1  0  0  0  0  0  0   0   0
##          b12 b13 b14 b15
## a10        0   0   0   0
## a10a11     0   0   0   0
## a10a11a4   0   0   0   0
## a11        0   0   0   0
## a12        0   0   0   0
## a12a14     0   0   0   0

现在我们创建一个函数，用于说明小鼠的样本信息，如下所示：

1
2
3
4
5
6
7
8
9
10
11
library(rafalib)
mypar()
flipt <- function(m) t(m[nrow(m):1,])
myimage <- function(m,...) {
  image(flipt(m),xaxt="n",yaxt="n",...)
  }
myimage(as.matrix(pData(maPooling)),col=c("white","black"),
        xlab="experiments",
        ylab="individuals",
        main="phenoData")

plot of chunk unnamed-chunk-3

需要注意的是，最终我们只对两个品系小鼠的差异基因感兴趣，也就是数据中标为0的品系和1的品系。我们可以对合并样本的技术重复进行检验或来源于12只小鼠的数据进行检验。我们可以提取出这些合并的样本，因为来自于每个品系的所有小鼠都都可以用0或1来表示，因此那些来源于某品系（例如用1表示的这个品系）代表了这些样本，因此我们提取出矩阵行之和为12的基因，这段话不太理解，原文如下：

Note that ultimately we are interested in detecting genes that are differentially expressed between the two strains of mice which we will refer to as strain 0 and 1. We can apply tests to the technical replicates of pooled samples or the data from 12 individual mice. We can identify these pooled samples because all mice from each strain were represented in these samples and thus the sum of the rows of experimental design matrix add up to 12:

如下所示：

1
2
3
data(maPooling)
pd=pData(maPooling)
pooled=which(rowSums(pd)==12)

我们可以从列名确定品系，如下所示：

1
2
3
factor(as.numeric(grepl("b",names(pooled))))
## [1] 0 0 0 0 1 1 1 1
## Levels: 0 1

如果我们要比较每个基因在不同品系之间的差异，我们就会发现有几个基因出现了比较一致的差异，如下所示：

1
2
3
4
5
6
7
8
9
###look at 2 pre-selected genes for illustration
i=11425;j=11878
pooled_y=exprs(maPooling[,pooled])
pooled_g=factor(as.numeric(grepl("b",names(pooled))))
mypar(1,2)
stripchart(split(pooled_y[i,],pooled_g),vertical=TRUE,method="jitter",col=c(1,2),
           main="Gene 1",xlab="Group",pch=15)
stripchart(split(pooled_y[j,],pooled_g),vertical=TRUE,method="jitter",col=c(1,2),
           main="Gene 2",xlab="Group",pch=15)

plot of chunk unnamed-chunk-6

这里要注意一下，如果我们根据这些值来进行t检验，那么我们就会得到非常显著的结果，也就是说p值极低，如下所示：

1
2
3
4
5
6
library(genefilter)
pooled_tt=rowttests(pooled_y,pooled_g)
pooled_tt$p.value[i]
## [1] 2.075617e-07
pooled_tt$p.value[j]
## [1] 3.400476e-07

但是，如果我们再次从两个品系的小鼠中各选择12个小鼠的话，这样的结论会成立吗？请注意，t检验的定义包括要比较的总体的标准差。我们这里做的计算结果是否能够反映出总体差异？（注：这里要回顾一下t检验的定义，样本与方差的差异等概念）。

需要注意的是，我们现在重复计算的过程就是在重复原来的实验流程。我们为每个合并的样本创建了4个技术重复。或许Gene 1在小鼠品系内存在着很高的变异，而Gene 2的变异则比较低，但是我们不能看到这一点，因为小鼠之间的变异被合并后的数据所掩盖了（注：这里作者只是假设了Gene 1可能存在着比较高的组内变异）。

我们还有每只小鼠的微阵列数据。对于每个品系的小鼠来说，我们有12个生物学重复。我们可以通过查询1来找到这个品系，如下所示：

1
individuals=which(rowSums(pd)==1)

事实证明，某些单独的小鼠本身就含有技术重复，因此现在我们将这种情况剔除掉，从而用于说明仅用生物学重复的分析结果，如下所示： 

1
2
3
4
##remove replicates
individuals=individuals[-grep("tr",names(individuals))]
y=exprs(maPooling)[,individuals]
g=factor(as.numeric(grepl("b",names(individuals))))

我们可以计算每个基因的在某个品系中的方差，并与通过技术重复获得的标准差进行比较，如下所示：

1
2
3
4
5
technicalsd <- rowSds(pooled_y[,pooled_g==0])
biologicalsd <- rowSds(y[,g==0])
LIM=range(c(technicalsd,biologicalsd))
mypar(1,1)
boxplot(technicalsd,biologicalsd,names=c("technical","biological"),ylab="standard deviation")

plot of chunk unnamed-chunk-10

从上图上我们可以看到，生物学方差(variance)远大于技术方差。并且对于方差的变异而言，生物学方差更大。这只是我们前面提到的2个基因的情况，现在我们展示每只小鼠上测得的这2个基因的表达值，如下所示：

1
2
3
4
5
mypar(1,2)
stripchart(split(y[i,],g),vertical=TRUE,method="jitter",col=c(1,2),xlab="Gene 1",pch=15)
points(c(1,2),tapply(y[i,],g,mean),pch=4,cex=1.5)
stripchart(split(y[j,],g),vertical=TRUE,method="jitter",col=c(1,2),xlab="Gene 2",pch=15)
points(c(1,2),tapply(y[j,],g,mean),pch=4,cex=1.5)

plot of chunk unnamed-chunk-11

现在p值就发生了一点变化，如下所示：

1
2
3
4
5
6
library(genefilter)
tt=rowttests(y,g)
tt$p.value[i]
## [1] 0.0898726
tt$p.value[j]
## [1] 1.979172e-07

当我们在不同品系的小鼠之间比较这2个基因的差异时，我们对哪个基因的差异更有信心呢（Gene 1还是Gene 2）？如果另一位研究人员另外一批随机样本来做实验，你认识他会发现哪个基因有差异（Gene 1还是Gene 2）？如果我们希望我们的结论是关于小鼠品系的通用结论（通用结论这里我的理解就是，无论放哪个实验室做，都能做出这个结果，得到这个结论），而不仅仅是我们这次实验的结论（也就是说，有可能就是只有自己实验室能做出这个结果，那么这个结果就不太可信，有可能是假阳性），那么检测生物学的变异则非常重要。

在这一节中，我们使用了两种分析方法。第一种分析方法是将这两个品系的生物学重复做为总体总体来进行分析析（也就是说把某个品系的12只小鼠数据汇总起来分析）。第二种则是将合并后的数据来作为总体进行分析（我的理解就是把12只小鼠的RNA混合在一起，作为一个样本，然后使用4个技术重复来进行分析，这种分析方法其实检测的就是技术变异）。在科学研究中，我们通常关注的是总体（也就是第一种分析方法）。作为一个非常实际的例子，这里需要注意的是，如果另外一个实验室重复该实验，他们会有另外一组12个小鼠的样本，因此关于我们实验中对总体的推断也会在另外的实验室中得到重复。

An analysis with biological replicates has as a population these two strains of mice. An analysis with technical replicates has as a population the twelve mice we selected and the variability is related to
the measurement technology. In science we typically are concerned with populations. As a very practical example, note that if another lab performs this experiment they will have another set of twelve mice and thus inferences about populations are more likely to be reproducible.

参考资料

  1. Biological versus technical variability

发表于 | 分类于
| 字数统计: 2,258 | 阅读时长 ≈ 11

GSEA图

这篇教程翻译自生信博客DAVE TANGE’S BLOG(https://davetang.org/muse/)上的一篇教程，参考资料中已经列出了地址。

fgsea这个包用于做GSEA分析，先来看一下使用这个包做的图，如下所示：

现在简单解释一下这个图形：

x轴——排序后的基因列表L位置对应的坐标，也就是我们自己通过RNA-seq，芯片，qPCR等手段获得的基因表达值倍数变化，或p值排序，总之，这是一个有序列表。

垂直的黑色细胞——上图中类似条形码的图形，这是指的是某一个基因集S中基因对应于L基因中，的位置，在上图中，这个基因集是细胞周期(Cell Cycle)，明天看到S中的成员在L的左侧比较密集。

y轴——富集分布，从上面我们可以看到，细胞周期(Cell Cycle)这个基因集在左侧富集，也就是绿色曲线表示的位置。

fgsea使用

安装

先安装fgsea包，如下所示：

1
2
3
install.packages("BiocManager")
BiocManager::install("fgsea")
library(fgsea)

数据

fgsea包中内置的有数据集examplePathways，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
library(fgsea)
data(examplePathways)
# examplePathways是从'reactome.db'包中
# 提取的信息，并以列表的形式存了这些信息，
# 这些信息主要是小鼠的基因
help("examplePathways")
# 这个数据集是一个列表
class(examplePathways)
# 一共有1,457 "通路"
length(examplePathways)
# 第1列含有Meiotic_Synapsis pathway的EntrezID
examplePathways[1]
# exampleRanks储存的是排序信息
data(exampleRanks)
# exampleRnak是一个数字型的向量
class(exampleRanks)
# exampleRanks中的向量名称则是对称的Entrez ID
head(exampleRanks)
tail(unname(exampleRanks))

运行结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
> library(fgsea)
> 
> data(examplePathways)
> 
> # examplePathways是从'reactome.db'包中
# 提取的信息，并以列表的形式
> # 储存了这些信息，这些信息主要是小鼠的基因
> help("examplePathways")
> 
> # 这个数据集是一个列表
> class(examplePathways)
[1] "list"
> 
> # 一共有1,457 "通路"
> length(examplePathways)
[1] 1457
> 
> # 第1列含有Meiotic_Synapsis pathway的EntrezID
> examplePathways[1]
$`1221633_Meiotic_Synapsis`
 [1] "12189"     "13006"     "15077"     "15078"     "15270"     "15512"    
 [7] "16905"     "16906"     "19357"     "20842"     "20843"     "20957"    
[13] "20962"     "21749"     "21750"     "22196"     "23856"     "24061"    
[19] "28113"     "50878"     "56739"     "57321"     "64009"     "66654"    
[25] "69386"     "71846"     "74075"     "77053"     "94244"     "97114"    
[31] "97122"     "97908"     "101185"    "140557"    "223697"    "260423"   
[37] "319148"    "319149"    "319150"    "319151"    "319152"    "319153"   
[43] "319154"    "319155"    "319156"    "319157"    "319158"    "319159"   
[49] "319160"    "319161"    "319565"    "320332"    "320558"    "326619"   
[55] "326620"    "360198"    "497652"    "544973"    "625328"    "667250"   
[61] "100041230" "102641229" "102641751" "102642045"
> 
> # exampleRanks储存的是排序信息
> data(exampleRanks)
> 
> # exampleRnak是一个数字型的向量
> class(exampleRanks)
[1] "numeric"
> 
> # exampleRanks中的向量名称则是对称的Entrez ID
> head(exampleRanks)
   170942    109711     18124     12775     72148     16010 
-63.33703 -49.74779 -43.63878 -41.51889 -33.26039 -32.77626 
> tail(unname(exampleRanks))
[1] 47.58235 49.87543 50.25179 50.86532 51.16110 53.28400

分析

现在我们使用上面的数据进行GSEA分析，进行GSEA分析时，我们的通路文件（也就是GSEA官网中的GMT文件，对应的就是基因集S）储存在pathways参数中，用户自己的数据（排过序的数据）放在stats数据集中，剩下的参数不用管，如下所示：

1
2
3
4
5
fgseaRes <- fgsea(pathways = examplePathways, 
                  stats = exampleRanks,
                  minSize=15,
                  maxSize=500,
                  nperm=100000)

分析的结果fgseaRes是一个data.table格式的文件，使用plotEnrichment函数可以绘制出GSEA分析的结果，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
library(ggplot2)
# 参考 https://github.com/Rdatatable/data.table/wiki
# 来查看data.table相关的信息
class(fgseaRes)
# 查看p值小于0.01的通路的数目
sum(fgseaRes[, padj < 0.01])
# 绘制出最显著的富集通路
plotEnrichment(examplePathways[[head(fgseaRes[order(pval), ], 1)$pathway]],
               exampleRanks) + 
               labs(title=head(fgseaRes[order(pval), ], 1)$pathway)

结果如下所示：

1
2
3
4
5
6
> class(fgseaRes)
[1] "data.table" "data.frame"
> 
> # 查看p值小于0.01的通路的数目
> sum(fgseaRes[, padj < 0.01])
[1] 78

可以画出具体的某条通路，如下所示；

1
2
3
plotEnrichment(examplePathways[["5991130_Programmed_Cell_Death"]],
               exampleRanks) + 
  labs(title="Programmed Cell Death")

还可以在一张图中绘制出前10个富集通路，以及后10个富集通路，一共20个，如下所示：

1
2
3
4
5
topPathwaysUp <- fgseaRes[ES > 0][head(order(pval), n=10), pathway]
topPathwaysDown <- fgseaRes[ES < 0][head(order(pval), n=10), pathway]
topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
plotGseaTable(examplePathways[topPathways], exampleRanks, fgseaRes, 
              gseaParam = 0.5)

Reactome

也可以使用Reactome通路来进行GSEA分析，此时需要安装reactome.db包，这个包很大，600多M，下载的时间很长，安装过程如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# BiocManager::install("reactome.db")
library(reactome.db)
my_pathways <- reactomePathways(names(exampleRanks))
 
# Reactome pathways have a median of 11 genes
summary(sapply(my_pathways, length))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0     4.0    11.0    30.3    30.0  1140.0
 
fgsea_reactome <- fgsea(pathways = my_pathways, 
                        stats = exampleRanks,
                        minSize=15,
                        maxSize=500,
                        nperm=100000)
 
head(fgsea_reactome[order(pval), ])
sum(fgsea_reactome[, padj < 0.01])

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
> head(fgsea_reactome[order(pval), ])
                    pathway         pval         padj
1:               Cell Cycle 1.212856e-05 0.0002571089
2:      Cell Cycle, Mitotic 1.236552e-05 0.0002571089
3: Neutrophil degranulation 1.274974e-05 0.0002571089
4: Signaling by Rho GTPases 1.291339e-05 0.0002571089
5:                  M Phase 1.303509e-05 0.0002571089
6:   Cell Cycle Checkpoints 1.350749e-05 0.0002571089
          ES      NES nMoreExtreme size
1: 0.5324037 2.674889            0  414
2: 0.5475346 2.719737            0  363
3: 0.4258088 2.074081            0  296
4: 0.4073481 1.966039            0  271
5: 0.5022189 2.408762            0  255
6: 0.6085986 2.836881            0  200
                                  leadingEdge
1:   66336,66977,15366,12442,107995,66442,...
2:   66336,66977,15366,12442,107995,66442,...
3:    11676,14190,53381,12306,20430,12505,...
4: 66336,66977,20430,104215,233406,107995,...
5:   66336,66977,12442,107995,66442,52276,...
6:   66336,66977,12442,107995,66442,12428,...
> sum(fgsea_reactome[, padj < 0.01])
[1] 103

Leading edge

在GSEA分析中，我们通常会提取那些构成高得分的核心基因。我们对高得分的核心基因的定义就是，基因集S中位于排序基因列表L位置中的得分最大处之前或之后的那些基因集（也就是GSEA结果中绿色曲线最高点的前面或后面）。

前面我们注意到，GSEA分析的结果中有1列被命名为leadingEdge。这一列包含了高得分的基因。我们使用Reactome通路的富集结果来提取这些基因，如下所示：

1
2
3
4
5
6
7
8
9
10
11
# the most significant pathway
fgsea_reactome[order(pval),][1,]
# list of Entrez gene IDs that contributed to the enrichment score
fgsea_reactome[order(pval),][1,]$leadingEdge
# how many genes are in the leading edge?
length(fgsea_reactome[order(pval),][1,]$leadingEdge[[1]])
# how many genes are in the Cell Cycle pathway?
length(my_pathways[['Cell Cycle']])

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
> # the most significant pathway
> fgsea_reactome[order(pval),][1,]
      pathway         pval         padj        ES
1: Cell Cycle 1.212856e-05 0.0002571089 0.5324037
        NES nMoreExtreme size
1: 2.674889            0  414
                                leadingEdge
1: 66336,66977,15366,12442,107995,66442,...
> 
> # list of Entrez gene IDs that contributed to the enrichment score
> fgsea_reactome[order(pval),][1,]$leadingEdge
[[1]]
  [1] "66336"  "66977"  "15366"  "12442"  "107995"
  [6] "66442"  "12571"  "12428"  "52276"  "54392" 
 [11] "66311"  "215387" "67629"  "12649"  "72415" 
 [16] "56150"  "57441"  "20877"  "67121"  "12615" 
 [21] "11799"  "66468"  "67849"  "19053"  "73804" 
 [26] "76044"  "20878"  "15270"  "13555"  "60411" 
 [31] "12580"  "17219"  "69270"  "12575"  "69263" 
 [36] "12448"  "14211"  "20873"  "18005"  "72119" 
 [41] "71988"  "12189"  "17215"  "12534"  "66156" 
 [46] "208628" "237911" "22390"  "68240"  "228421"
 [51] "68014"  "269582" "19348"  "12236"  "72151" 
 [56] "18817"  "21781"  "18968"  "217653" "66934" 
 [61] "272551" "227613" "67141"  "67951"  "68612" 
 [66] "68298"  "108000" "23834"  "106344" "56452" 
 [71] "242705" "18141"  "223921" "26886"  "13557" 
 [76] "26909"  "72787"  "268697" "72155"  "56371" 
 [81] "17535"  "107823" "12531"  "327762" "12567" 
 [86] "229841" "67052"  "16319"  "66634"  "171567"
 [91] "26931"  "67203"  "12235"  "19891"  "74470" 
 [96] "72083"  "381318" "66570"  "17216"  "76308" 
[101] "19687"  "17218"  "102920" "29870"  "18973" 
[106] "16881"  "17463"  "75786"  "19645"  "19075" 
[111] "26417"  "69736"  "19357"  "76816"  "70385" 
[116] "70645"  "22628"  "225182" "22627"  "52683" 
[121] "19076"  "18972"  "231863" "26932"  "12544" 
[126] "17997"  "51788"  "26440"  "68549"  "12445" 
[131] "19088"  "269113" "26444"  "19324"  "103733"
[136] "59001"  "107976" "19179"  "12579"  "232987"
[141] "17420"  "228769" "219072" "26445"  "105988"
[146] "69745"  "18538"  "69928"  "11651"  "235559"
[151] "68097"  "57296"  "63955"  "14235"  "19170" 
[156] "17246"  "17220"  "12144"  "50793"  "77605" 
[161] "18392"  "236930" "67151"  "70024"  "59126" 
[166] "66296"  "16906"  "109145" "71819"  "67733" 
[171] "50883"  "12447"  "12532"  "14156"  "26442" 
[176] "19177"  "230376" "245688"
> 
> # how many genes are in the leading edge?
> length(fgsea_reactome[order(pval),][1,]$leadingEdge[[1]])
[1] 178
> 
> # how many genes are in the Cell Cycle pathway?
> length(my_pathways[['Cell Cycle']])
[1] 414

总结

GSEA是于2005年首次提出来的，现在已经成了基因表达分析中的常规分析手段，它不同于GO分析，GO分析只关注差异基因，而GSEA分析则关注所有的基因。fgsea包可以使用预先排列好的基因一R中进行GSEA分析。p值的计算结果是基于置换检验(permutation test)，这种方法并不是十分精我，因为它忽略了基因之间的相关性，有可能会导致假阳性。但是，在这种方法在研究上调与上调基因方面还是很有用的。即使你计算出的GSEA结果中，p值大于0.05，但是是也可以参考leading edge基因集，为你的实验进行指导。

案例分析

作者提供了用于生成类似于exampleRanks文件的R脚本，不过使用的GEO的数据，平时自己利用fgsea包进行GSEA分析时，生成就好，现在看一下如何将GEO的数据生成类似于exampleRanks文件的排序信息，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
source("https://raw.githubusercontent.com/assaron/r-utils/master/R/exprs.R")
library(GEOquery)
library(limma)
gse14308 <- getGEO("GSE14308")[[1]]
pData(gse14308)$condition <- sub("-.*$", "", gse14308$title)
es <- collapseBy(gse14308, fData(gse14308)$ENTREZ_GENE_ID, FUN=median)
es <- es[!grepl("///", rownames(es)), ]
es <- es[rownames(es) != "", ]
exprs(es) <- normalizeBetweenArrays(log2(exprs(es)+1), method="quantile")
es.design <- model.matrix(~0+condition, data=pData(es))
fit <- lmFit(es, es.design)
fit2 <- contrasts.fit(fit, makeContrasts(conditionTh1-conditionNaive,
                                         levels=es.design))
fit2 <- eBayes(fit2)
de <- data.table(topTable(fit2, adjust.method="BH", number=12000, sort.by = "B"), keep.rownames = T)
ranks <- de[order(t), list(rn, t)]
head(ranks)
tail(ranks)

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
> head(ranks)
       rn         t
1: 109711 -50.99316
2:  18124 -44.40559
3:  12775 -43.22109
4:  72148 -33.74441
5:  16010 -33.34034
6:  16206 -30.67962
> tail(ranks)
      rn        t
1: 80901 47.55618
2: 58801 48.53756
3: 15937 49.88068
4: 13730 50.17693
5: 12772 50.32302
6: 80876 51.72987

现在绘制一下上述数据中6个上升与6个下降的基因热图：

1
2
3
4
5
6
7
8
9
10
11
12
library(pheatmap)
 
my_group <- data.frame(group = pData(es)$condition)
row.names(my_group) <- colnames(exprs(es))
 
pheatmap(
  mat = es[c(head(de[order(t), 1])$rn, tail(de[order(t), 1])$rn),],
  annotation_col = my_group,
  cluster_rows = FALSE,
  cellwidth=25,
  cellheight=15
)

参考资料

  1. Using fgsea package
  2. Comparison of clusterProfiler and GSEA-P
  3. Using the fast preranked gene set enrichment analysis (fgsea) package

发表于 | 分类于
| 字数统计: 6,582 | 阅读时长 ≈ 31

前言

这一部分内容涉及R中使用人类基因且，内含子，外显子，转录本，AnnotationHub，基因组的注释包，GO分析，KEGG分析等，笔记末尾的参考文献是原文。

基础注释资源与发现

在这一部分里，我们将回顾Bioconductor中用于处理和注释基因组序列的一些工具。我们将研究参考基因组序列，转录本和基因，并以基因通路(gene pathway)作为结束。我们学习这一部分的最终目标就是使用注释信息来帮助我们对基因组实验进行可靠的解释。Bioconductor的基本目标就是更加方便地有关基因组结构和功能的信息统计统计分析程序。

注释概念的层次结构

Bioconductor包括许多不同类型的基因组注释。我们可以在层次结构中来理解这些注释资源。

  • 最基因的注释就是某个物种的参考基因组序列。它总是按照核苷酸的线性方式排列成染色体（例如参考基因组。
  • 在此之上就是将染色体序列排列到感兴趣的区域中。最感兴趣的区域就是基因，但是注释中也含有其它的信息，例如SNP或CpG位点。基因具有内部结构，即被转录的部分和未被转录的部分。“基因模式”定义了在基因组坐标中的标记和布置这些结构的方式。
    • 在感兴趣的区域(regions of interest)的理念下，我们还定义了面向平台的注释(platform-oriented annotation)。这处类型的注释通常首先是由厂家提供的，但随着研究的进行，对这些平台中最初有歧义信息进行了确认和更新，从而完善了这些注释内容。密歇根大学的brainarray project 说明了affymetrix阵列注释的过程。我们将在本节最后讨论面向平台注释的问题。
  • 在此之是是将区域（通常是基因或基因的产物）组成成具有共同结构或功能特性的组。例如在细胞中共同被发现的，或者是被鉴定为在生物学过程中协同作用的基因组（我的理解就是GO分析，KEGG分析这一类）。

发现可用的参考基因组

Bioconductor已经包含了注释包的合成，将它这一层次结构上的所有元素都带了可编程环境中。参考基因组序列是使用Biostrings和BSgenome包中的工具进行管理的，available.genomes 函数能够列出构建好的人和现在各种模式生物的参考基因组，如下所示：

1
2
3
4
5
6
7
8
9
10
11
library(Biostrings)
ag = available.genomes()
length(ag)
## [1] 87
head(ag)
## [1] "BSgenome.Alyrata.JGI.v1"                
## [2] "BSgenome.Amellifera.BeeBase.assembly4"  
## [3] "BSgenome.Amellifera.UCSC.apiMel2"       
## [4] "BSgenome.Amellifera.UCSC.apiMel2.masked"
## [5] "BSgenome.Athaliana.TAIR.04232008"       
## [6] "BSgenome.Athaliana.TAIR.TAIR9"

参考基因组的版本很重要

不同物种的参考基因组是从头构建的，然后随着算法和测序数据的不断改进而进一步完善。对人类而言，基因组研究联盟（Genome Research Consortium）于2009年构建了37号版本，并于2013年构建了38号版本。

一旦参考基因组构建完成，就哦可以很轻松地对某个物种进行信息丰富的基因组序列分析，因为人们可以专注于那引起已知含有等位基因多样性的区域。

The reference build for an organism is created de novo and then refined as algorithms and sequenced data improve. For humans, the Genome Research Consortium signed off on build 37 in 2009, and on build 38 in 2013.

需要注意的是，基因组序列包含有很长的名称，这个名称里包括版本信息。这样命名的方式就是为了避免与不同版本的参考基因组混淆。在LiftOver这节视频里，我们就展示了如何使用UCSC的liftOver工具与rtracklayer包中的接口对接，从而实现不同版本的基因组坐标转化的过程。

为了帮助用户避免混淆从不同参考基因组坐标上收集分析来的数据，我们提供了一个”基因组“标签，这个标签填充了大多关于序列的信息。在随后的部分里，我们会看到一些案例。用于序列比对的软件可以检查被比对上的序列的兼容标签，从而有助于确保有意义的结果。

H. sapiens的参考基因序列

通过安装并添加一个单独的R包就能获取智人(Homo sapiens)的参考序列。这个程序包定义了一个Hsapiens对象，试剂公司对象是染色体序列的来源，但是当对其进行单独显示时，它会提供相关序列数据来源的信息，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
library(BSgenome.Hsapiens.UCSC.hg19)
Hsapiens
## Human genome:
## # organism: Homo sapiens (Human)
## # provider: UCSC
## # provider version: hg19
## # release date: Feb. 2009
## # release name: Genome Reference Consortium GRCh37
## # 93 sequences:
## #   chr1                  chr2                  chr3                 
## #   chr4                  chr5                  chr6                 
## #   chr7                  chr8                  chr9                 
## #   chr10                 chr11                 chr12                
## #   chr13                 chr14                 chr15                
## #   ...                   ...                   ...                  
## #   chrUn_gl000235        chrUn_gl000236        chrUn_gl000237       
## #   chrUn_gl000238        chrUn_gl000239        chrUn_gl000240       
## #   chrUn_gl000241        chrUn_gl000242        chrUn_gl000243       
## #   chrUn_gl000244        chrUn_gl000245        chrUn_gl000246       
## #   chrUn_gl000247        chrUn_gl000248        chrUn_gl000249       
## # (use 'seqnames()' to see all the sequence names, use the '$' or '[['
## # operator to access a given sequence)
head(genome(Hsapiens))  # see the tag
##   chr1   chr2   chr3   chr4   chr5   chr6 
## "hg19" "hg19" "hg19" "hg19" "hg19" "hg19"

我们使用 $ 符号来获取17号染色体的序列，如下所示：

1
2
3
Hsapiens$chr17
##   81195210-letter "DNAString" instance
## seq: AAGCTTCTCACCCTGTTCCTGCATAGATAATTGC...GGTGTGGGTGTGGTGTGTGGGTGTGGGTGTGGT

参考序列的转录本和基因

UCSC注释

TxDb包家族和数据对象管理了转录本和基因模式信息。我们可以认为这些信息来源于UCSC基因组浏览器的注释表，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene # abbreviate
txdb
## TxDb object:
## # Db type: TxDb
## # Supporting package: GenomicFeatures
## # Data source: UCSC
## # Genome: hg19
## # Organism: Homo sapiens
## # Taxonomy ID: 9606
## # UCSC Table: knownGene
## # Resource URL: http://genome.ucsc.edu/
## # Type of Gene ID: Entrez Gene ID
## # Full dataset: yes
## # miRBase build ID: GRCh37
## # transcript_nrow: 82960
## # exon_nrow: 289969
## # cds_nrow: 237533
## # Db created by: GenomicFeatures package from Bioconductor
## # Creation time: 2015-10-07 18:11:28 +0000 (Wed, 07 Oct 2015)
## # GenomicFeatures version at creation time: 1.21.30
## # RSQLite version at creation time: 1.0.0
## # DBSCHEMAVERSION: 1.1

我们使用 genes() 来获取Entrez Gene ID的地址，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
ghs = genes(txdb)
ghs
## GRanges object with 23056 ranges and 1 metadata column:
##         seqnames                 ranges strand |     gene_id
##            <Rle>              <IRanges>  <Rle> | <character>
##       1    chr19 [ 58858172,  58874214]      - |           1
##      10     chr8 [ 18248755,  18258723]      + |          10
##     100    chr20 [ 43248163,  43280376]      - |         100
##    1000    chr18 [ 25530930,  25757445]      - |        1000
##   10000     chr1 [243651535, 244006886]      - |       10000
##     ...      ...                    ...    ... .         ...
##    9991     chr9 [114979995, 115095944]      - |        9991
##    9992    chr21 [ 35736323,  35743440]      + |        9992
##    9993    chr22 [ 19023795,  19109967]      - |        9993
##    9994     chr6 [ 90539619,  90584155]      + |        9994
##    9997    chr22 [ 50961997,  50964905]      - |        9997
##   -------
##   seqinfo: 93 sequences (1 circular) from hg19 genome

我们也可以使用合适的标识符进行信息过滤。现在我们提取两个不同基因的外显子，这些外显子由其Entrez基因ID标明，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
exons(txdb, columns=c("EXONID", "TXNAME", "GENEID"),
                  filter=list(gene_id=c(100, 101)))
## GRanges object with 39 ranges and 3 metadata columns:
##        seqnames                 ranges strand |    EXONID
##           <Rle>              <IRanges>  <Rle> | <integer>
##    [1]    chr10 [135075920, 135076737]      - |    144421
##    [2]    chr10 [135077192, 135077269]      - |    144422
##    [3]    chr10 [135080856, 135080921]      - |    144423
##    [4]    chr10 [135081433, 135081570]      - |    144424
##    [5]    chr10 [135081433, 135081622]      - |    144425
##    ...      ...                    ...    ... .       ...
##   [35]    chr20   [43254210, 43254325]      - |    256371
##   [36]    chr20   [43255097, 43255240]      - |    256372
##   [37]    chr20   [43257688, 43257810]      - |    256373
##   [38]    chr20   [43264868, 43264929]      - |    256374
##   [39]    chr20   [43280216, 43280376]      - |    256375
##                                  TXNAME          GENEID
##                         <CharacterList> <CharacterList>
##    [1] uc009ybi.3,uc010qva.2,uc021qbe.1             101
##    [2]            uc009ybi.3,uc021qbe.1             101
##    [3] uc009ybi.3,uc010qva.2,uc021qbe.1             101
##    [4]                       uc009ybi.3             101
##    [5]            uc010qva.2,uc021qbe.1             101
##    ...                              ...             ...
##   [35]                       uc002xmj.3             100
##   [36]                       uc002xmj.3             100
##   [37]                       uc002xmj.3             100
##   [38]                       uc002xmj.3             100
##   [39]                       uc002xmj.3             100
##   -------
##   seqinfo: 93 sequences (1 circular) from hg19 genome

ENSEMBL注释

Ensembl home主页上写道：Ensembl创建，整合和发布研究基因组的参考数据库和工具。该项目位于 欧洲分子生物学实验室，该实验室的数据库支持其注释资源可以与Bioconductor兼容。

ensembldb 包含有一个简要说明，其内容如下所示：

ensembldb包提供了一些函数，这些函数用于创建和使用以转录本为中心的注释数据库/包。使用注释数据库的Perl API可以从Ensembl 1中直接获取这些数据。TxDb 包的功能和数据类似于GenomicFeatures包，另外，除了从数据库检索所有的基因/转录本模型和注释外，ensembldb包还提供了一个过滤框架，用于检索特定条目的注释，例如位于染色体区域上的某编码基因或某LincRNA转录模式的特定条目。从1.7版本开始，由ensembldb创建的EnsDb数据库还包含蛋白质注释数据库（参考第11节：数据库而已和可用属性/列的概述）。有关蛋白质注释的信息请参考蛋白质的vignette，如下所示：

1
2
3
4
5
6
library(ensembldb)
library(EnsDb.Hsapiens.v75)
names(listTables(EnsDb.Hsapiens.v75))
##  [1] "gene"           "tx"             "tx2exon"        "exon"          
##  [5] "chromosome"     "protein"        "uniprot"        "protein_domain"
##  [9] "entrezgene"     "metadata"

举例说明如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
edb = EnsDb.Hsapiens.v75  # abbreviate
txs <- transcripts(edb, filter = GenenameFilter("ZBTB16"),
                   columns = c("protein_id", "uniprot_id", "tx_biotype"))
txs
## GRanges object with 20 ranges and 5 metadata columns:
##                   seqnames                 ranges strand |      protein_id
##                      <Rle>              <IRanges>  <Rle> |     <character>
##   ENST00000335953       11 [113930315, 114121398]      + | ENSP00000338157
##   ENST00000335953       11 [113930315, 114121398]      + | ENSP00000338157
##   ENST00000335953       11 [113930315, 114121398]      + | ENSP00000338157
##   ENST00000335953       11 [113930315, 114121398]      + | ENSP00000338157
##   ENST00000335953       11 [113930315, 114121398]      + | ENSP00000338157
##               ...      ...                    ...    ... .             ...
##   ENST00000392996       11 [113931229, 114121374]      + | ENSP00000376721
##   ENST00000539918       11 [113935134, 114118066]      + | ENSP00000445047
##   ENST00000545851       11 [114051488, 114118018]      + |            <NA>
##   ENST00000535379       11 [114107929, 114121279]      + |            <NA>
##   ENST00000535509       11 [114117512, 114121198]      + |            <NA>
##                     uniprot_id              tx_biotype           tx_id
##                    <character>             <character>     <character>
##   ENST00000335953  ZBT16_HUMAN          protein_coding ENST00000335953
##   ENST00000335953 Q71UL7_HUMAN          protein_coding ENST00000335953
##   ENST00000335953 Q71UL6_HUMAN          protein_coding ENST00000335953
##   ENST00000335953 Q71UL5_HUMAN          protein_coding ENST00000335953
##   ENST00000335953 F5H6C3_HUMAN          protein_coding ENST00000335953
##               ...          ...                     ...             ...
##   ENST00000392996 F5H5Y7_HUMAN          protein_coding ENST00000392996
##   ENST00000539918         <NA> nonsense_mediated_decay ENST00000539918
##   ENST00000545851         <NA>    processed_transcript ENST00000545851
##   ENST00000535379         <NA>    processed_transcript ENST00000535379
##   ENST00000535509         <NA>         retained_intron ENST00000535509
##                     gene_name
##                   <character>
##   ENST00000335953      ZBTB16
##   ENST00000335953      ZBTB16
##   ENST00000335953      ZBTB16
##   ENST00000335953      ZBTB16
##   ENST00000335953      ZBTB16
##               ...         ...
##   ENST00000392996      ZBTB16
##   ENST00000539918      ZBTB16
##   ENST00000545851      ZBTB16
##   ENST00000535379      ZBTB16
##   ENST00000535509      ZBTB16
##   -------
##   seqinfo: 1 sequence from GRCh37 genome

你的数据将会成他人的注释：导入/导出

ENCODE项目很地说明了今天的实验是明天的注释。你应该以同样的方式考虑自己的实验（当然，要使实验成为可靠且持久的注释，它必须解决有关基因组结构或功能的重要问题，并且必须使用适当的，能正确执行的实验流程。需要注意，ENCODE能够非常明确地将实验流程链接到数据）。

例如，我们来看一个雌激素受体结合数据，它是由ENCODE发布的一个narrowPeak 数据。它的碱基是用ascii文本表示的，因此可以很容易地导入为一组文本数据。如果记录的字段有一定的规律性，则可以将文件作为表格导入。

但是，我们不仅是想导入数据，还想将导入的数据作为可计算的对象。我们认识到arrowePeak和bedGraph格式之间的联系后，我们就可以立即将其导入GRanges中。

为了说明这一点，我们在ERBS包中找到narrowPeak原始数据文件的路径，如下所示：

1
2
3
4
5
6
f1 = dir(system.file("extdata",package="ERBS"), full=TRUE)[1]
readLines(f1, 4) # look at a few lines
## [1] "chrX\t1509354\t1512462\t5\t0\t.\t157.92\t310\t32.000000\t1991"    
## [2] "chrX\t26801421\t26802448\t6\t0\t.\t147.38\t310\t32.000000\t387"   
## [3] "chr19\t11694101\t11695359\t1\t0\t.\t99.71\t311.66\t32.000000\t861"
## [4] "chr19\t4076892\t4079276\t4\t0\t.\t84.74\t310\t32.000000\t1508"

使用import命令非常简单，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
library(rtracklayer)
imp = import(f1, format="bedGraph")
imp
## GRanges object with 1873 ranges and 7 metadata columns:
##          seqnames               ranges strand |     score       NA.
##             <Rle>            <IRanges>  <Rle> | <numeric> <integer>
##      [1]     chrX [ 1509355,  1512462]      * |         5         0
##      [2]     chrX [26801422, 26802448]      * |         6         0
##      [3]    chr19 [11694102, 11695359]      * |         1         0
##      [4]    chr19 [ 4076893,  4079276]      * |         4         0
##      [5]     chr3 [53288568, 53290767]      * |         9         0
##      ...      ...                  ...    ... .       ...       ...
##   [1869]    chr19 [11201120, 11203985]      * |      8701         0
##   [1870]    chr19 [ 2234920,  2237370]      * |       990         0
##   [1871]     chr1 [94311336, 94313543]      * |      4035         0
##   [1872]    chr19 [45690614, 45691210]      * |     10688         0
##   [1873]    chr19 [ 6110100,  6111252]      * |      2274         0
##               NA.1      NA.2      NA.3      NA.4      NA.5
##          <logical> <numeric> <numeric> <numeric> <integer>
##      [1]      <NA>    157.92       310        32      1991
##      [2]      <NA>    147.38       310        32       387
##      [3]      <NA>     99.71    311.66        32       861
##      [4]      <NA>     84.74       310        32      1508
##      [5]      <NA>      78.2   299.505        32      1772
##      ...       ...       ...       ...       ...       ...
##   [1869]      <NA>      8.65     7.281   0.26576      2496
##   [1870]      <NA>      8.65    26.258  1.995679      1478
##   [1871]      <NA>      8.65    12.511   1.47237      1848
##   [1872]      <NA>      8.65     6.205         0       298
##   [1873]      <NA>      8.65    17.356  2.013228       496
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths
genome(imp)  # genome identifier tag not set, but you should set it
##  chrX chr19  chr3 chr17  chr8 chr11 chr16  chr1  chr2  chr6  chr9  chr7 
##    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA 
##  chr5 chr12 chr20 chr21 chr22 chr18 chr10 chr14 chr15  chr4 chr13 
##    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA

我们可以通过一次获取GRanges。元数据列中还有一些其他字段用于指定名称，但是如果我们只对范围感兴趣，除了添加基因组元数据以防止与不兼容的坐标中记录的数据非法组合外，我们就完成了这个任务（这一段不太理解，原文如下）：

We obtain a GRanges in one stroke. There are some additional fields in the metadata columns whose names should be specified, but if we are interested only in the ranges, we are done, with the exception of adding the genome metadata to protect against illegitimate combination with data recorded in an incompatible coordinate system.

为了与其他得养家或系统进行交流，我们有两个主要选择。我们可以将GRanges保存为RData对象，轻松地传递给另外一个R用户使用。或者，我们歌词采用其他标准格式进行导出。例如，如果我们仅对间隔地址和绑定的得分感兴趣，则仅保存为“bed”格式就足够了，如下所示：

1
2
3
4
5
6
7
export(imp, "demoex.bed")  # implicit format choice
cat(readLines("demoex.bed", n=5), sep="\n")
## chrX	1509354	1512462	.	5	.
## chrX	26801421	26802448	.	6	.
## chr19	11694101	11695359	.	1	.
## chr19	4076892	4079276	.	4	.
## chr3	53288567	53290767	.	9	.

我们已经进行了导入，建模和导入实验数据之间的“往返”，该实验数据可以与其他数据集成在一起，从而增进生物学的理解。

我们需要注意的是，注释在某种程度上是永久正确的，它与在知识边界上的研究进展乏味地隔离开来。我们已经看到了，甚至人类染色体的参考序列也受到了修订。在使用ERBS包时，我们将未知的实验结果视为定义ER结合位点从而进入潜在的生物学解释。不确定性，峰鉴定的可变质量，尚未得到明确估计，但应该是这个样子。

Bioconductor已经尽力致力于这种情况的多个方面。我们维护软件先前版本和注释的存档，以便可以检查或修改过去的工作。我们每年会两次更新中疏注释资源，以确保正在进行的工作以及获得新知识的稳定性。而且，我们已经简化了导入和创建实验数据和注释数据的表示形式。

AnnotationHub

AnnotationHub包用于获取GRanges或其它的合适设计的容器，用于机构设计的容器，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
library(AnnotationHub)
## 
## Attaching package: 'AnnotationHub'
## The following object is masked from 'package:Biobase':
## 
##     cache
ah = AnnotationHub()
## snapshotDate(): 2017-10-27
ah
## AnnotationHub with 42282 records
## # snapshotDate(): 2017-10-27 
## # $$dataprovider: BroadInstitute, Ensembl, UCSC, ftp://ftp.ncbi.nlm.nih....
## # $$species: Homo sapiens, Mus musculus, Drosophila melanogaster, Bos ta...
## # $$rdataclass: GRanges, BigWigFile, FaFile, TwoBitFile, Rle, ChainFile,...
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass,
## #   tags, rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["AH2"]]' 
## 
##             title                                                         
##   AH2     | Ailuropoda_melanoleuca.ailMel1.69.dna.toplevel.fa             
##   AH3     | Ailuropoda_melanoleuca.ailMel1.69.dna_rm.toplevel.fa          
##   AH4     | Ailuropoda_melanoleuca.ailMel1.69.dna_sm.toplevel.fa          
##   AH5     | Ailuropoda_melanoleuca.ailMel1.69.ncrna.fa                    
##   AH6     | Ailuropoda_melanoleuca.ailMel1.69.pep.all.fa                  
##   ...       ...                                                           
##   AH58988 | org.Flavobacterium_piscicida.eg.sqlite                        
##   AH58989 | org.Bacteroides_fragilis_YCH46.eg.sqlite                      
##   AH58990 | org.Pseudomonas_mendocina_ymp.eg.sqlite                       
##   AH58991 | org.Salmonella_enterica_subsp._enterica_serovar_Typhimurium...
##   AH58992 | org.Acinetobacter_baumannii.eg.sqlite

我们可以通过AnnotationHub获得许多与HepG2细胞系相关的实验数据对象，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
query(ah, "HepG2")
## AnnotationHub with 440 records
## # snapshotDate(): 2017-10-27 
## # $$dataprovider: UCSC, BroadInstitute, Pazar
## # $$species: Homo sapiens, NA
## # $$rdataclass: GRanges, BigWigFile
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass,
## #   tags, rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["AH22246"]]' 
## 
##             title                                                         
##   AH22246 | pazar_CEBPA_HEPG2_Schmidt_20120522.csv                        
##   AH22249 | pazar_CTCF_HEPG2_Schmidt_20120522.csv                         
##   AH22273 | pazar_HNF4A_HEPG2_Schmidt_20120522.csv                        
##   AH22309 | pazar_STAG1_HEPG2_Schmidt_20120522.csv                        
##   AH22348 | wgEncodeAffyRnaChipFiltTransfragsHepg2CytosolLongnonpolya.b...
##   ...       ...                                                           
##   AH41564 | E118-H4K5ac.imputed.pval.signal.bigwig                        
##   AH41691 | E118-H4K8ac.imputed.pval.signal.bigwig                        
##   AH41818 | E118-H4K91ac.imputed.pval.signal.bigwig                       
##   AH46971 | E118_15_coreMarks_mnemonics.bed.gz                            
##   AH49484 | E118_RRBS_FractionalMethylation.bigwig

query 方法可以使用过滤字符串的向量。要限制对寻址组蛋白H4K5的注释资源的响应，只需要添加该标签，如下所示(To limit response to annotation resources addressing the histone H4K5, simply add that tag)：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
query(ah, c("HepG2", "H4K5"))
## AnnotationHub with 1 record
## # snapshotDate(): 2017-10-27 
## # names(): AH41564
## # $$dataprovider: BroadInstitute
## # $$species: Homo sapiens
## # $$rdataclass: BigWigFile
## # $$rdatadateadded: 2015-05-08
## # $$title: E118-H4K5ac.imputed.pval.signal.bigwig
## # $$description: Bigwig File containing -log10(p-value) signal tracks fr...
## # $$taxonomyid: 9606
## # $$genome: hg19
## # $$sourcetype: BigWig
## # $$sourceurl: http://egg2.wustl.edu/roadmap/data/byFileType/signal/cons...
## # $$sourcesize: 226630905
## # $$tags: c("EpigenomeRoadMap", "signal", "consolidatedImputed",
## #   "H4K5ac", "E118", "ENCODE2012", "LIV.HEPG2.CNCR", "HepG2
## #   Hepatocellular Carcinoma Cell Line") 
## # retrieve record with 'object[["AH41564"]]'

The OrgDb基因注释图

那些命名为org.*.ge.db 的包含在基因水平上链接到位置，蛋白产物标识符，KEGG途径和GO term，PMIDs以及其它注释资源的标识符的信息，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db) # columns() gives same answer
##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT" 
##  [5] "ENSEMBLTRANS" "ENTREZID"     "ENZYME"       "EVIDENCE"    
##  [9] "EVIDENCEALL"  "GENENAME"     "GO"           "GOALL"       
## [13] "IPI"          "MAP"          "OMIM"         "ONTOLOGY"    
## [17] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
## [21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
## [25] "UNIGENE"      "UNIPROT"
head(select(org.Hs.eg.db, keys="ORMDL3", keytype="SYMBOL", 
   columns="PMID"))
## 'select()' returned 1:many mapping between keys and columns
##   SYMBOL     PMID
## 1 ORMDL3 11042152
## 2 ORMDL3 12093374
## 3 ORMDL3 12477932
## 4 ORMDL3 14702039
## 5 ORMDL3 15489334
## 6 ORMDL3 16169070

基因集和通路资源

基因本体论

Gene Ontology (GO)是一种广泛使用的结构化词汇，它组织了基因和基因产物在以下方面的内容：

  • 生物过程
  • 分子功能
  • 细胞组分。

这套词汇本身旨在与所有生物有关。它采用有向无环图的形式，其中term作为节点，使用is-a和part-of关系作构成了大多数链接。

将生物体特定基因链接到基因本体中的术语的注释与词汇表本身是分开的，并且涉及不同类型的证据。这些记录都在Bioconductor的注释包中。

我们可以使用GO.db包来快速地访问GO词汇，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
library(GO.db)
GO.db # metadata
## GODb object:
## | GOSOURCENAME: Gene Ontology
## | GOSOURCEURL: ftp://ftp.geneontology.org/pub/go/godatabase/archive/latest-lite/
## | GOSOURCEDATE: 2017-Nov01
## | Db type: GODb
## | package: AnnotationDbi
## | DBSCHEMA: GO_DB
## | GOEGSOURCEDATE: 2017-Nov6
## | GOEGSOURCENAME: Entrez Gene
## | GOEGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
## | DBSCHEMAVERSION: 2.1
## 
## Please see: help('select') for usage information

使用AnnotationDbi包中的keys，columns和select函数也很容易在地id与不同terms之间进行映射，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
k5 = keys(GO.db)[1:5]
cgo = columns(GO.db)
select(GO.db, keys=k5, columns=cgo[1:3])
## 'select()' returned 1:1 mapping between keys and columns
##         GOID
## 1 GO:0000001
## 2 GO:0000002
## 3 GO:0000003
## 4 GO:0000006
## 5 GO:0000007
##                                                                                                                                                                                                                                                                                      DEFINITION
## 1                                                                                                       The distribution of mitochondria, including the mitochondrial genome, into daughter cells after mitosis or meiosis, mediated by interactions between mitochondria and the cytoskeleton.
## 2                                                                                                                                             The maintenance of the structure and integrity of the mitochondrial genome; includes replication and segregation of the mitochondrial chromosome.
## 3                                                                                                                                                                  The production of new individuals that contain some portion of genetic material inherited from one or more parent organisms.
## 4                                      Enables the transfer of zinc ions (Zn2+) from one side of a membrane to the other, probably powered by proton motive force. In high-affinity transport the transporter is able to bind the solute even if it is only present at very low concentrations.
## 5 Enables the transfer of a solute or solutes from one side of a membrane to the other according to the reaction: Zn2+ = Zn2+, probably powered by proton motive force. In low-affinity transport the transporter is able to bind the solute only if it is present at very high concentrations.
##   ONTOLOGY
## 1       BP
## 2       BP
## 3       BP
## 4       MF
## 5       MF

词汇表的图形结构被编码在SQLite数据库的表中。我们可以使用RSQLite接口对此进行查询，如下所示：

1
2
3
4
5
6
7
con = GO_dbconn()
dbListTables(con)
##  [1] "go_bp_offspring" "go_bp_parents"   "go_cc_offspring"
##  [4] "go_cc_parents"   "go_mf_offspring" "go_mf_parents"  
##  [7] "go_obsolete"     "go_ontology"     "go_synonym"     
## [10] "go_term"         "map_counts"      "map_metadata"   
## [13] "metadata"        "sqlite_stat1"

以下查询提示了一些内部标识符：

1
2
3
4
5
6
7
dbGetQuery(con, "select _id, go_id, term from go_term limit 5")
##   _id      go_id                                        term
## 1  30 GO:0000001                   mitochondrion inheritance
## 2  32 GO:0000002            mitochondrial genome maintenance
## 3  33 GO:0000003                                reproduction
## 4  37 GO:0042254                         ribosome biogenesis
## 5  38 GO:0044183 protein binding involved in protein folding

我们可以将 mitochondrion inheritance term追溯到父项和祖父母项，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
dbGetQuery(con, "select * from go_bp_parents where _id=30")
##   _id _parent_id relationship_type
## 1  30      26537              is_a
## 2  30      26540              is_a
dbGetQuery(con, "select _id, go_id, term from go_term where _id=26616")
##     _id      go_id
## 1 26616 GO:0048387
##                                                              term
## 1 negative regulation of retinoic acid receptor signaling pathway
dbGetQuery(con, "select * from go_bp_parents where _id=26616")
##     _id _parent_id    relationship_type
## 1 26616       8389                 is_a
## 2 26616      26614                 is_a
## 3 26616      26613 negatively_regulates
dbGetQuery(con, "select _id, go_id, term from go_term where _id=5932")
##    _id      go_id                    term
## 1 5932 GO:0019237 centromeric DNA binding

将 “mitochondrion inheritance” 视为过程“mitochondrion distribution”和 “organelle inheritance”在概念上的精练是有意义的，这两个term在数据库中被为父项。

可以使用 GO_dbschema()来查看整个数据库模式。

KEGG

自Bioconductor诞生以来，KEGG的注释就能在Bioconductor中人使用了，但KEGG的数据库使用权限已经进行了更改。当我们使用KEGG.db加载后会出现以下信息，如下所示：

1
2
3
4
5
6
7
> library(KEGG.db)
KEGG.db contains mappings based on older data because the original
  resource was removed from the the public domain before the most
  recent update was produced. This package should now be considered
  deprecated and future versions of Bioconductor may not have it
  available.  Users who want more current data are encouraged to look
  at the KEGGREST or reactome.db packages

因此我们可以关注KEGGREST这个包，它需要联网。这是一个非常有用的，基于Entrez标识符的工具。现在我们查询一下BRCA2的信息（它的EntrezID为675），如下所示：

1
2
3
4
5
6
library(KEGGREST)
brca2K = keggGet("hsa:675")
names(brca2K[[1]])
##  [1] "ENTRY"      "NAME"       "DEFINITION" "ORTHOLOGY"  "ORGANISM"  
##  [6] "PATHWAY"    "DISEASE"    "BRITE"      "POSITION"   "MOTIF"     
## [11] "DBLINKS"    "STRUCTURE"  "AASEQ"      "NTSEQ"

我们也可以通过keggGet函数来获取构成通路模式的基因列表，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
brpat = keggGet("path:hsa05212")
names(brpat[[1]])
##  [1] "ENTRY"       "NAME"        "DESCRIPTION" "CLASS"       "PATHWAY_MAP"
##  [6] "DISEASE"     "DRUG"        "ORGANISM"    "GENE"        "COMPOUND"   
## [11] "KO_PATHWAY"  "REFERENCE"
brpat[[1]]$GENE[seq(1,132,2)] # entrez gene ids
##  [1] "3845"  "5290"  "5293"  "5291"  "5295"  "5296"  "8503"  "9459" 
##  [9] "5879"  "5880"  "5881"  "4790"  "5970"  "207"   "208"   "10000"
## [17] "1147"  "3551"  "8517"  "572"   "598"   "842"   "369"   "673"  
## [25] "5894"  "5604"  "5594"  "5595"  "5599"  "5602"  "5601"  "5900" 
## [33] "5898"  "5899"  "10928" "998"   "7039"  "1950"  "1956"  "2064" 
## [41] "2475"  "6198"  "6199"  "3716"  "6774"  "6772"  "7422"  "1029" 
## [49] "1019"  "1021"  "595"   "5925"  "1869"  "1870"  "1871"  "7157" 
## [57] "1026"  "1647"  "4616"  "10912" "581"   "578"   "1643"  "51426"
## [65] "7040"  "7042"

KEGGREST还有许多值得研究的地方，例如还可以查询BRCA2（人类）关于胰腺癌途径的静态图像，如下所示：

1
2
3
4
library(png)
library(grid)
brpng = keggGet("hsa05212", "image")
grid.raster(brpng)

plot of chunk getp

其它本体

rols包含有与EMBL-EBI连接的接口 Ontology Lookup Service.

1
2
3
4
5
6
7
8
9
10
11
12
library(rols)
oo = Ontologies()
oo
## Object of class 'Ontologies' with 198 entries
##    GENEPIO, MP ... SEPIO, SIBO
oo[[1]]
## Ontology: Genomic Epidemiology Ontology (genepio)  
##   The Genomic Epidemiology Ontology (GenEpiO) covers vocabulary
##   necessary to identify, document and research foodborne pathogens
##   and associated outbreaks.
##    Loaded: 2017-04-10 Updated: 2017-10-20 Version: 2017-04-09 
##    4351 terms  137 properties  38 individuals

为了控制查询检索中涉及的网络流量，搜索分为几个阶段，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
glis = OlsSearch("glioblastoma")
glis
## Object of class 'OlsSearch':
##   query: glioblastoma 
##   requested: 20 (out of 502)
##   response(s): 0
res = olsSearch(glis)
dim(res)
## NULL
resdf = as(res, "data.frame") # get content
resdf[1:4,1:4]
##                                                     id
## 1 ncit:class:http://purl.obolibrary.org/obo/NCIT_C3058
## 2     omit:http://purl.obolibrary.org/obo/OMIT_0007102
## 3    ordo:class:http://www.orpha.net/ORDO/Orphanet_360
## 4   hp:class:http://purl.obolibrary.org/obo/HP_0100843
##                                           iri   short_form        label
## 1   http://purl.obolibrary.org/obo/NCIT_C3058   NCIT_C3058 Glioblastoma
## 2 http://purl.obolibrary.org/obo/OMIT_0007102 OMIT_0007102 Glioblastoma
## 3      http://www.orpha.net/ORDO/Orphanet_360 Orphanet_360 Glioblastoma
## 4   http://purl.obolibrary.org/obo/HP_0100843   HP_0100843 Glioblastoma
resdf[1,5]  # full description for one instance
## [[1]]
## [1] "The most malignant astrocytic tumor (WHO grade IV).  It is composed of poorly differentiated neoplastic astrocytes and it is characterized by the presence of cellular polymorphism, nuclear atypia, brisk mitotic activity, vascular thrombosis, microvascular proliferation and necrosis. It typically affects adults and is preferentially located in the cerebral hemispheres. It may develop from diffuse astrocytoma WHO grade II or anaplastic astrocytoma (secondary glioblastoma, IDH-mutant), but more frequently, it manifests after a short clinical history de novo, without evidence of a less malignant precursor lesion (primary glioblastoma, IDH- wildtype). (Adapted from WHO)"

ontologyIndex 包支持导入开放生物本体(OBO, Open Biological Ontologies)格式的数据，并含有用于查询和可视化本体系统高效的工具。

通用基因集管理

GSEABase 包有一个用于管理基因集和集合的优秀工具。我们可以从MSigDb中导入胶质母细胞瘤相关的基因集来说明一下，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
library(GSEABase)
glioG = getGmt(system.file("gmt/glioSets.gmt", package="ph525x"))
## Warning in readLines(con, ...): incomplete final line found on '/
## Library/Frameworks/R.framework/Versions/3.4/Resources/library/ph525x/gmt/
## glioSets.gmt'
glioG
## GeneSetCollection
##   names: BALDWIN_PRKCI_TARGETS_UP, BEIER_GLIOMA_STEM_CELL_DN, ..., ZHENG_GLIOBLASTOMA_PLASTICITY_UP (47 total)
##   unique identifiers: ADA, AQP9, ..., ZFP28 (3671 total)
##   types in collection:
##     geneIdType: NullIdentifier (1 total)
##     collectionType: NullCollection (1 total)
head(geneIds(glioG[[1]]))
## [1] "ADA"      "AQP9"     "ATP2B4"   "ATP6V1G1" "CBX6"     "CCDC165"

模式生物的统一，自我描述方法

OrganismDb包简化了对注释的访问。还可以针对TxDb和org.[Nn].eg.db进行直接查询，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
library(Homo.sapiens)
class(Homo.sapiens)
## [1] "OrganismDb"
## attr(,"package")
## [1] "OrganismDbi"
Homo.sapiens
## OrganismDb Object:
## # Includes GODb Object:  GO.db 
## # With data about:  Gene Ontology 
## # Includes OrgDb Object:  org.Hs.eg.db 
## # Gene data about:  Homo sapiens 
## # Taxonomy Id:  9606 
## # Includes TxDb Object:  TxDb.Hsapiens.UCSC.hg19.knownGene 
## # Transcriptome data about:  Homo sapiens 
## # Based on genome:  hg19 
## # The OrgDb gene id ENTREZID is mapped to the TxDb gene id GENEID .
tx = transcripts(Homo.sapiens)
## 'select()' returned 1:1 mapping between keys and columns
keytypes(Homo.sapiens)
##  [1] "ACCNUM"       "ALIAS"        "CDSID"        "CDSNAME"     
##  [5] "DEFINITION"   "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
##  [9] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL" 
## [13] "EXONID"       "EXONNAME"     "GENEID"       "GENENAME"    
## [17] "GO"           "GOALL"        "GOID"         "IPI"         
## [21] "MAP"          "OMIM"         "ONTOLOGY"     "ONTOLOGYALL" 
## [25] "PATH"         "PFAM"         "PMID"         "PROSITE"     
## [29] "REFSEQ"       "SYMBOL"       "TERM"         "TXID"        
## [33] "TXNAME"       "UCSCKG"       "UNIGENE"      "UNIPROT"
columns(Homo.sapiens)
##  [1] "ACCNUM"       "ALIAS"        "CDSCHROM"     "CDSEND"      
##  [5] "CDSID"        "CDSNAME"      "CDSSTART"     "CDSSTRAND"   
##  [9] "DEFINITION"   "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
## [13] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL" 
## [17] "EXONCHROM"    "EXONEND"      "EXONID"       "EXONNAME"    
## [21] "EXONRANK"     "EXONSTART"    "EXONSTRAND"   "GENEID"      
## [25] "GENENAME"     "GO"           "GOALL"        "GOID"        
## [29] "IPI"          "MAP"          "OMIM"         "ONTOLOGY"    
## [33] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
## [37] "PROSITE"      "REFSEQ"       "SYMBOL"       "TERM"        
## [41] "TXCHROM"      "TXEND"        "TXID"         "TXNAME"      
## [45] "TXSTART"      "TXSTRAND"     "TXTYPE"       "UCSCKG"      
## [49] "UNIGENE"      "UNIPROT"

面向平台的注释

通过在NCBI GEO的GPL信息页面 上对信息进行排序，我们就可以看到最常用的寡核苷阵列平台（数据库中有4760个系列）就是Affy Human Genome U133 plus 2.0 array (GPL 570)。我们可以使用hgu133plus2.db对这些数据进行注释，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
library(hgu133plus2.db)
## 
hgu133plus2.db
## ChipDb object:
## | DBSCHEMAVERSION: 2.1
## | Db type: ChipDb
## | Supporting package: AnnotationDbi
## | DBSCHEMA: HUMANCHIP_DB
## | ORGANISM: Homo sapiens
## | SPECIES: Human
## | MANUFACTURER: Affymetrix
## | CHIPNAME: Human Genome U133 Plus 2.0 Array
## | MANUFACTURERURL: http://www.affymetrix.com/support/technical/byproduct.affx?product=hg-u133-plus
## | EGSOURCEDATE: 2015-Sep27
## | EGSOURCENAME: Entrez Gene
## | EGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
## | CENTRALID: ENTREZID
## | TAXID: 9606
## | GOSOURCENAME: Gene Ontology
## | GOSOURCEURL: ftp://ftp.geneontology.org/pub/go/godatabase/archive/latest-lite/
## | GOSOURCEDATE: 20150919
## | GOEGSOURCEDATE: 2015-Sep27
## | GOEGSOURCENAME: Entrez Gene
## | GOEGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
## | KEGGSOURCENAME: KEGG GENOME
## | KEGGSOURCEURL: ftp://ftp.genome.jp/pub/kegg/genomes
## | KEGGSOURCEDATE: 2011-Mar15
## | GPSOURCENAME: UCSC Genome Bioinformatics (Homo sapiens)
## | GPSOURCEURL: ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19
## | GPSOURCEDATE: 2010-Mar22
## | ENSOURCEDATE: 2015-Jul16
## | ENSOURCENAME: Ensembl
## | ENSOURCEURL: ftp://ftp.ensembl.org/pub/current_fasta
## | UPSOURCENAME: Uniprot
## | UPSOURCEURL: http://www.uniprot.org/
## | UPSOURCEDATE: Thu Oct  1 23:31:58 2015
## 
## Please see: help('select') for usage information

这个资源（以及ChipDb类的所有实例）的基本目的是在探针集(probeset)识别符和更高层次的基因组注释之间进行映射。

有关探针的详细信息（探针集的组成部分）已经由那些后缀为probe的文件包提供，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
library(hgu133plus2probe)
head(hgu133plus2probe)
##                    sequence    x   y Probe.Set.Name
## 1 CACCCAGCTGGTCCTGTGGATGGGA  718 317      1007_s_at
## 2 GCCCCACTGGACAACACTGATTCCT 1105 483      1007_s_at
## 3 TGGACCCCACTGGCTGAGAATCTGG  584 901      1007_s_at
## 4 AAATGTTTCCTTGTGCCTGCTCCTG  192 205      1007_s_at
## 5 TCCTTGTGCCTGCTCCTGTACTTGT  844 979      1007_s_at
## 6 TGCCTGCTCCTGTACTTGTCCTCAG  537 971      1007_s_at
##   Probe.Interrogation.Position Target.Strandedness
## 1                         3330           Antisense
## 2                         3443           Antisense
## 3                         3512           Antisense
## 4                         3563           Antisense
## 5                         3570           Antisense
## 6                         3576           Antisense
dim(hgu133plus2probe)
## [1] 604258      6

将探针集标识符映射到基因水平的信息可以提示一些有意思的歧视，如下所示：

1
2
3
4
5
6
7
8
9
select(hgu133plus2.db, keytype="PROBEID", 
  columns=c("SYMBOL", "GENENAME", "PATH", "MAP"), keys="1007_s_at")
## 'select()' returned 1:many mapping between keys and columns
##     PROBEID  SYMBOL                                    GENENAME PATH
## 1 1007_s_at    DDR1 discoidin domain receptor tyrosine kinase 1 <NA>
## 2 1007_s_at MIR4640                               microRNA 4640 <NA>
##       MAP
## 1 6p21.33
## 2 6p21.33

显然，该探针集合可以用于mRNA和miRNA丰度的定量。作为稳定的检查，我们可以看到，不同的符号映射到了相同的细胞带（最后一句不懂，原文为： As a sanity check we see that the distinct symbols map to the same cytoband）。

总结

我们现在已经拥有了含有从核酸到通路水平的许多数据。通过Bioconductor.org上的View就可以查看现有的一些资源。

参考资料

  1. Genomic annotation in Bioconductor: The general situation

发表于 | 分类于
| 字数统计: 5,164 | 阅读时长 ≈ 24

前言

原始的Rmd文档可以参考Github上。

IRanges 和 GRanges 对象是Bioconductor基础数据类型的核心成员，其中 IRanges 定义了 integer ranges ，它用于解决基因组上的位置问题，而 GRanges 则包含了染色体和DNA链的信息。这里我们简介绍一下这两个对象，以及如何操作 IRanges 和 GRanges 。

IRanges

第一步是加载IRanges包，如下所示：

1
library(IRanges)

IRanges 函数定义了区间范围(interval ranges)。如果你输入2个数据，那么就表示一个区间，例如 [5,10]={5,6,7,8,9,10}，它的宽度是6。另外，如果我们要查看一个 IRanges 定义的区间宽度，我们使用的width()函数，而非 length()函数，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
ir <- IRanges(5,10)
ir
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         5        10         6
start(ir)
## [1] 5
end(ir)
## [1] 10
width(ir)
## [1] 6
# for detailed information on the IRanges class:
# ?IRanges
length(ir)
# [1] 1

一个单独的IRanges对象可以含有多个范围。现在我们可以在一个向量中指定start，在另外一个向量中指定end，如下所示：

1
2
3
4
5
6
7
IRanges(start=c(3,5,17), end=c(10,8,20))
## IRanges object with 3 ranges and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         3        10         8
##   [2]         5         8         4
##   [3]        17        20         4

内部范围(Intra-range)操作

现在我们继续使用单个范围[5,10]，我们查看内部范围(intra-range){…,−2,−1,0,1,2,…}的一个数，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
# full details on the intra-range methods:
# ?"intra-range-methods"
ir
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         5        10         6
shift(ir, -2)
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         3         8         6

这里我们展示了应用于ir的多个不同操作的结果，结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
ir
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         5        10         6
narrow(ir, start=2)
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         6        10         5
narrow(ir, end=5)
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         5         9         5
flank(ir, width=3, start=TRUE, both=FALSE)
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         2         4         3
flank(ir, width=3, start=FALSE, both=FALSE)
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]        11        13         3
flank(ir, width=3, start=TRUE, both=TRUE)
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         2         7         6
ir * 2
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         6         8         3
ir * -2
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         2        13        12
ir + 2
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         3        12        10
ir - 2
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         7         8         2
resize(ir, 1)
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         5         5         1

我们使用下图展示了上面操作的结果。红色柱状部分显示了最初的ir范围。掌握这些操作的最好方法是在自己定义的范围内在控制台中亲自操作一下。

plot of chunk unnamed-chunk-6

内部范围方法

针对intra-range有着一系列的方法。还有一些方法是针对一系列范围的操作，它的输出取决于所有的范围，因此这些方法与intra-range方法不同，后者不会改变输出。现在我们一些案例来说明一下。range()函数提供了从最左侧到最右侧的整数输出，如下所示：

注：原文我不太理解，翻译的不精准，原文如下：

There are also a set of inter-range methods. These are
operations which work on a set of ranges, and the output depends on all
the ranges, thus distinguishes these methods from the intra-range methods, for which the other ranges in the set do not change the output. This is best explained with some examples. The range function gives the integer range from the start of the leftmost range to the end of the rightmost range:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# full details on the inter-range methods:
# ?"inter-range-methods"
(ir <- IRanges(start=c(3,5,17), end=c(10,8,20)))
## IRanges object with 3 ranges and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         3        10         8
##   [2]         5         8         4
##   [3]        17        20         4
range(ir)
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         3        20        18

reduce 函数能折叠这些范围，以便于整个可以只被输出范围中的一个范围覆盖，如下所示：

1
2
3
4
5
6
reduce(ir)
## IRanges object with 2 ranges and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         3        10         8
##   [2]        17        20         4

注：ir原来有3个范围，其中[5,8]范围是在[3,10]之间，因此前者被折叠到了后者之中。

gaps 函数能够给出上面两个范围的区间，但是它不覆盖任何ir中的范围，如下所示：

1
2
3
4
5
gaps(ir)
## IRanges object with 1 range and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]        11        16         6

disjoin 会打断 ir的范围，并它们分成不同的范围，如下所示：

disjoin返回一个disjoint对象，它会查看参数中端点的并集，返回不相交的对象。换话句讲，结果中由区间的最大长度范围组成，在该范围上，对象x中的重叠范围集合是相同的，并且至少为1。

1
2
3
4
5
6
7
8
disjoin(ir)
## IRanges object with 4 ranges and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         3         4         2
##   [2]         5         8         4
##   [3]         9        10         2
##   [4]        17        20         4

现在我们使用示意图来说明一下：

第一二行是原始的范围，第四行是经过disjoin处理后的范围。

GRanges

GRanges对象包含了Iranges，它含有2个方面的重要信息：

  • 染色体信息（Bioconductor称为 seqnames）。
  • DNA链。

链的信息可以指定为 + 号或 - ，或者是保留未指定的 * 。正义链的特征则是在编号线上从左到右的生物学方面，而负链特征则是从右到左。就IRanges而言，正链特征是从start到 end ,负链则是从 end 到 start 。这初看起来比较混乱，但这是必需的，因此 width 被定义为 end - start + 1，并且不允许有负的宽度范围。因为DNA有两条链，这两条链具有相反的方向性，因此当我们说DNA时，可以只指明其中的一条链。

通过IRnage、染色体名称和一条链，我们可以确定我们与另外一个研究人员所指的DNA分子具有相同的区域和链，因为我们使用的相同的基因组构建的。GRanges对象中还可以包括其它部分的信息，不过上述的两个信息则是最重要的 。

1
library(GenomicRanges)

现在我们创建一个虚构染色体chrZ 上的两个范围。我们会说这两个范围指的基因组hg19上的范围。因为我们没有将我们的基因组链接到数据库，因此我们可以指定一个在hg19中并不存在的染色体，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
gr <- GRanges("chrZ", IRanges(start=c(5,10),end=c(35,45)),
              strand="+", seqlengths=c(chrZ=100L))
gr
## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrZ  [ 5, 35]      +
##   [2]     chrZ  [10, 45]      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome
genome(gr) <- "hg19"
gr
## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrZ  [ 5, 35]      +
##   [2]     chrZ  [10, 45]      +
##   -------
##   seqinfo: 1 sequence from hg19 genome

我们可以调用 seqnames 和 seqlengths 来查看上面的信息：

1
2
3
4
5
6
7
8
seqnames(gr)
## factor-Rle of length 2 with 1 run
##   Lengths:    2
##   Values : chrZ
## Levels(1): chrZ
seqlengths(gr)
## chrZ 
##  100

我们可以使用 shift 函数（就像在IRanges中使用的一样）来处理GRanges。但是，当我们尝试超过染色体长度之外的范围时，就会出现警告信息，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
shift(gr, 10)
## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrZ  [15, 45]      +
##   [2]     chrZ  [20, 55]      +
##   -------
##   seqinfo: 1 sequence from hg19 genome
shift(gr, 80)
## Warning in valid.GenomicRanges.seqinfo(x, suggest.trim = TRUE): GRanges object contains 2 out-of-bound ranges located on sequence
##   chrZ. Note that only ranges located on a non-circular sequence
##   whose length is not NA can be considered out-of-bound (use
##   seqlengths() and isCircular() to get the lengths and circularity
##   flags of the underlying sequences). You can use trim() to trim
##   these ranges. See ?`trim,GenomicRanges-method` for more
##   information.
## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrZ [85, 115]      +
##   [2]     chrZ [90, 125]      +
##   -------
##   seqinfo: 1 sequence from hg19 genome

如果我们使用 trim 处理这个范围，我们就能获取剩下的范围，而不考虑超出染色体长的部分：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
trim(shift(gr, 80))
## Warning in valid.GenomicRanges.seqinfo(x, suggest.trim = TRUE): GRanges object contains 2 out-of-bound ranges located on sequence
##   chrZ. Note that only ranges located on a non-circular sequence
##   whose length is not NA can be considered out-of-bound (use
##   seqlengths() and isCircular() to get the lengths and circularity
##   flags of the underlying sequences). You can use trim() to trim
##   these ranges. See ?`trim,GenomicRanges-method` for more
##   information.
## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrZ [85, 100]      +
##   [2]     chrZ [90, 100]      +
##   -------
##   seqinfo: 1 sequence from hg19 genome

我们可以使用 mcols 函数（链表示元数据列）向每个范围添加列信息。需要注意的是，这个操作对IRanges也是可行的。我们还可以通过指定 NULL 来移除列信息，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
mcols(gr)
## DataFrame with 2 rows and 0 columns
mcols(gr)$value <- c(-1,4)
gr
## GRanges object with 2 ranges and 1 metadata column:
##       seqnames    ranges strand |     value
##          <Rle> <IRanges>  <Rle> | <numeric>
##   [1]     chrZ  [ 5, 35]      + |        -1
##   [2]     chrZ  [10, 45]      + |         4
##   -------
##   seqinfo: 1 sequence from hg19 genome
mcols(gr)$value <- NULL
> gr
# GRanges object with 2 ranges and 0 metadata columns:
#      seqnames    ranges strand
#         <Rle> <IRanges>  <Rle>
#  [1]     chrZ      5-35      +
#  [2]     chrZ     10-45      +
#  -------
#  seqinfo: 1 sequence from hg19 genome

GRangesList

当我们涉及基因时，我们通过创建一个GRanges列表是很有用的，例如用于表示分组信息（比如每个基因的外显子）。该列表的元素是基因，并且在每个元素中，外显子的范围被定义为GRanges。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
gr2 <- GRanges("chrZ",IRanges(11:13,51:53))
grl <- GRangesList(gr, gr2)
grl
## GRangesList object of length 2:
## [[1]] 
## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrZ  [ 5, 35]      +
##   [2]     chrZ  [10, 45]      +
## 
## [[2]] 
## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames   ranges strand
##   [1]     chrZ [11, 51]      *
##   [2]     chrZ [12, 52]      *
##   [3]     chrZ [13, 53]      *
## 
## -------
## seqinfo: 1 sequence from hg19 genome

GRangesList 对象中的长度就是其中的GRanges 对象个数。为了获取每个GRanges的长度，我们可以调用 elementNROWS。我们可以通过两个方括号的典型列表索引方法来建立索引，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
length(grl)
## [1] 2
elementNROWS(grl)
## [1] 2 3
grl[[1]]
## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrZ  [ 5, 35]      +
##   [2]     chrZ  [10, 45]      +
##   -------
##   seqinfo: 1 sequence from hg19 genome

width函数会返回一个 integerList。如果使用sum，我们会得到列表中每个GRanges对象的宽度向量，如下所示：

1
2
3
4
5
6
width(grl)
## IntegerList of length 2
## [[1]] 31 36
## [[2]] 41 41 41
sum(width(grl))
## [1]  67 123

我们可以像以前那样添加元数据(metadata)，现在每个GRange对象都有一行元数据。当我们输出GRangesList时它们并不业显示出来，但我们可以通过 mcols 进行访问，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
mcols(grl)$value <- c(5,7)
grl
## GRangesList object of length 2:
## [[1]] 
## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrZ  [ 5, 35]      +
##   [2]     chrZ  [10, 45]      +
## 
## [[2]] 
## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames   ranges strand
##   [1]     chrZ [11, 51]      *
##   [2]     chrZ [12, 52]      *
##   [3]     chrZ [13, 53]      *
## 
## -------
## seqinfo: 1 sequence from hg19 genome
mcols(grl)
## DataFrame with 2 rows and 1 column
##       value
##   <numeric>
## 1         5
## 2         7

findOverlaps与%over%

现在我们来看两个比较GRanges对象的常用方法，首先我们要创建两组范围如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
(gr1 <- GRanges("chrZ",IRanges(c(1,11,21,31,41),width=5),strand="*"))
## GRanges object with 5 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrZ  [ 1,  5]      *
##   [2]     chrZ  [11, 15]      *
##   [3]     chrZ  [21, 25]      *
##   [4]     chrZ  [31, 35]      *
##   [5]     chrZ  [41, 45]      *
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
(gr2 <- GRanges("chrZ",IRanges(c(19,33),c(38,35)),strand="*"))
## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrZ  [19, 38]      *
##   [2]     chrZ  [33, 35]      *
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

findOverlaps 会返回一个 Hits 对象，这个对象含有关于检索范围的信息（第一个参数），这个范围与主题（第二个参数）有哪些地方重叠。还有许多参数用于指定计算哪种类型的重叠，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
fo <- findOverlaps(gr1, gr2)
fo
## Hits object with 3 hits and 0 metadata columns:
##       queryHits subjectHits
##       <integer>   <integer>
##   [1]         3           1
##   [2]         4           1
##   [3]         4           2
##   -------
##   queryLength: 5 / subjectLength: 2
queryHits(fo)
## [1] 3 4 4
subjectHits(fo)
## [1] 1 1 2

得到重叠信息的另外一种方式就是使用 %over% ，它返回一个逻辑向量，这个向量表示了第一个参数中的范围与第二个参数中任何范围重叠的信息，如下所示：

1
2
3
4
5
6
7
8
9
10
gr1 %over% gr2
## [1] FALSE FALSE  TRUE  TRUE FALSE
gr1[gr1 %over% gr2]
## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrZ  [21, 25]      *
##   [2]     chrZ  [31, 35]      *
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

需要注意的是，这两种操作方式都是链特异性的，虽然 findOverlaps 有一个 ignore.strand 选择参数，如下所示：

1
2
3
4
gr1 <- GRanges("chrZ",IRanges(1,10),strand="+")
gr2 <- GRanges("chrZ",IRanges(1,10),strand="-")
gr1 %over% gr2
## [1] FALSE

Rle 和Views

最后，我们来简单了解一下由IRanges定义的两个相关类，即Rle类和Views类。Rle全称为run-length encoding，它表示的是重复数据的压缩工，用于代替类似于[1, 1, 1, 1]这样的储存。

我们只用储存数据1，以及重复次数4即可。数据越重复，使用Rle对象的压缩效果就越明显。

我们使用 str 来查看 Rle的内部结构，用于显示它的储存数据和重复次数，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
(r <- Rle(c(1,1,1,0,0,-2,-2,-2,rep(-1,20))))
## numeric-Rle of length 28 with 4 runs
##   Lengths:  3  2  3 20
##   Values :  1  0 -2 -1
str(r)
## Formal class 'Rle' [package "S4Vectors"] with 4 slots
##   ..@ values         : num [1:4] 1 0 -2 -1
##   ..@ lengths        : int [1:4] 3 2 3 20
##   ..@ elementMetadata: NULL
##   ..@ metadata       : list()
as.numeric(r)
##  [1]  1  1  1  0  0 -2 -2 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
## [24] -1 -1 -1 -1 -1

Views对象可以被为一个“窗口”用于查看一个序列信息，如下所示：

A Views object can be thought of as “windows” looking into a sequence.

1
2
3
4
5
6
7
(v <- Views(r, start=c(4,2), end=c(7,6)))
## Views on a 28-length Rle subject
## 
## views:
##     start end width
## [1]     4   7     4 [ 0  0 -2 -2]
## [2]     2   6     5 [ 1  1  0  0 -2]

需要注意的是，Veiws对象的内部结构只是原始对象，以及指定窗口的IRangs。使用Views的最好好处就是，当原始对象没有储存在内存中时，在这种情况下，View对象就是一个轻量级的类，它能帮助我们引用子序列，而不必将整个序列都加载到内存中，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
str(v)
## Formal class 'RleViews' [package "IRanges"] with 5 slots
##   ..@ subject        :Formal class 'Rle' [package "S4Vectors"] with 4 slots
##   .. .. ..@ values         : num [1:4] 1 0 -2 -1
##   .. .. ..@ lengths        : int [1:4] 3 2 3 20
##   .. .. ..@ elementMetadata: NULL
##   .. .. ..@ metadata       : list()
##   ..@ ranges         :Formal class 'IRanges' [package "IRanges"] with 6 slots
##   .. .. ..@ start          : int [1:2] 4 2
##   .. .. ..@ width          : int [1:2] 4 5
##   .. .. ..@ NAMES          : NULL
##   .. .. ..@ elementType    : chr "integer"
##   .. .. ..@ elementMetadata: NULL
##   .. .. ..@ metadata       : list()
##   ..@ elementType    : chr "ANY"
##   ..@ elementMetadata: NULL
##   ..@ metadata       : list()

基因组元件的应用：链特异性操作

在这一部分，我们使用一个小型的范围来说明一下基因的范围内(intra-range)操作，包括reduce, disjoin与gaps。我们添加链和seqname信息，并显示如想调整大小和侧翼(flank)用于识别TSS和启动子区域，如下所示：

一组简单的范围

1
2
ir <- IRanges(c(3, 8, 14, 15, 19, 34, 40),
  width = c(12, 6, 6, 15, 6, 2, 7))

现在我们使用图片形式查看一下ir，以及内部范围(intra-range)操作，如下所示：

1
2
3
4
5
par(mfrow=c(4,1), mar=c(4,2,2,2))
plotRanges(ir, xlim=c(0,60))
plotRanges(reduce(ir), xlim=c(0,60))
plotRanges(disjoin(ir), xlim=c(0,60))
plotRanges(gaps(ir), xlim=c(0,60))

plot of chunk lkir

注：plotRange()函数似乎不在这个包中，网上找到了它的原代码，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
plotRanges <- function(x, xlim = x, main = deparse(substitute(x)),
                       col = "black", sep = 0.5, ...)
{
  height <- 1
  if (is(xlim, "Ranges"))
    xlim <- c(min(start(xlim)), max(end(xlim)))
  bins <- disjointBins(IRanges(start(x), end(x) + 1))
  plot.new()
  plot.window(xlim, c(0, max(bins)*(height + sep)))
  ybottom <- bins * (sep + height) - height
  rect(start(x)-0.5, ybottom, end(x)+0.5, ybottom + height, col = col, ...)
  title(main)
  axis(1)
}

reduce(x)函数生成一组不重叠的范围，它覆盖了x中的所有位点。这种操作用于生了你许多转录本的基因模式的复杂性，其中我们可能只想要转录出已知区间的地址，不考虑transcript of residence。

disjoin(x)生成一组覆盖x的所有位置的范围，使得分离输出中的任何范围都不会与x中的任何间隔的端点重叠。这为我们提供了连续间隔最大可能的集合，这些连续间隔在原始间隔集中具有端点的地方被分割开。

gaps(x)产生一组覆盖[start(X)，end(X)]中未被x中任何范围覆盖的位置的范围。它能给出编码序列位置和外显子间隔，这可用于枚举内含子。

GRanges的扩展

现在我们添加上染色体与链信息，如下所示： 

1
2
library(GenomicRanges)
gir = GRanges(seqnames="chr1", ir, strand=c(rep("+", 4), rep("-",3)))

让我们假设间隔代表基因。下图说明了转录起始位点(绿色)、上游启动子区域(紫色)、下游启动子区域(棕色)的位置，如下所示：

1
2
3
4
5
6
par(mfrow=c(4,1), mar=c(4,2,2,2))
plotGRanges(gir, xlim=c(0,60))
# 
plotGRanges(resize(gir,1), xlim=c(0,60),col="green")
plotGRanges(flank(gir,3), xlim=c(0,60), col="purple")
plotGRanges(flank(gir,2,start=FALSE), xlim=c(0,60), col="brown")

plot of chunk dopr

注：plotGRanges的源代码如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
plotGRanges = function (x, xlim = x, col = "black", sep = 0.5, xlimits = c(0, 
    60), ...) 
{
    main = deparse(substitute(x))
    ch = as.character(seqnames(x)[1])
    x = ranges(x)
    height <- 1
    if (is(xlim, "Ranges")) 
        xlim <- c(min(start(xlim)), max(end(xlim)))
    bins <- disjointBins(IRanges(start(x), end(x) + 1))
    plot.new()
    plot.window(xlim = xlimits, c(0, max(bins) * (height + sep)))
    ybottom <- bins * (sep + height) - height
    rect(start(x) - 0.5, ybottom, end(x) + 0.5, ybottom + height, 
        col = col, ...)
    title(main, xlab = ch)
    axis(1)
}

请注意，我们不需要采取特殊步骤来处理链中的差异

甲基化微阵列数据的可视化

在我们讨论 SummarizedExperiment applications时，我们输入了由Illumina 450k甲基化微阵列生成的数据。

在这一部分中，我们将介绍如何使用GenomicRanges和Gviz来研究一下在基因水平注释下的甲基化模型。这个思路非常简单：只需从所有样本中提取M值（甲基化与总DNA 比值的基因座特定的估计的log转换），并在选定基因的基因模型中绘制出来，如下所示：

我们回顾一下数据是如何获取和导入的：

1
2
3
4
5
6
7
8
9
library(ArrayExpress)
if (!file.exists("E-MTAB-5797.sdrf.txt")) nano = getAE("E-MTAB-5797")
library(minfi)
pref = unique(substr(dir(patt="idat"),1,17)) # find the prefix strings
raw = read.metharray(pref)
glioMeth = preprocessQuantile(raw) # generate SummarizedExperiment
## [preprocessQuantile] Mapping to genome.
## [preprocessQuantile] Fixing outliers.
## [preprocessQuantile] Quantile normalizing.

这些步骤需要网联，并且需要花几分钟。

一旦我们的会话中有 glioMeth ，需要将以下代码添加到会话中。我们将会介绍它是如何工作的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
MbyGene = function(mset, symbol="TP53", rad=5000) {
# phase 1: annotated GRanges for the gene
 require(erma)
 require(Gviz)
 gmod = suppressMessages(genemodel(symbol))     # erma utility
 gseq = as.character(seqnames(gmod)[1])
 gmod$transcript = symbol
# phase 2: filter down to the region of interest
 mlim = mset[which(seqnames(mset)==gseq),] # restrict to chromosome
 # now focus the methylation data to vicinity of gene
 d1 = subsetByOverlaps(GRanges(rowRanges(mlim),,, getM(mlim)), 
         range(gmod)+rad)
# phase 3: use Gviz
 plotTracks(list(DataTrack(d1), 
              GeneRegionTrack(gmod, 
               transcriptAnnotation="transcript", name=gseq), 
              GenomeAxisTrack(name=gseq, showTitle=TRUE)))
}

代码表明了三个阶段：获取基因区域，并添加转录本(transcript)注释用于绘制所有外显子的结合；将GenomicRatioSet（继承自

对代码的注释表明了三个阶段：获取基因区并添加用于所有外显子联合的信息性绘图的转录注释；将GenomicRatioSet（继承自RangedSummarizedExperiment）缩小到感兴趣的间隔，这个间隔是由基因模式和rad参数确定；使用Gviz来构建可用于绘图的对象，并绘制出来。

Gviz的详细信息在该包的文档中。现在我们绘制出以基因为中心的图形，如下所示：

1
MbyGene(glioMeth, symbol="TERT")

从上图中我们可以发现这些信息：

  • 基因组背景（以兆(megabase)为单位的染色体和区域）。
  • 感兴趣的基因在单行中的结构，它们表示了表达间隔的联合。
  • 450k针近的位置（蓝色数据点的x坐标）。
  • 甲基化估计中的样本间变异。
  • 所选基因组区域的甲基化变异。

在6x模块中，我们将学习如何使用额外的软件包来创建这种类型的交互式展示，允许我们选择基因并随意放大感兴趣的区域。

脚注

更多的信息可以查看 GenomicRanges 包，以及查阅PLOS Comp Bio上的文献，它是GenomicRanges的作者写的：

http://www.ploscompbiol.org/article/info%3Adoi%2F10.1371%2Fjournal.pcbi.1003118

使用vignettes也能查看这个包的众多功能文档，如下所示：

1
browseVignettes("GenomicRanges")

或者是使用help函数：

1
help(package="GenomicRanges", help_type="html")

对于Bed Tools的用户来说，HelloRanges包对于在BEd和GRanges概念框架之间转换概念非常有用。

参考资料

  1. bioc1_grangeOps.Rmd
  2. IRanges and GRanges