前言

  说到富集分析，大家肯定第一时间会想到GO、KEGG、GSEA等常见的基因富集分析。那什么是基因富集分析呢？简单来说就是把我们筛选出来的基因按照先验的知识分一个类，以便看看哪些基因的功能和我们的研究相关，然后再做个检验，根据P值的显著性来判断这个事件的发生是否为偶然事件。对于基因富集，我们都很熟悉就是看看差异是否富集到特定的GO条目或者特定的Pathway通路中，以便找到想要关注的基因。再或者通过GSEA找到表型与特征相关。
  今天我们想说的是genomic region富集，就拿ChIP-seq为例来说，我们都知道该技术用来研究转录因子或组蛋白在基因组上的结合位点，同时基因组按功能可以划分为很多特征区域如intron、centromere、ncRNA、pseudogene、snoRNA、telomeric_repeat等，如果我们想知道ChIP-seq的结果在基因组哪些特征区域中富集应该怎么办呢？此时，我们可以使用LOLA包来做富集分析，该包基于fisher精确检验来校验genomic region集在哪些特征区域富集。

代码

下面以R包自带的测试数据演示代码，具体如下：

source("http://bioconductor.org/biocLite.R")
biocLite("LOLA")

library("LOLA")
# load region database
dbPath = system.file("extdata", "hg19", package="LOLA")
regionDB = loadRegionDB(dbPath)

# load userSets
data("sample_input", package="LOLA") 
# load userUniverse
data("sample_universe", package="LOLA") 
head(userSets)
GRangesList object of length 2:
$setA
GRanges object with 3142 ranges and 0 metadata columns:
         seqnames               ranges strand
            <Rle>            <IRanges>  <Rle>
     [1]     chr1   [ 437151,  438164]      *
     [2]     chr1   [ 875730,  878363]      *
     [3]     chr1   [ 933387,  937410]      *
     [4]     chr1   [ 967966,  970238]      *
     [5]     chr1   [1016863, 1017439]      *
     ...      ...                  ...    ...
  [3138]     chrY [ 9364545,  9364859]      *
  [3139]     chrY [ 9385471,  9385777]      *
  [3140]     chrY [14532115, 14533600]      *
  [3141]     chrY [23696580, 23696878]      *
  [3142]     chrY [26959489, 26959716]      *

...
<1 more element>
-------
seqinfo: 69 sequences from an unspecified genome; no seqlengths
head(userUniverse)
GRanges object with 6 ranges and 0 metadata columns:
      seqnames           ranges strand
         <Rle>        <IRanges>  <Rle>
  [1]     chr1 [ 28735,  29810]      *
  [2]     chr1 [135124, 135563]      *
  [3]     chr1 [327790, 328229]      *
  [4]     chr1 [437151, 438164]      *
  [5]     chr1 [449273, 450544]      *
  [6]     chr1 [533219, 534114]      *
  -------
  seqinfo: 69 sequences from an unspecified genome; no seqlengths

# run the analysis
locResults = runLOLA(userSets, userUniverse, regionDB, cores=1)
head(locResults)
   userSet dbSet   collection   pValueLog logOddsRatio support rnkPV rnkLO
1:    setB     2 ucsc_example 264.4863407    7.7578283     850     1     1
2:    setA     2 ucsc_example 254.6188080    8.6487312     632     1     1
3:    setB     1 ucsc_example  34.6073689    3.3494078    5747     2     2
4:    setA     4 ucsc_example   1.7169689    1.2377725     124     2     2
5:    setA     5 ucsc_example   1.7169689    1.2377725     124     2     2
6:    setA     3 ucsc_example   0.1877354    0.9135696       8     4     4
   rnkSup maxRnk meanRnk     b    c     d                      description
1:      2      2    1.33   452 4981 20546                     ucsc_example
2:      2      2    1.33   670 2510 23017                     ucsc_example
3:      1      2    1.67 20018   84   980 CpG islands from UCSC annotation
4:      3      3    2.33   761 3018 22926                     ucsc_example
5:      3      3    2.33   761 3018 22926                     ucsc_example
6:      5      5    4.33    66 3134 23621                     ucsc_example
   cellType tissue antibody treatment dataSource                    filename
1:     <NA>   <NA>     <NA>      <NA>       <NA>             laminB1Lads.bed
2:     <NA>   <NA>     <NA>      <NA>       <NA>             laminB1Lads.bed
3:     <NA>   <NA>     <NA>      <NA>       <NA>            cpgIslandExt.bed
4:     <NA>   <NA>     <NA>      <NA>       <NA>          vistaEnhancers.bed
5:     <NA>   <NA>     <NA>      <NA>       <NA> vistaEnhancers_colNames.bed
6:     <NA>   <NA>     <NA>      <NA>       <NA>          numtSAssembled.bed
          qValue  size
1: 3.263317e-264  1302
2: 1.202713e-254  1302
3:  8.232086e-35 28691
4:  3.837612e-02  1339
5:  3.837612e-02  1340
6:  1.000000e+00    78

# Write out results
writeCombinedEnrichment(locResults, outFolder= "lolaResults")

可以看到这个R包用起来还是很方便的，关于用法就不再多说了。用的时候最主要分清userSets, userUniverse, regionDB分别代表什么就可以了，其中关于regionDB软件自带一些特征区域库，但是不一定是我们感兴趣的区域，所以我们需要自己构建数据库来使用。下面介绍一下构建自定义数据库的过程。
自定义数据库目录结构如下：

hg19
└── ucsc_example
    ├── collection.txt
    ├── index.txt
    ├── regions
    │   ├── cpgIslandExt.bed
    │   ├── laminB1Lads.bed
    │   ├── numtSAssembled.bed
    │   ├── vistaEnhancers.bed
    │   └── vistaEnhancers_colNames.bed
    └── sizes.txt

hg19就是regionDB的最上级目录，在该目录下可以有很多个平行的collection目录(例如这里的ucsc_example目录)，ucsc_example就是region相关的目录，包括以下几个方面：

1. regions：子目录，用于存放感兴趣区域的bed文件，目录结构如上所示；
2. collection.txt file：区域集的描述信息，内容格式如下：

collector       date    source  description
Nathan  2015-04-01      UCSC Genome Browser     Selected bed tracks from the UCSC database.

3. index.txt：区域bed文件的描述信息，内容格式如下：

filename        description
cpgIslandExt.bed        CpG islands from UCSC annotation

4. 脚本或描述如何制作区域集的文件(可选项)，例如这里的size.txt描述了每个bed文件包含的区间数量，内容如下：

filename        size
cpgIslandExt.bed        28691
laminB1Lads.bed 1302
numtSAssembled.bed      78
vistaEnhancers.bed      1339
vistaEnhancers_colNames.bed     1340

可以看到准备自定义的regionDB库也是相当的简单，只需按照要求准备好相应的文件即可，构建好感兴趣的区域集后就可以愉快地进行富集分析了。。。

最后

  由于该R包使用起来很方便快捷，就不再过多赘述了，最后附上该包的链接：http://www.bioconductor.org/packages/release/bioc/html/LOLA.html。。。。(ô‿ô)。。。