参考学习资料：https://www.bioconductor.org/packages/release/bioc/vignettes/GeneOverlap/inst/doc/GeneOverlap.pdf
1.数据准备
Overlapping的基因列表可以揭示生物学意义，并可能形成新的假说。例如，组蛋白修饰是一个重要的细胞机制，可以包装和重新包装染色质。通过使染色质结构更密集或更疏松，起到开启或关闭基因表达的作用。组蛋白H3(H3K4me3)第4位赖氨酸(H3K4me3)上的三甲基化与基因激活有关，其全基因组富集可以通过ChIP-seq实验进行定位。由于H3K4me3的激活作用，如果我们将H3K4me3结合的基因与高表达的基因重叠，我们应该预计会有积极的关联。类似地，我们可以在不同组蛋白修饰的基因列表与不同表达组的基因列表之间进行这种Overlap，从而确定每个组蛋白修饰在基因调控中的作用。
这个问题可以表示为超几何分布或结合表(这可以通过费舍尔精确检验来解决；请参阅GeneOverlap文档)。这项任务在生物学研究中非常常见，GeneOverlap是为了操作数据、执行测试并可视化Overlapping的结果。
GeneOverlap包由两个类组成：GeneOverlap和GeneOverlipMatrix。GeneOverlap类执行测试，并充当GeneOverlipMatrix类的构建块。
首先是安装和加载包，因为收录在bioconductor，需要用到相应的安装函数：
rm(list = ls())
#BiocManager::install("GeneOverlap")
library(GeneOverlap)
data(GeneOverlap)
sapply(hESC.ChIPSeq.list, length)
sapply(hESC.RNASeq.list, length)
gs.RNASeq
若要使用GeneOverlay类，请创建两个表示基因名称的字符向量。一种简单的方法是使用read.table("filename.txt")从文本文件中读取它们。
该包有内置数据集包括几个基因列表数据,可以用data()加载。现在，让我们看看数据中包含的内容：有三个对象。hESC.ChIPSeq.list和hESC.RNASeq.list对象包含来自Chip-seq和RNA-seq实验的基因列表。gs.RNASeq变量包含基因组背景中的基因数量。后面有介绍它们是如何创建的。先来看看基因列表中有多少个基因。
> sapply(hESC.ChIPSeq.list, length)
H3K4me3 H3K9me3 H3K27me3 H3K36me3
13448 297 3853 4533
> sapply(hESC.RNASeq.list, length)
Exp High Exp Medium Exp Low
6540 6011 8684
> gs.RNASeq
[1] 21196
在GeneOverlap中，我们将基因列表集合称为表示为命名列表的基因集。在这里我们可以看到，ChIP-seq基因集包含四个不同组蛋白标记的基因列表：H3K4me3、H3K9me3、H3K27me3和H3K36me3；RNA-seq基因集包含三个不同表达水平的基因列表：高、中、低表达水平。有两个组蛋白标记与基因激活相关：H3K4me3和H3K36me3，另外两个与基因抑制相关：H3K9me3和H3K27me3。
2.检验两个基因列表之间的GeneOverlap
我们想要探索激活标记-H3K4me3是否与高表达的基因相关。首先，让我们构造一个GeneOverlap对象
go.obj <- newGeneOverlap(hESC.ChIPSeq.list$H3K4me3,
hESC.RNASeq.list$"Exp High",
genome.size=gs.RNASeq)
go.obj
显示如下：
> go.obj
GeneOverlap object:
listA size=13448
listB size=6540
Intersection size=5916
Overlap testing has not been performed yet.
正如我们所看到的，GeneOverlap构造函数已经为我们做了一些基本的数据统计，比如交叉点的数量。为了检验关联性的统计显著性，我们做了GeneOverlap类，将问题描述为检验两个变量是否独立，这可以表示为一个列联表(contingency table)，然后使用Fisher的精确检验来寻找统计显著性。这里的P值为零，这意味着overlap非常重要。费舍尔精确检验还给出了代表关联强度的odds ratio。如果odds ratio等于或小于1，则两个列表之间没有关联。如果odds ratio远大于1，那么这种关联就很强。该类还计算Jaccard指数，该指数衡量两个列表之间的相似性。Jaccard指数在0和1之间变化，0表示两者之间没有相似性，1表示两者相同。要显示更多详细信息，请使用print函数。
go.obj <- testGeneOverlap(go.obj)
go.obj
> go.obj
GeneOverlap object:
listA size=13448
listB size=6540
Intersection size=5916
Overlapping p-value=0e+00
Jaccard Index=0.4
> print(go.obj)
Detailed information about this GeneOverlap object:
listA size=13448, e.g. ENSG00000187634 ENSG00000188976 ENSG00000187961
listB size=6540, e.g. ENSG00000000003 ENSG00000000419 ENSG00000000460
Intersection size=5916, e.g. ENSG00000188976 ENSG00000187961 ENSG00000187608
Union size=14072, e.g. ENSG00000187634 ENSG00000188976 ENSG00000187961
Genome size=21196
# Contingency Table:
notA inA
notB 7124 7532
inB 624 5916
Overlapping p-value=0e+00
Odds ratio=9.0
Overlap tested using Fisher's exact test (alternative=greater)
Jaccard Index=0.4
显示odds ratio为9.0，Jaccard指数为0.4。print函数还显示列联表（contingency table），该表说明了问题是如何形成的。
此外，我们还想知道H3K4me3是否与低表达的基因有关。修改参数：
go.obj <- newGeneOverlap(hESC.ChIPSeq.list$H3K4me3,
hESC.RNASeq.list$"Exp Low",
genome.size=gs.RNASeq)
go.obj <- testGeneOverlap(go.obj)
print(go.obj)
结果为
> print(go.obj)
Detailed information about this GeneOverlap object:
listA size=13448, e.g. ENSG00000187634 ENSG00000188976 ENSG00000187961
listB size=8684, e.g. ENSG00000000005 ENSG00000000938 ENSG00000000971
Intersection size=2583, e.g. ENSG00000187583 ENSG00000187642 ENSG00000215014
Union size=19549, e.g. ENSG00000187634 ENSG00000188976 ENSG00000187961
Genome size=21196
# Contingency Table:
notA inA
notB 1647 10865
inB 6101 2583
Overlapping p-value=1
Odds ratio=0.1
Overlap tested using Fisher's exact test (alternative=greater)
Jaccard Index=0.1
与高表达的基因相比，P值现在是1，odds ratio为0.1。因此，H3K4me3与基因抑制无关。
一旦创建了一个GeneOverlap对象，就可以使用几个访问器来提取它的节点值slots。例如
> head(getIntersection(go.obj))
[1] "ENSG00000187583" "ENSG00000187642" "ENSG00000215014" "ENSG00000142609" "ENSG00000203301"
[6] "ENSG00000215912"
> getOddsRatio(go.obj)
[1] 0.06418943
> getContbl(go.obj)
notA inA
notB 1647 10865
inB 6101 2583
> getContbl(go.obj)
notA inA
notB 1647 10865
inB 6101 2583
> getGenomeSize(go.obj)
[1] 21196
也可以在创建对象后更改slots "listA", "listB" 和 "genome.size"。例如
setListA(go.obj) <- hESC.ChIPSeq.list"Exp Medium"
go.obj
> go.obj
GeneOverlap object:
listA size=3853
listB size=6011
Intersection size=1549
Overlap testing has not been performed yet.
更换上述任一slots后，对象进入未检验状态。所以我们需要重新检验一下：
> go.obj <- testGeneOverlap(go.obj)
> go.obj
GeneOverlap object:
listA size=3853
listB size=6011
Intersection size=1549
Overlapping p-value=2.8e-69
Jaccard Index=0.2
我们还可以改变基因组大小，看看p值是如何随之变化的。
v.gs <- c(12e3, 14e3, 16e3, 18e3, 20e3)
setNames(sapply(v.gs, function(g) {
setGenomeSize(go.obj) <- g
go.obj <- testGeneOverlap(go.obj)
getPval(go.obj)}), v.gs)
结果如下：
> setNames(sapply(v.gs, function(g) {
+ setGenomeSize(go.obj) <- g
+ go.obj <- testGeneOverlap(go.obj)
+ getPval(go.obj)}), v.gs)
12000 14000 16000 18000 20000
1.000000e+00 9.999746e-01 6.039536e-05 9.006016e-24 5.589067e-51
正如预期的那样，基因组大小越大，P值越显著。这是因为，随着 "universe"的增长，两个基因列表偶然重叠的可能性变得越来越小。
3.Visualizing all pairwise overlaps
当需要比较两个基因集合，每个都有一个或多个列表时，手动比较它们的效率会相当低。可以构造矩阵，其中行表示来自一个集合的列表，而列表示来自另一个集合的列表。矩阵的每个元素都是一个GeneOverlay对象。要完全可视化这些重叠信息，可以使用热图。为了说明这一点，我们希望将Chip-seq的所有基因列表与rna-seq的所有基因列表进行比较。
gom.obj <- newGOM(hESC.ChIPSeq.list, hESC.RNASeq.list,
gs.RNASeq)
drawHeatmap(gom.obj)
newGOM构造函数使用两个命名列表创建一个新的GeneOverlipMatrix对象。然后，我们需要做的就是使用draHeatmap函数来显示它。在热图中，颜色键表示odds ratios，重要的p值叠加在网格上。
要从GeneOverlipMatrix对象检索信息，可以使用两个重要的访问器，分别称为getMatrix和getNestedList。例如，getMatrix accessor以矩阵形式获取诸如p值之类的信息
> getMatrix(gom.obj, name="pval")
Exp High Exp Medium Exp Low
H3K4me3 0.0000000 0.000000e+00 1.000000e+00
H3K9me3 0.9996654 7.071538e-13 9.999497e-01
H3K27me3 1.0000000 2.797009e-69 1.061420e-09
H3K36me3 0.0000000 9.999998e-01 1.000000e+00
或者提取odds ratios：
> getMatrix(gom.obj, "odds.ratio")
Exp High Exp Medium Exp Low
H3K4me3 8.9672775 3.8589147 0.06418943
H3K9me3 0.6366248 2.3460391 0.62254173
H3K27me3 0.3338513 1.9408115 1.24113618
H3K36me3 10.9219512 0.8265542 0.02145576
getNestedListaccessor可以以嵌套列表的形式获取每个比较的基因列表：外部列表表示列，内部列表表示行:
inter.nl <- getNestedList(gom.obj, name="intersection")
str(inter.nl)
另一个重要的访问器是方法"["，它允许以类似于矩阵的方式检索GeneOverlap对象。例如,
> go.k4.high <- gom.obj[1, 1]
> go.k4.high
GeneOverlap object:
listA size=13448
listB size=6540
Intersection size=5916
Overlapping p-value=0e+00
Jaccard Index=0.4
获取表示H3K4me3和高表达基因之间的比较的GeneOverlap对象。还可以使用标签获取GeneOverlap对象，例如
> gom.obj["H3K9me3", "Exp Medium"]
GeneOverlap object:
listA size=297
listB size=6011
Intersection size=142
Overlapping p-value=7.1e-13
Jaccard Index=0.0
GeneOverlapMatrix 还可以对一组基因进行自我比较。例如，如果我们想知道Chip-seq基因列表之间是如何关联的，我们可以这样做
gom.self <- newGOM(hESC.ChIPSeq.list,
genome.size=gs.RNASeq)
drawHeatmap(gom.self)
仅使用上三角矩阵。
还可以更改颜色数量和用于热图的颜色。例如：
drawHeatmap(gom.self, ncolused=5, grid.col="Blues", note.col="black")
或者改为显示Jaccard指数值
drawHeatmap(gom.self, what="Jaccard", ncolused=3, grid.col="Oranges",
note.col="black")
4.数据源及其处理
这里使用的实验数据是从ENCODE [ENCODE Consortium, ]项目的网站下载的。ChIP-seq和RNA-seq样品均来源于人胚胎干细胞(hESC)。使用Bowtie[Langmead et al., 2009]和Tophat [Trapnell et al., 2009]将原始读取文件与参考基因组进行比对。
Cufflinks [Trapnell et al., 2010]被用来从RNA-seq数据中估计基因表达水平。只保留了蛋白质编码基因以供进一步分析。用FPKM状态对基因进行筛选，只保留“OK”状态的基因。这给我们留下了21,196个编码基因，它们的FPKM值得到了可靠的估计。然后将这些基因分为高(FPKM>10)、中(FPKM>1和<=10)和低(FPKM<=1)三组。每组中重复的基因ID在测试前被删除。
对于ChIP-seq，使用来自IP处理的一个复制和来自输入控制的一个复制。
使用MACS2(v2.0.10)使用默认参数值执行Peak calling。然后，Peak列表使用diffReps包进行区域注释[Shen et al., 2013]。
然后提取在基因区和启动子区域出现Peak的基因。使用上面获得的RNA-seq基因列表进一步筛选基因。
5 SessionInfo
> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] GeneOverlap_1.22.0 ggthemes_4.2.0 patchwork_1.0.0 forcats_0.4.0 stringr_1.4.0
[6] dplyr_0.8.3 purrr_0.3.3 readr_1.3.1 tidyr_1.0.0 tibble_2.1.3
[11] ggplot2_3.3.0 tidyverse_1.3.0
loaded via a namespace (and not attached):
[1] nlme_3.1-141 bitops_1.0-6 fs_1.3.1 usethis_1.5.1 lubridate_1.7.4
[6] devtools_2.2.2 RColorBrewer_1.1-2 httr_1.4.1 rprojroot_1.3-2 tools_3.6.1
[11] backports_1.1.5 utf8_1.1.4 R6_2.4.1 KernSmooth_2.23-16 DBI_1.0.0
[16] mgcv_1.8-29 colorspace_1.4-1 withr_2.1.2 tidyselect_0.2.5 prettyunits_1.0.2
[21] processx_3.4.1 curl_4.2 compiler_3.6.1 cli_2.0.1 rvest_0.3.5
[26] xml2_1.2.2 desc_1.2.0 isoband_0.2.0 labeling_0.3 caTools_1.17.1.2
[31] scales_1.1.0 callr_3.3.2 digest_0.6.23 pkgconfig_2.0.3 sessioninfo_1.1.1
[36] dbplyr_1.4.2 rlang_0.4.4.9000 readxl_1.3.1 rstudioapi_0.10 farver_2.0.1
[41] generics_0.0.2 jsonlite_1.6 gtools_3.8.1 magrittr_1.5 Matrix_1.2-17
[46] Rcpp_1.0.3 munsell_0.5.0 fansi_0.4.0 lifecycle_0.1.0 stringi_1.4.3
[51] yaml_2.2.0 MASS_7.3-51.4 pkgbuild_1.0.6 gplots_3.0.1.1 grid_3.6.1
[56] gdata_2.18.0 crayon_1.3.4 lattice_0.20-38 haven_2.2.0 splines_3.6.1
[61] hms_0.5.2 zeallot_0.1.0 ps_1.3.0 pillar_1.4.2 pkgload_1.0.2
[66] reprex_0.3.0 glue_1.3.1 remotes_2.1.1 BiocManager_1.30.9 modelr_0.1.5
[71] vctrs_0.2.0 testthat_2.3.0 cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
[76] broom_0.5.4 viridisLite_0.3.0 memoise_1.1.0 ellipsis_0.3.0
