參考學(xué)習(xí)資料：https://www.bioconductor.org/packages/release/bioc/vignettes/GeneOverlap/inst/doc/GeneOverlap.pdf

1.數(shù)據(jù)準(zhǔn)備

Overlapping的基因列表可以揭示生物學(xué)意義，并可能形成新的假說耳胎。例如，組蛋白修飾是一個(gè)重要的細(xì)胞機(jī)制胰默，可以包裝和重新包裝染色質(zhì)憎瘸。通過使染色質(zhì)結(jié)構(gòu)更密集或更疏松荣暮，起到開啟或關(guān)閉基因表達(dá)的作用硬纤。組蛋白H3(H3K4me3)第4位賴氨酸(H3K4me3)上的三甲基化與基因激活有關(guān)敞峭，其全基因組富集可以通過ChIP-seq實(shí)驗(yàn)進(jìn)行定位。由于H3K4me3的激活作用球订，如果我們將H3K4me3結(jié)合的基因與高表達(dá)的基因重疊后裸，我們應(yīng)該預(yù)計(jì)會有積極的關(guān)聯(lián)。類似地冒滩，我們可以在不同組蛋白修飾的基因列表與不同表達(dá)組的基因列表之間進(jìn)行這種Overlap微驶，從而確定每個(gè)組蛋白修飾在基因調(diào)控中的作用。
這個(gè)問題可以表示為超幾何分布或結(jié)合表(這可以通過費(fèi)舍爾精確檢驗(yàn)來解決开睡；請參閱GeneOverlap文檔)因苹。這項(xiàng)任務(wù)在生物學(xué)研究中非常常見较店，GeneOverlap是為了操作數(shù)據(jù)、執(zhí)行測試并可視化Overlapping的結(jié)果容燕。
GeneOverlap包由兩個(gè)類組成：GeneOverlap和GeneOverlipMatrix。GeneOverlap類執(zhí)行測試婚度，并充當(dāng)GeneOverlipMatrix類的構(gòu)建塊蘸秘。
首先是安裝和加載包，因?yàn)槭珍浽赽ioconductor蝗茁，需要用到相應(yīng)的安裝函數(shù)：

rm(list = ls())
#BiocManager::install("GeneOverlap")
library(GeneOverlap)
data(GeneOverlap)
sapply(hESC.ChIPSeq.list, length)
sapply(hESC.RNASeq.list, length)
gs.RNASeq

若要使用GeneOverlay類醋虏，請創(chuàng)建兩個(gè)表示基因名稱的字符向量。一種簡單的方法是使用read.table("filename.txt")從文本文件中讀取它們哮翘。
該包有內(nèi)置數(shù)據(jù)集包括幾個(gè)基因列表數(shù)據(jù),可以用data()加載【苯溃現(xiàn)在，讓我們看看數(shù)據(jù)中包含的內(nèi)容：有三個(gè)對象饭寺。hESC.ChIPSeq.list和hESC.RNASeq.list對象包含來自Chip-seq和RNA-seq實(shí)驗(yàn)的基因列表阻课。gs.RNASeq變量包含基因組背景中的基因數(shù)量。后面有介紹它們是如何創(chuàng)建的艰匙。先來看看基因列表中有多少個(gè)基因限煞。

> 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基因集包含四個(gè)不同組蛋白標(biāo)記的基因列表：H3K4me3、H3K9me3健霹、H3K27me3和H3K36me3旺上；RNA-seq基因集包含三個(gè)不同表達(dá)水平的基因列表：高、中糖埋、低表達(dá)水平宣吱。有兩個(gè)組蛋白標(biāo)記與基因激活相關(guān)：H3K4me3和H3K36me3，另外兩個(gè)與基因抑制相關(guān)：H3K9me3和H3K27me3阶捆。

2.檢驗(yàn)兩個(gè)基因列表之間的GeneOverlap

我們想要探索激活標(biāo)記-H3K4me3是否與高表達(dá)的基因相關(guān)凌节。首先，讓我們構(gòu)造一個(gè)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構(gòu)造函數(shù)已經(jīng)為我們做了一些基本的數(shù)據(jù)統(tǒng)計(jì)倍奢，比如交叉點(diǎn)的數(shù)量。為了檢驗(yàn)關(guān)聯(lián)性的統(tǒng)計(jì)顯著性垒棋，我們做了GeneOverlap類卒煞，將問題描述為檢驗(yàn)兩個(gè)變量是否獨(dú)立，這可以表示為一個(gè)列聯(lián)表(contingency table)叼架，然后使用Fisher的精確檢驗(yàn)來尋找統(tǒng)計(jì)顯著性畔裕。這里的P值為零衣撬，這意味著overlap非常重要。費(fèi)舍爾精確檢驗(yàn)還給出了代表關(guān)聯(lián)強(qiáng)度的odds ratio扮饶。如果odds ratio等于或小于1具练，則兩個(gè)列表之間沒有關(guān)聯(lián)。如果odds ratio遠(yuǎn)大于1甜无，那么這種關(guān)聯(lián)就很強(qiáng)扛点。該類還計(jì)算Jaccard指數(shù)，該指數(shù)衡量兩個(gè)列表之間的相似性岂丘。Jaccard指數(shù)在0和1之間變化陵究，0表示兩者之間沒有相似性，1表示兩者相同奥帘。要顯示更多詳細(xì)信息铜邮，請使用print函數(shù)。

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指數(shù)為0.4松蒜。print函數(shù)還顯示列聯(lián)表（contingency table），該表說明了問題是如何形成的已旧。

此外牍鞠，我們還想知道H3K4me3是否與低表達(dá)的基因有關(guān)。修改參數(shù)：

go.obj <- newGeneOverlap(hESC.ChIPSeq.list$H3K4me3,
                         hESC.RNASeq.list$"Exp Low",
                         genome.size=gs.RNASeq)
go.obj <- testGeneOverlap(go.obj)
print(go.obj)

結(jié)果為

> 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

與高表達(dá)的基因相比评姨，P值現(xiàn)在是1难述，odds ratio為0.1。因此吐句，H3K4me3與基因抑制無關(guān)胁后。
一旦創(chuàng)建了一個(gè)GeneOverlap對象，就可以使用幾個(gè)訪問器來提取它的節(jié)點(diǎn)值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

也可以在創(chuàng)建對象后更改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后，對象進(jìn)入未檢驗(yàn)狀態(tài)文虏。所以我們需要重新檢驗(yàn)一下：

> 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)

結(jié)果如下：

> 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 

正如預(yù)期的那樣氧秘，基因組大小越大年鸳，P值越顯著。這是因?yàn)橥柘啵S著 "universe"的增長搔确，兩個(gè)基因列表偶然重疊的可能性變得越來越小。

3.Visualizing all pairwise overlaps

當(dāng)需要比較兩個(gè)基因集合，每個(gè)都有一個(gè)或多個(gè)列表時(shí)膳算，手動比較它們的效率會相當(dāng)?shù)妥丁？梢詷?gòu)造矩陣涕蜂，其中行表示來自一個(gè)集合的列表华匾，而列表示來自另一個(gè)集合的列表。矩陣的每個(gè)元素都是一個(gè)GeneOverlay對象机隙。要完全可視化這些重疊信息瘦真，可以使用熱圖。為了說明這一點(diǎn)黍瞧，我們希望將Chip-seq的所有基因列表與rna-seq的所有基因列表進(jìn)行比較。

gom.obj <- newGOM(hESC.ChIPSeq.list, hESC.RNASeq.list,
                  gs.RNASeq)
drawHeatmap(gom.obj)

圖1

newGOM構(gòu)造函數(shù)使用兩個(gè)命名列表創(chuàng)建一個(gè)新的GeneOverlipMatrix對象原杂。然后印颤，我們需要做的就是使用draHeatmap函數(shù)來顯示它。在熱圖中穿肄，顏色鍵表示odds ratios年局，重要的p值疊加在網(wǎng)格上。
要從GeneOverlipMatrix對象檢索信息咸产，可以使用兩個(gè)重要的訪問器矢否，分別稱為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可以以嵌套列表的形式獲取每個(gè)比較的基因列表：外部列表表示列僵朗，內(nèi)部列表表示行:

inter.nl <- getNestedList(gom.obj, name="intersection")
str(inter.nl)

另一個(gè)重要的訪問器是方法"["，它允許以類似于矩陣的方式檢索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和高表達(dá)基因之間的比較的GeneOverlap對象验庙。還可以使用標(biāo)簽獲取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 還可以對一組基因進(jìn)行自我比較社牲。例如粪薛，如果我們想知道Chip-seq基因列表之間是如何關(guān)聯(lián)的，我們可以這樣做

gom.self <- newGOM(hESC.ChIPSeq.list,
                   genome.size=gs.RNASeq)
drawHeatmap(gom.self) 

圖2

僅使用上三角矩陣搏恤。
還可以更改顏色數(shù)量和用于熱圖的顏色违寿。例如：

drawHeatmap(gom.self, ncolused=5, grid.col="Blues", note.col="black")

圖3

或者改為顯示Jaccard指數(shù)值

drawHeatmap(gom.self, what="Jaccard", ncolused=3, grid.col="Oranges",
            note.col="black")

圖4

4.數(shù)據(jù)源及其處理

這里使用的實(shí)驗(yàn)數(shù)據(jù)是從ENCODE [ENCODE Consortium, ]項(xiàng)目的網(wǎng)站下載的。ChIP-seq和RNA-seq樣品均來源于人胚胎干細(xì)胞(hESC)熟空。使用Bowtie[Langmead et al., 2009]和Tophat [Trapnell et al., 2009]將原始讀取文件與參考基因組進(jìn)行比對藤巢。
Cufflinks [Trapnell et al., 2010]被用來從RNA-seq數(shù)據(jù)中估計(jì)基因表達(dá)水平。只保留了蛋白質(zhì)編碼基因以供進(jìn)一步分析息罗。用FPKM狀態(tài)對基因進(jìn)行篩選菌瘪，只保留“OK”狀態(tài)的基因。這給我們留下了21,196個(gè)編碼基因，它們的FPKM值得到了可靠的估計(jì)俏扩。然后將這些基因分為高(FPKM>10)糜工、中(FPKM>1和<=10)和低(FPKM<=1)三組。每組中重復(fù)的基因ID在測試前被刪除录淡。
對于ChIP-seq捌木，使用來自IP處理的一個(gè)復(fù)制和來自輸入控制的一個(gè)復(fù)制。
使用MACS2(v2.0.10)使用默認(rèn)參數(shù)值執(zhí)行Peak calling嫉戚。然后刨裆，Peak列表使用diffReps包進(jìn)行區(qū)域注釋[Shen et al., 2013]。
然后提取在基因區(qū)和啟動子區(qū)域出現(xiàn)Peak的基因彬檀。使用上面獲得的RNA-seq基因列表進(jìn)一步篩選基因帆啃。

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