背景
GREAT(Genomic Regions Enrichment of Annotations Tool)是一种广泛使用的基因组区域功能富集工具，该工具于2010年由斯坦福大学开发。然而，作为在线工具，其存在注释数据过时、支持的物种和功能基因集数量少以及用户不可扩展等局限性。因此，有人就开发了一个本地实现GREAT算法的R包rGREAT，其默认支持600多个物种和大量的功能基因集，同时也支持用户自备基因集和物种基因组。此外，该包还实现了处理背景区域的通用方法。
基因组学和表观组学研究通常会生成许多感兴趣的基因组区域列表，例如来自全基因组或外显子组测序数据的单核苷酸变异(SNV)、ChIP-seq 的特定染色质修饰的peak或基因组甲基化测序的差异甲基化区域 (DMR)。下一步的分析自然是将生物学功能与这些基因组区域联系起来。一种广泛使用的方法是首先将基因组区域注释到最近的基因，然后针对特定功能的基因集做富集分析。
GREAT在做基因组区域富集时，考虑了基因在基因组上的位置分布与长度，采用了不同的策略，如上图所示。对于给定功能基因集中的基因，首先，将基因的TSS上下游分别延伸5kb和1kb来建立一个基础结构域；然后，将基础结构域的上下游再继续延申至最大1mb，或者到达它邻近基因的基础结构域，如此每个基因相当于被转化成了一个区间；最后，将这些转化的区间进行合并形成一个没有overlap的区间集，相当于将特定生物学功能相关的基因集转化为了“功能区间集”，然后使用二项分布来计算输入区域集是否在功能区间集中富集。
分析
- 使用R包里基因集和注释数据
一般常见物种如人、小鼠等，在R里面都可以找到相应基因注释数据库，如TxDb.Hsapiens.UCSC.hg38.knownGene。功能基因集也可直接使用R数据包里面的集合，如GO.db、msigdbr，使用格式为GO:BP、GO:CC、GO:MP、msigdb:H等。
library(rGREAT)
library(ChIPseeker)
peak <- readPeakFile('GSM1233959_peaks.narrowPeak')
res <- great(peak, "GO:BP", "txdb:hg38")
tab <- getEnrichmentTable(res)
head(tab)
id description
1 GO:1904464 regulation of matrix metallopeptidase secretion
2 GO:1904465 negative regulation of matrix metallopeptidase secretion
3 GO:1990773 matrix metallopeptidase secretion
4 GO:0043615 astrocyte cell migration
5 GO:2000405 negative regulation of T cell migration
6 GO:2000321 positive regulation of T-helper 17 cell differentiation
genome_fraction observed_region_hits fold_enrichment p_value p_adjust
1 0.0002162107 73 5.103134 0 0
2 0.0002162107 73 5.103134 0 0
3 0.0002162107 73 5.103134 0 0
4 0.0001482271 50 5.098400 0 0
5 0.0004188049 136 4.908158 0 0
6 0.0002992001 95 4.799027 0 0
mean_tss_dist observed_gene_hits gene_set_size fold_enrichment_hyper
1 73676 5 5 1.307880
2 73676 5 5 1.307880
3 73676 5 5 1.307880
4 38163 4 5 1.046304
5 79190 6 6 1.307880
6 35132 9 9 1.307880
p_value_hyper p_adjust_hyper
1 0.26128140 0.5159287
2 0.26128140 0.5159287
3 0.26128140 0.5159287
4 0.66357982 0.8754404
5 0.19976279 0.4453912
6 0.08926921 0.2738572
- 手动提供基因集或注释数据
如果想提供自定的功能基因集或者在R里面没有现成的注释数据包可用时，可使用此方法：
gs <- read_gmt(url("http://dsigdb.tanlab.org/Downloads/D2_LINCS.gmt"))
df <- read.table(url("https://jokergoo.github.io/rGREAT_suppl/data/GREATv4.genes.hg19.tsv"))
tss <- GRanges(seqnames = df[, 2], ranges = IRanges(df[, 3], df[, 3]), strand = df[, 4], gene_id = df[, 5])
head(tss)
GRanges object with 6 ranges and 1 metadata column:
seqnames ranges strand | gene_id
<Rle> <IRanges> <Rle> | <character>
[1] chr1 69090 + | OR4F5
[2] chr1 367639 + | OR4F29
[3] chr1 622053 - | OR4F16
[4] chr1 861117 + | SAMD11
[5] chr1 894670 - | NOC2L
[6] chr1 895966 + | KLHL17
tss_ext <- extendTSS(tss, gene_id_type = "SYMBOL")
res <- great(gr, gs, tss)
tab <- getEnrichmentTable(res)
- 设定背景集
对于背景的设定，rGREAt提供了两种模式，分别由background和exclude两个参数来控制设定最终的背景集，前者可以直接设定需要的背景集，而后者是从当前使用的背景集中排除一些不想考虑的区间，即exclude参数需要配合tss_source参数来使用。
当然，无论使用哪一种方法都可以重新定义背景集，这从上面的原理图也可以直观地看到，这意味着富集的结果将会受到影响。所以，当不知道如何设定合适的背景集时，默认的参数就是最好的选择。
gap <- getGapFromUCSC("hg38", paste0("chr", c(1:22, "X", "Y")))
# 直接设定背景集
res1 <- great(peak, "GO:BP", background = paste0("chr", 1:22))
# 去除背景集里面的gap区域
res2 <- great(peak, "MSigDB:H", "hg38", exclude = gap)
结束语
目前对于基因组区间做富集的软件，大多都是先基于线性距离将区间注释到基因，然后利用超几何分布来做富集检验，这个原理的前提假设是每个基因独立且被选取到的概率相同。而GREAT则采用不同的方式，直接从基因组区间层面来考虑，则前提假设就是基因组区间均匀分布在基因组上。然而，由于长度的不同，基因在基因组上并不符合均匀分布。所以，从基因组区间到基因的转换会导致基因不会以相等的概率被挑选出来。例如，当所有区间都远离一个基因时，该基因被挑选的可能性就会很低；然而，当一个基因附近有一组区间时，它就更有可能被挑选出来；当基因长度较长时，也更有可能被挑选出来。因而，GREAT采用先将基因转换为特定的区间，然后使用二项分布来做富集检验。
友情提示：version1.6.0及之前本版是基于在线网站做的分析，不仅基因组本版受限，也不支持自定义背景集。如果想使用rGREAT，可以安装最新版本来使用。从上面的示例可以看出，用rGREAT来做富集还是挺简单的，且不仅有二项分布的检验结果，也保留了超几何分布的检验结果。虽然GREAT的最初目标是将生物学功能与顺式调控元件相关联，例如转录因子结合位点(TFBS)，但是它的算法允许其扩展到任何类型的基因组区间。一句话形容，great！
往期回顾
利用UCSC预测启动子序列可能结合的转录因子
Vision | scRNA细胞相似性 + 差异signature
HiC | contacts vs distance
hdWGCNA | 单细胞数据共表达网络分析
bed基因注释
