背景
??GREAT(Genomic Regions Enrichment of Annotations Tool)是一種廣泛使用的基因組區(qū)域功能富集工具，該工具于2010年由斯坦福大學(xué)開發(fā)。然而，作為在線工具，其存在注釋數(shù)據(jù)過時、支持的物種和功能基因集數(shù)量少以及用戶不可擴展等局限性。因此，有人就開發(fā)了一個本地實現(xiàn)GREAT算法的R包rGREAT，其默認(rèn)支持600多個物種和大量的功能基因集，同時也支持用戶自備基因集和物種基因組。此外，該包還實現(xiàn)了處理背景區(qū)域的通用方法。
??基因組學(xué)和表觀組學(xué)研究通常會生成許多感興趣的基因組區(qū)域列表，例如來自全基因組或外顯子組測序數(shù)據(jù)的單核苷酸變異(SNV)、ChIP-seq 的特定染色質(zhì)修飾的peak或基因組甲基化測序的差異甲基化區(qū)域 (DMR)。下一步的分析自然是將生物學(xué)功能與這些基因組區(qū)域聯(lián)系起來。一種廣泛使用的方法是首先將基因組區(qū)域注釋到最近的基因，然后針對特定功能的基因集做富集分析。
??GREAT在做基因組區(qū)域富集時，考慮了基因在基因組上的位置分布與長度，采用了不同的策略，如上圖所示。對于給定功能基因集中的基因，首先，將基因的TSS上下游分別延伸5kb和1kb來建立一個基礎(chǔ)結(jié)構(gòu)域；然后，將基礎(chǔ)結(jié)構(gòu)域的上下游再繼續(xù)延申至最大1mb，或者到達(dá)它鄰近基因的基礎(chǔ)結(jié)構(gòu)域，如此每個基因相當(dāng)于被轉(zhuǎn)化成了一個區(qū)間；最后，將這些轉(zhuǎn)化的區(qū)間進行合并形成一個沒有overlap的區(qū)間集，相當(dāng)于將特定生物學(xué)功能相關(guān)的基因集轉(zhuǎn)化為了“功能區(qū)間集”，然后使用二項分布來計算輸入?yún)^(qū)域集是否在功能區(qū)間集中富集。
分析
- 使用R包里基因集和注釋數(shù)據(jù)
??一般常見物種如人、小鼠等，在R里面都可以找到相應(yīng)基因注釋數(shù)據(jù)庫，如TxDb.Hsapiens.UCSC.hg38.knownGene。功能基因集也可直接使用R數(shù)據(jù)包里面的集合，如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
- 手動提供基因集或注釋數(shù)據(jù)
??如果想提供自定的功能基因集或者在R里面沒有現(xiàn)成的注釋數(shù)據(jù)包可用時，可使用此方法：
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)
- 設(shè)定背景集
??對于背景的設(shè)定，rGREAt提供了兩種模式，分別由background和exclude兩個參數(shù)來控制設(shè)定最終的背景集，前者可以直接設(shè)定需要的背景集，而后者是從當(dāng)前使用的背景集中排除一些不想考慮的區(qū)間，即exclude參數(shù)需要配合tss_source參數(shù)來使用。
??當(dāng)然，無論使用哪一種方法都可以重新定義背景集，這從上面的原理圖也可以直觀地看到，這意味著富集的結(jié)果將會受到影響。所以，當(dāng)不知道如何設(shè)定合適的背景集時，默認(rèn)的參數(shù)就是最好的選擇。
gap <- getGapFromUCSC("hg38", paste0("chr", c(1:22, "X", "Y")))
# 直接設(shè)定背景集
res1 <- great(peak, "GO:BP", background = paste0("chr", 1:22))
# 去除背景集里面的gap區(qū)域
res2 <- great(peak, "MSigDB:H", "hg38", exclude = gap)
結(jié)束語
??目前對于基因組區(qū)間做富集的軟件，大多都是先基于線性距離將區(qū)間注釋到基因，然后利用超幾何分布來做富集檢驗，這個原理的前提假設(shè)是每個基因獨立且被選取到的概率相同。而GREAT則采用不同的方式，直接從基因組區(qū)間層面來考慮，則前提假設(shè)就是基因組區(qū)間均勻分布在基因組上。然而，由于長度的不同，基因在基因組上并不符合均勻分布。所以，從基因組區(qū)間到基因的轉(zhuǎn)換會導(dǎo)致基因不會以相等的概率被挑選出來。例如，當(dāng)所有區(qū)間都遠(yuǎn)離一個基因時，該基因被挑選的可能性就會很低；然而，當(dāng)一個基因附近有一組區(qū)間時，它就更有可能被挑選出來；當(dāng)基因長度較長時，也更有可能被挑選出來。因而，GREAT采用先將基因轉(zhuǎn)換為特定的區(qū)間，然后使用二項分布來做富集檢驗。
??友情提示：version1.6.0及之前本版是基于在線網(wǎng)站做的分析，不僅基因組本版受限，也不支持自定義背景集。如果想使用rGREAT，可以安裝最新版本來使用。從上面的示例可以看出，用rGREAT來做富集還是挺簡單的，且不僅有二項分布的檢驗結(jié)果，也保留了超幾何分布的檢驗結(jié)果。雖然GREAT的最初目標(biāo)是將生物學(xué)功能與順式調(diào)控元件相關(guān)聯(lián)，例如轉(zhuǎn)錄因子結(jié)合位點(TFBS)，但是它的算法允許其擴展到任何類型的基因組區(qū)間。一句話形容，great！
往期回顧
利用UCSC預(yù)測啟動子序列可能結(jié)合的轉(zhuǎn)錄因子
Vision | scRNA細(xì)胞相似性 + 差異signature
HiC | contacts vs distance
hdWGCNA | 單細(xì)胞數(shù)據(jù)共表達(dá)網(wǎng)絡(luò)分析
bed基因注釋
