ChIPseq_In_Bioconductor3.knit

class: center, middle, inverse, title-slide

# ChIPseq In Bioconductor (part3)
<html>
<div style="float:left">

</div>
<hr color='#EB811B' size=1px width=796px>
</html>
### Rockefeller University, Bioinformatics Resource Center
### <a href="http://rockefelleruniversity.github.io/RU_ChIPseq/" class="uri">http://rockefelleruniversity.github.io/RU_ChIPseq/</a>

---

## Data

In todays session we will continue to review the Myc ChIPseq we were working on in our last sessions.

This include Myc ChIPseq for MEL and Ch12 celllines.

Information and files for the [Myc ChIPseq in MEL cell line can be found here](https://www.encodeproject.org/experiments/ENCSR000EUA/)

Information and files for the [Myc ChIPseq in Ch12 cell line can be found here](https://www.encodeproject.org/experiments/ENCSR000ERN/)

---

## Data

In the data directory we have provided peak calls from MACS2 following the processing steps outlined in our last session.

Peak calls for Myc in MEL and Ch12 cellines can be found in

**data/peaks/**

* **data/peaks/Mel_1_peaks.xls**
* **data/peaks/Mel_2_peaks.xls**
* **data/peaks/Ch12_1_peaks.xls**
* **data/peaks/Ch12_1_peaks.xls**

---

## ChIP peaks in R

In our last session we reviewed how we can [identify putative transciption factor binding sites using peak calling programs such as MACS2.](https://rockefelleruniversity.github.io/RU_ChIPseq/presentations/slides/ChIPseq_In_Bioconductor2.html#38)

```r
library(GenomicRanges)
macsPeaks <- "data/peaks/Mel_1_peaks.xls"
macsPeaks_DF <- read.delim(macsPeaks,comment.char="#")
macsPeaks_GR <- GRanges(seqnames=macsPeaks_DF[,"chr"],
                        IRanges(macsPeaks_DF[,"start"],macsPeaks_DF[,"end"]))
mcols(macsPeaks_GR) <- macsPeaks_DF[,c("abs_summit", "fold_enrichment")]
macsPeaks_GR[1:5,]
```

```
## GRanges object with 5 ranges and 2 metadata columns:
##       seqnames          ranges strand | abs_summit fold_enrichment
##          <Rle>       <IRanges>  <Rle> |  <integer>       <numeric>
##   [1]     chr1 4785370-4785641      * |    4785565         5.23296
##   [2]     chr1 5082993-5083216      * |    5083122         4.26417
##   [3]     chr1 7397557-7398093      * |    7397835         9.84580
##   [4]     chr1 9545120-9545629      * |    9545400        11.85079
##   [5]     chr1 9700895-9701037      * |    9700978         5.14291
##   -------
##   seqinfo: 21 sequences from an unspecified genome; no seqlengths
```

---
## Annotation of peaks to genes

Since transcription factors, as suggested in name, may regulate the transcription of their target genes, we used the **ChIPseeker package** to associate our peaks, representing potential transcription factor binding events, to their overlapping or closest mm10 genes.

```r
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(ChIPseeker)
peakAnno <- annotatePeak(macsPeaks_GR, tssRegion=c(-1000, 1000), 
                         TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene, 
                         annoDb="org.Mm.eg.db")
```

```
## >> preparing features information...		 2021-08-02 03:22:33 PM 
## >> identifying nearest features...		 2021-08-02 03:22:33 PM 
## >> calculating distance from peak to TSS...	 2021-08-02 03:22:34 PM 
## >> assigning genomic annotation...		 2021-08-02 03:22:34 PM 
## >> adding gene annotation...			 2021-08-02 03:22:37 PM 
## >> assigning chromosome lengths			 2021-08-02 03:22:37 PM 
## >> done...					 2021-08-02 03:22:37 PM
```

---

## Annotation of peaks to genes

This allowed us to produce a GRanges or data.frame of our peaks and their predicted target genes.

```r
annotatedPeaksGR <- as.GRanges(peakAnno)
annotatedPeaksDF <- as.data.frame(peakAnno)
annotatedPeaksDF[1:2, ]
```

```
##   seqnames   start     end width strand abs_summit fold_enrichment annotation
## 1     chr1 4785370 4785641   272      *    4785565         5.23296   Promoter
## 2     chr1 5082993 5083216   224      *    5083122         4.26417   Promoter
##   geneChr geneStart geneEnd geneLength geneStrand geneId          transcriptId
## 1       1   4783572 4785692       2121          2  27395  ENSMUST00000132625.1
## 2       1   5083080 5162529      79450          1 108664 ENSMUST00000044369.12
##   distanceToTSS            ENSEMBL  SYMBOL
## 1            51 ENSMUSG00000033845  Mrpl15
## 2             0 ENSMUSG00000033793 Atp6v1h
##                                          GENENAME
## 1             mitochondrial ribosomal protein L15
## 2 ATPase, H+ transporting, lysosomal V1 subunit H
```

---
class: inverse, center, middle

# Gene Set Enrichment

---

## Gene Set testing

Transcription factors or epigenetic marks may act on specific sets of genes grouped by a common biological feature (shared Biological function, common regulation in RNAseq experiment etc).

A frequent step in ChIPseq analysis is to test whether common gene sets are enriched for transcription factor binding or epigenetic marks.

Sources of well curated gene sets include [GO consortium](http://geneontology.org/) (gene's function, biological process and cellular localisation), [REACTOME](http://www.reactome.org/) (Biological Pathways) and [MsigDB](http://software.broadinstitute.org/gsea/msigdb/) (Computationally and Experimentally derived).

---
## Gene Set testing for ChIPseq

Gene set enrichment testing may be performed on the sets of genes with peaks associated to them. In this example we will consider genes with peaks within 1000bp of a gene's TSS.

We will not access these database libraries directly in testing but will use other R/Bioconductor libraries which make extensive use of them.

---
## Gene ontology and gene set testing

To perform gene set testing here, we will use the **clusterProfiler package**.

clusterProfiler provides multiple enrichment functions that allow for comparison of your gene list to known (e.g. GO, KEGG) or custom gene sets. Detailed information about all of the functionality within this package is available [here](http://yulab-smu.top/clusterProfiler-book/).

In this example we use all TSS sites we found to be overlapped by Myc peaks.

The peaks landing in TSS regions will be marked as "Promoter" in the **annotation** column of our annotated GRanges object.

```r
annotatedPeaksGR[1, ]
```

```
## GRanges object with 1 range and 14 metadata columns:
##       seqnames          ranges strand | abs_summit fold_enrichment  annotation
##          <Rle>       <IRanges>  <Rle> |  <integer>       <numeric> <character>
##   [1]     chr1 4785370-4785641      * |    4785565         5.23296    Promoter
##         geneChr geneStart   geneEnd geneLength geneStrand      geneId
##       <integer> <integer> <integer>  <integer>  <integer> <character>
##   [1]         1   4783572   4785692       2121          2       27395
##               transcriptId distanceToTSS            ENSEMBL      SYMBOL
##                <character>     <numeric>        <character> <character>
##   [1] ENSMUST00000132625.1            51 ENSMUSG00000033845      Mrpl15
##                     GENENAME
##                  <character>
##   [1] mitochondrial riboso..
##   -------
##   seqinfo: 21 sequences from mm10 genome
```

---
## Gene ontology and gene set testing

We can extract the unique names of genes with peaks in their TSS by subsetting the annotated GRanges and retrieving gene names from the **geneId** column.

```r
annotatedPeaksGR_TSS <- annotatedPeaksGR[annotatedPeaksGR$annotation == "Promoter",
    ]
genesWithPeakInTSS <- unique(annotatedPeaksGR_TSS$geneId)
genesWithPeakInTSS[1:2]
```

```
## [1] "27395"  "108664"
```

---
## Gene ontology and functional testing

Next we can extract all genes which are included in the TxDb object to use as our universe of genes for pathway enrichment.

```r
allGeneGR <- genes(TxDb.Mmusculus.UCSC.mm10.knownGene)
allGeneGR[1:2, ]
```

```
## GRanges object with 2 ranges and 1 metadata column:
##             seqnames            ranges strand |     gene_id
##                <Rle>         <IRanges>  <Rle> | <character>
##   100009600     chr9 21062393-21073096      - |   100009600
##   100009609     chr7 84935565-84964115      - |   100009609
##   -------
##   seqinfo: 66 sequences (1 circular) from mm10 genome
```

```r
allGeneIDs <- allGeneGR$gene_id
```

---
## Gene ontology and functional testing

Once we have our gene list and the universe of genes in the same format, we can use them in the **enrichGO** function to perform gene ontology analysis

For the **ont** argument, we can choose between the "BP", "MF", and "CC" subontologies, or "ALL" for all three.

```r
library(clusterProfiler)
library(org.Mm.eg.db)
GO_result <- enrichGO(gene = genesWithPeakInTSS, universe = allGeneIDs, OrgDb = org.Mm.eg.db,
    ont = "BP")
```

---
## Gene ontology and functional testing

We now have an **enrichResult** instance. From this object, we can extract a data frame of the most highly enriched gene ontology categories.

```r
GO_result_df <- data.frame(GO_result)
GO_result_df[1:5, ]
```

```
##                    ID                          Description GeneRatio   BgRatio
## GO:0022613 GO:0022613 ribonucleoprotein complex biogenesis  288/4898 405/20741
## GO:0034660 GO:0034660              ncRNA metabolic process  293/4898 444/20741
## GO:0034470 GO:0034470                     ncRNA processing  248/4898 354/20741
## GO:0042254 GO:0042254                  ribosome biogenesis  212/4898 291/20741
## GO:0006364 GO:0006364                      rRNA processing  153/4898 206/20741
##                  pvalue     p.adjust       qvalue
## GO:0022613 1.050936e-92 6.443286e-89 5.183878e-89
## GO:0034660 7.292474e-82 2.235508e-78 1.798554e-78
## GO:0034470 1.300402e-77 2.657589e-74 2.138135e-74
## GO:0042254 2.150571e-71 3.296288e-68 2.651994e-68
## GO:0006364 1.328419e-53 1.628908e-50 1.310520e-50
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             geneID
## GO:0022613                                      59014/19989/27993/73674/55989/57750/22596/17975/67949/66942/20643/13211/15278/67459/67223/215193/353258/67239/56200/67332/108707/110816/18114/216156/20054/20222/103136/70683/69878/52705/432508/56351/18813/64934/67465/52513/57316/103573/18148/52530/216766/216767/19069/12261/18432/67493/104662/192159/276919/216987/217011/78394/56321/217038/110809/12812/67684/107732/19942/68979/67671/72344/217431/20115/71952/15519/217869/69237/217995/30946/208366/70930/28126/94066/105372/108143/66254/72198/218629/20088/52535/27374/75751/27275/72662/14000/223499/78581/239528/67959/109075/12181/55944/27367/20826/74778/72193/75705/21357/66409/106298/224092/57808/68193/67014/106338/207932/224481/67045/213773/66978/67097/20630/19896/72462/54127/27756/53817/15516/224823/66229/53414/18519/72515/12607/67242/20641/108013/69241/225348/67936/67222/20044/68731/20227/67674/107071/28000/59028/19211/54364/107094/66583/18150/18572/13669/227522/20935/27176/227656/118451/64930/227674/227715/66489/75422/13046/75786/98956/66181/271842/320632/20638/67134/66580/67204/16418/68295/67996/228889/66481/68879/69126/69639/229317/56758/97112/74200/56417/70767/56195/68147/229937/70285/100019/230082/20104/12578/72181/20116/69072/67205/75062/230737/214150/14105/67025/69902/50912/74035/66448/57741/20817/433864/75416/17089/65961/231452/94061/100503670/100608/117109/11837/71990/19988/208144/19384/22384/100609/74097/66596/192196/28035/68011/57905/56505/72171/213895/110109/60321/20901/68988/20103/68077/107686/13871/69731/20085/27998/14113/20055/22121/67973/20068/67148/101867/67920/67276/234344/234374/50783/19656/78651/24015/234549/214627/214572/66164/67619/72544/235036/215194/71989/382137/77591/11920/102462/67941/225215/245000/22258/72831/16785/67115/66446/13205/72584/68995/72554/19655/110954/245474/26905/76130/237107/69499/22184/237221
## GO:0034660 59014/17864/19989/27993/73674/212679/57750/23874/69163/17975/67949/66942/226414/15278/67459/381314/107508/67223/77065/215193/73681/67239/56200/108707/110816/18114/216156/20054/28088/70683/52705/17299/216443/18813/64934/52513/57316/74467/71701/18148/104458/66628/52530/194908/68626/54351/67493/104662/216987/78394/217038/70422/20848/69684/19942/52686/68979/76265/72344/68730/217431/20115/328092/66132/211064/328162/74691/69237/217995/30946/68916/208366/69955/70930/105148/28114/100715/105372/66254/72198/105083/67053/68045/20088/67011/105638/56335/18130/72662/110960/14000/67724/223499/78581/68260/239528/109075/12181/20826/74778/72026/21357/66409/74108/70120/57808/106248/52575/106338/207932/224481/67045/213773/19896/57773/72462/54127/22321/68915/224823/66229/53414/18519/224907/72515/72167/240174/76781/58523/225348/67936/15115/107045/68857/67222/20044/70223/68731/20227/19697/67674/107071/109077/70044/59028/54364/107094/66583/18150/18572/72133/227522/13434/74455/27176/227656/68205/227715/66489/20823/75422/98956/66181/271842/69185/66580/66044/16418/228889/66481/67005/69639/229487/56417/229543/71807/70560/20226/229780/28036/68147/108943/229937/70285/100019/109093/230082/69934/74753/230233/20104/12578/100087/72181/68276/20116/69072/66966/67205/214150/107271/67025/69902/50912/74035/71957/78697/75416/78890/17089/211006/65961/231413/94061/100503670/56361/117109/71990/70650/208144/100929/100609/83701/74097/231803/56463/53379/353172/69786/72171/70047/14911/381802/213895/110109/60321/66078/68077/13871/20085/27998/14113/20055/233189/67973/20068/244141/101861/101867/233802/71941/67920/67276/234344/234374/70359/78651/212528/66590/68544/214627/66369/67619/234734/72544/30947/85305/66965/212728/215194/71989/67049/382137/77591/102462/69882/69606/68291/75686/97541/67789/72341/67115/66446/102436/54632/13205/69942/72554/245474/76130/69499/15108
## GO:0034470                                                                                                                                                                                                                                                                                                    59014/19989/27993/73674/57750/69163/17975/67949/66942/15278/67459/67223/77065/215193/73681/67239/56200/108707/110816/18114/216156/20054/28088/70683/52705/17299/18813/64934/52513/57316/74467/18148/66628/52530/68626/54351/67493/104662/216987/78394/217038/70422/20848/19942/52686/68979/76265/72344/68730/217431/20115/66132/211064/328162/69237/217995/30946/68916/208366/69955/70930/28114/105372/66254/72198/67053/68045/20088/67011/105638/56335/18130/72662/14000/67724/223499/78581/68260/239528/109075/12181/20826/74778/72026/21357/66409/74108/57808/106248/52575/106338/207932/224481/67045/213773/19896/57773/72462/54127/224823/66229/53414/18519/224907/72515/72167/240174/58523/225348/67936/68857/67222/20044/68731/20227/67674/107071/109077/70044/59028/54364/107094/66583/18150/18572/72133/227522/13434/74455/27176/227656/68205/227715/66489/20823/75422/98956/66181/271842/69185/66580/16418/228889/66481/67005/69639/56417/229543/229780/28036/68147/108943/229937/70285/100019/230082/69934/74753/230233/20104/12578/100087/72181/68276/20116/69072/66966/67205/214150/67025/69902/50912/74035/71957/78697/75416/78890/17089/65961/231413/94061/100503670/56361/117109/71990/70650/208144/100929/100609/83701/74097/53379/69786/72171/70047/14911/381802/213895/110109/60321/66078/68077/13871/20085/27998/14113/20055/233189/67973/20068/101861/101867/233802/67920/67276/234344/234374/70359/78651/212528/214627/66369/67619/234734/72544/30947/66965/212728/215194/71989/67049/382137/77591/102462/69882/69606/68291/72341/67115/66446/54632/13205/69942/72554/245474/76130/69499/15108
## GO:0042254                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              59014/19989/27993/73674/55989/57750/22596/17975/67949/66942/15278/67459/67223/215193/353258/67239/56200/108707/110816/18114/216156/20054/103136/70683/52705/18813/64934/52513/57316/103573/18148/52530/216767/19069/12261/18432/67493/104662/216987/217011/78394/56321/217038/107732/19942/68979/67671/72344/217431/20115/71952/69237/217995/30946/208366/70930/28126/94066/105372/66254/72198/218629/20088/52535/75751/72662/14000/223499/78581/109075/12181/27367/20826/74778/66409/106298/224092/57808/68193/67014/106338/207932/224481/67045/213773/67097/19896/72462/54127/224823/66229/53414/18519/72515/12607/225348/67936/67222/20044/68731/20227/67674/107071/59028/19211/54364/107094/66583/18150/18572/227522/20935/27176/227656/118451/64930/227674/227715/66489/75422/98956/66181/271842/67134/66580/16418/228889/66481/69126/69639/229317/97112/68147/229937/70285/100019/230082/20104/12578/72181/20116/69072/67205/230737/67025/69902/50912/74035/66448/57741/433864/75416/17089/65961/231452/94061/100503670/100608/117109/11837/71990/19988/208144/19384/100609/74097/66596/72171/213895/110109/60321/20103/68077/13871/20085/27998/14113/20055/67973/20068/101867/67920/67276/234344/234374/78651/24015/234549/214627/66164/67619/72544/235036/215194/71989/382137/77591/102462/67941/225215/72831/16785/67115/66446/13205/72584/72554/110954/245474/76130/237107/69499
## GO:0006364                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     59014/19989/27993/73674/57750/17975/67949/66942/15278/67459/67223/215193/67239/56200/108707/110816/18114/216156/20054/70683/52705/18813/64934/52513/57316/18148/52530/67493/104662/216987/78394/217038/19942/68979/72344/217431/20115/69237/217995/30946/208366/70930/105372/66254/72198/20088/72662/223499/78581/109075/12181/20826/74778/66409/57808/106338/207932/224481/67045/213773/19896/72462/54127/224823/66229/53414/18519/72515/225348/67936/67222/20044/68731/20227/67674/107071/59028/54364/107094/66583/18150/18572/227522/27176/227656/227715/66489/75422/98956/66181/271842/66580/16418/228889/66481/69639/68147/229937/70285/100019/230082/20104/12578/72181/20116/69072/67205/67025/69902/50912/74035/75416/17089/65961/94061/100503670/117109/71990/208144/100609/74097/72171/213895/110109/60321/68077/13871/20085/27998/14113/20055/67973/20068/101867/67920/67276/234344/234374/78651/214627/67619/72544/215194/71989/382137/77591/102462/67115/66446/72554/245474/76130/69499
##            Count
## GO:0022613   288
## GO:0034660   293
## GO:0034470   248
## GO:0042254   212
## GO:0006364   153
```

---
## Gene ontology and functional testing

Network plots can be generated from any enrichResult object using the **enrichplot** package.

We measure similarities between the various significant gene sets and group them accordingly. The **showCategory** argument specifies how many top gene ontology hits to show.

```r
library(enrichplot)
GO_result_plot <- pairwise_termsim(GO_result)
emapplot(GO_result_plot, showCategory = 20)
```

![](ChIPseq_In_Bioconductor3_files/figure-html/unnamed-chunk-12-1.png)

---
## Gene ontology and functional testing

In addition to gene ontology, we can test our gene list against custom gene sets that we import as gmt files using the clusterProfiler **enricher** function. Similar to the enrichGO function, this will produce an enrichResult object that can be used for visualization.

Here we will get gene sets from MSigDB using the **msigdbr** package. We can use this package to pull out specific gene sets, which can be seen [at this link](http://software.broadinstitute.org/gsea/msigdb/). We can also run the **msigdbr_collections** function to see the categories and subcategory codes that will be used for accessing the gene sets.

```r
library(msigdbr)
msigdbr_collections()
```

```
## # A tibble: 23 × 3
##    gs_cat gs_subcat         num_genesets
##    <chr>  <chr>                    <int>
##  1 C1     ""                         278
##  2 C2     "CGP"                     3368
##  3 C2     "CP"                        29
##  4 C2     "CP:BIOCARTA"              292
##  5 C2     "CP:KEGG"                  186
##  6 C2     "CP:PID"                   196
##  7 C2     "CP:REACTOME"             1604
##  8 C2     "CP:WIKIPATHWAYS"          615
##  9 C3     "MIR:MIR_Legacy"           221
## 10 C3     "MIR:MIRDB"               2377
## # … with 13 more rows
```

---

## Gene ontology and functional testing

From the data frame on the previous slide we can identify the category/subcategory we want, and use those in the **msigdbr** function. Here we will use "H" to access the Hallmark gene sets, and in the end we need to get a data frame where the first column contains the name of the gene sets and the second column contains the gene ID.

```r
library(msigdbr)
msig_t2g <- msigdbr(species = "Mus musculus", category = "H", subcategory = NULL)
msig_t2g <- msig_t2g[, colnames(msig_t2g) %in% c("gs_name", "entrez_gene")]
msig_t2g[1:3, ]
```

```
## # A tibble: 3 × 2
##   gs_name               entrez_gene
##   <chr>                       <int>
## 1 HALLMARK_ADIPOGENESIS       11303
## 2 HALLMARK_ADIPOGENESIS       74610
## 3 HALLMARK_ADIPOGENESIS       52538
```

---
## Gene ontology and functional testing

We then run the gene set enrichment, using the term to gene mapping we created as the **TERM2GENE** argument in the enricher function.

```r
hallmark <- enricher(gene = genesWithPeakInTSS, universe = allGeneIDs, TERM2GENE = msig_t2g)
hallmark_df <- data.frame(hallmark)
hallmark_df[1:3, ]
```

```
##                                              ID             Description
## HALLMARK_MYC_TARGETS_V1 HALLMARK_MYC_TARGETS_V1 HALLMARK_MYC_TARGETS_V1
## HALLMARK_E2F_TARGETS       HALLMARK_E2F_TARGETS    HALLMARK_E2F_TARGETS
## HALLMARK_MYC_TARGETS_V2 HALLMARK_MYC_TARGETS_V2 HALLMARK_MYC_TARGETS_V2
##                         GeneRatio  BgRatio       pvalue     p.adjust
## HALLMARK_MYC_TARGETS_V1  154/1425 200/4396 4.710378e-40 2.355189e-38
## HALLMARK_E2F_TARGETS     126/1425 200/4396 1.049344e-19 2.623360e-18
## HALLMARK_MYC_TARGETS_V2   50/1425  58/4396 1.878537e-17 3.130896e-16
##                               qvalue
## HALLMARK_MYC_TARGETS_V1 1.784985e-38
## HALLMARK_E2F_TARGETS    1.988231e-18
## HALLMARK_MYC_TARGETS_V2 2.372889e-16
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       geneID
## HALLMARK_MYC_TARGETS_V1 15510/18393/66942/17219/51810/107508/15182/56200/67332/66094/103136/18674/53605/12461/73192/12567/56351/18813/67465/18972/12464/103573/18148/14694/12330/52530/22333/27041/13681/12261/22627/110809/18102/18673/16211/26446/22123/16647/20382/18263/22630/26443/13665/13877/218236/105148/28126/30877/11792/15381/12465/18458/55944/14375/15382/66409/14852/22195/19385/20462/28185/21454/19826/67097/20383/109801/108121/27756/66973/15516/225027/20641/225363/27366/70769/11757/56086/20823/229279/100042807/110074/99138/67134/67204/26444/12428/15191/12462/68436/51788/13684/19359/20104/230484/107995/51797/230721/26445/74326/20810/230908/19934/69719/14208/13204/11991/11837/19988/19384/22384/17220/27979/381760/53379/56150/66870/12468/68011/23983/12785/12034/68988/20103/21849/107686/53607/14113/57296/50995/19899/16828/68981/27050/20133/13690/233870/19935/21781/17218/24015/14719/17463/101943/85305/26441/12858/56403/23918/67115/66446/18655/108062/53381/66235
## HALLMARK_E2F_TARGETS                                                                                                                                                                                   18393/17219/67196/98386/66403/12534/21917/53605/103468/12567/18813/18972/21379/15574/103573/15366/19360/20425/20877/78304/22059/114714/110809/18102/67684/21973/103733/16647/15374/20382/21877/11799/67241/18949/13665/218210/66197/69716/19355/70472/14375/26934/105988/12419/14852/19385/20462/22042/106618/15361/22154/18969/70769/20833/12236/110074/20638/67134/12531/69270/17865/110750/69639/70099/140917/51788/27354/12578/230484/50927/73804/107995/51797/11637/19891/100336/16765/21335/17089/19687/11991/50883/22256/12704/19384/19718/74097/17220/18861/12190/100710/71777/70699/14056/56150/66979/17768/68298/19362/77891/16881/12447/18971/108961/13626/233726/22390/18817/66422/109145/17218/68278/69724/17535/13433/16201/67629/12649/66131/56403/20842/12530/17350/19139/24061/53381
## HALLMARK_MYC_TARGETS_V2                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  27993/15510/66942/15278/338359/20509/70683/12567/18813/64934/21379/18148/18432/69071/18673/28126/30877/75751/53414/72515/21453/13340/107071/59028/107094/226169/70769/20322/13537/67134/68493/53608/69902/20810/100608/56361/22256/15277/110109/27998/101612/71974/67973/18817/66422/17218/66590/66164/235036/76130
##                         Count
## HALLMARK_MYC_TARGETS_V1   154
## HALLMARK_E2F_TARGETS      126
## HALLMARK_MYC_TARGETS_V2    50
```

---
## Gene ontology and functional testing

We learned about the **goseq** package during the RNAseq course, which is another functional annotation package similar to clusterProfiler,

Here we perform the same enrichment test for the MSigDB Hallmark gene sets.

For goseq, we need a named vector of all genes (the universe) with 1s or 0s representing whether a gene had peak in TSS or not. 
We can turn a logical vector into 1 for TRUE and 0 for FALSE simply using the **as.integer** function.

```r
allGenesForGOseq <- as.integer(allGeneIDs %in% genesWithPeakInTSS)
names(allGenesForGOseq) <- allGeneIDs
allGenesForGOseq[1:3]
```

```
## 100009600 100009609 100009614 
##         1         0         0
```

---
## Gene ontology and functional testing

First we must construct a **nullp** data.frame for use within **goseq** using the **nullp** function and supplying our named vector, genome to be used and gene identifier used.

The **nullp** function attempts to correct for gene length biases we may see in gene set testing. i.e. a longer gene may have more chance to have a peak within it.

```r
library(goseq)
pwf = nullp(allGenesForGOseq, "mm10", "knownGene", plot.fit = FALSE)
```

```
## Can't find mm10/knownGene length data in genLenDataBase...
```

```
## Found the annotation package, TxDb.Mmusculus.UCSC.mm10.knownGene
```

```
## Trying to get the gene lengths from it.
```

---
## Gene ontology and functional testing

We can use the same term to gene mapping we used for clusterProfiler (though it must be converted from a tibble to data frame for goseq) to run the gene set enrichment test.

```r
Myc_hallMarks <- goseq(pwf, "mm10", "knownGene", gene2cat = data.frame(msig_t2g))

Myc_hallMarks[1:3, ]
```

```
##                   category over_represented_pvalue under_represented_pvalue
## 32 HALLMARK_MYC_TARGETS_V1            6.653087e-41                        1
## 13    HALLMARK_E2F_TARGETS            4.136764e-20                        1
## 33 HALLMARK_MYC_TARGETS_V2            1.166240e-17                        1
##    numDEInCat numInCat
## 32        154      200
## 13        126      200
## 33         50       58
```

---

## GREAT for GO and functional testing

In addition to a standard enrichment tests, methods have been implemented specifically for ChIPseq. Many of these tools aim to incorporate peaks distal to genes into their enrichment testing such as the popular [GREAT](http://bejerano.stanford.edu/great/public/html/splash.php) toolset.

Incorporating distal peaks by rules such as nearest gene results in some genes having a higher chance of being selected and hence some gene sets as a whole having a higher chance of having its members selected.

[GREAT](http://bejerano.stanford.edu/great/public/html/splash.php) defines regulatory regions for each individual gene and compares the proportion of peaks mapping to a gene set's regulatory regions to the proportion of the genome occupied by gene set's regulatory regions.

i.e. If a gene set's regulatory regions account for 1% of the genome then one might expect 1% of peaks to overlap these regions by chance.

---
rGREAT - R interface to GREAT server
=========================================================

We can use the GREAT Bioconductor interface available in the rGREAT package.

```r
library(rGREAT)
```

```
## 
## ------------------
## Note: On Aug 19 2019 GREAT released version 4 where it supports `hg38`
## genome and removes some ontologies such pathways. `submitGreatJob()`
## still takes `hg19` as default. `hg38` can be specified by the `species
## = 'hg38'` argument. To use the older versions such as 3.0.0, specify as
## `submitGreatJob(..., version = '3.0.0')`.
## ------------------
```
---
rGREAT for GO and functional testing
=========================================================

To submit jobs we can use our GRanges of Myc peaks and specify a genome with the **submitGreatJob** function.

This function returns a GreatJob object containing a reference to our results on the GREAT server. To review the categories of results available we can use the availableCategories function on our GreatJob object.

```r
great_Job <- submitGreatJob(macsPeaks_GR, species = "mm10", version = "3.0.0", request_interval = 1)
availableCategories(great_Job)
```

```
## [1] "GO"                               "Phenotype Data and Human Disease"
## [3] "Pathway Data"                     "Gene Expression"                 
## [5] "Regulatory Motifs"                "Gene Families"
```
---
rGREAT for GO and functional testing
=========================================================

The results table can be retrieved using the getEnrichmentTables function and specifying which tables we wish to review.

Here we retrieve the results tables for the "Regulatory Motifs" gene sets which contains 2 separate database results.

```r
great_ResultTable = getEnrichmentTables(great_Job, category = "Regulatory Motifs")
names(great_ResultTable)
```

```
## [1] "MSigDB Predicted Promoter Motifs" "MSigDB miRNA Motifs"
```
---
rGREAT for GO and functional testin
=========================================================

Now we can review the enrichment of our genes with Myc peaks in their TSS for the "MSigDB Predicted Promoter Motifs" gene sets.

```r
msigProMotifs <- great_ResultTable[["MSigDB Predicted Promoter Motifs"]]
msigProMotifs[1:4, ]
```

```
##                     ID
## 1   SCGGAAGY_V$ELK1_02
## 2     GGGCGGR_V$SP1_Q6
## 3      CACGTG_V$MYC_Q2
## 4 RCGCANGCGY_V$NRF1_Q6
##                                                                              name
## 1                Motif SCGGAAGY matches ELK1: ELK1, member of ETS oncogene family
## 2                             Motif GGGCGGR matches SP1: Sp1 transcription factor
## 3 Motif CACGTG matches MYC: v-myc myelocytomatosis viral oncogene homolog (avian)
## 4                     Motif RCGCANGCGY matches NRF1: nuclear respiratory factor 1
##   Binom_Genome_Fraction Binom_Expected Binom_Observed_Region_Hits
## 1            0.05912876       814.6169                       1586
## 2            0.20132920      2773.7130                       3782
## 3            0.07555623      1040.9380                       1638
## 4            0.05163166       711.3293                       1195
##   Binom_Fold_Enrichment Binom_Region_Set_Coverage Binom_Raw_PValue
## 1              1.946928                0.11511940    1.738355e-136
## 2              1.363515                0.27451550     1.823467e-94
## 3              1.573580                0.11889380     1.110696e-71
## 4              1.679953                0.08673877     2.182812e-65
##   Binom_Adjp_BH Hyper_Total_Genes Hyper_Expected Hyper_Observed_Gene_Hits
## 1 1.069088e-133              1069       502.9698                      806
## 2  5.607161e-92              2689      1265.1880                     1812
## 3  2.276927e-69               941       442.7452                      695
## 4  3.356073e-63               822       386.7551                      606
##   Hyper_Fold_Enrichment Hyper_Gene_Set_Coverage Hyper_Term_Gene_Coverage
## 1              1.602482              0.07821446                0.7539757
## 2              1.432198              0.17583700                0.6738565
## 3              1.569752              0.06744299                0.7385760
## 4              1.566883              0.05880640                0.7372263
##   Hyper_Raw_PValue Hyper_Adjp_BH
## 1     1.739538e-83  5.349079e-81
## 2    5.523247e-114 3.396797e-111
## 3     1.936613e-65  3.970057e-63
## 4     1.614352e-56  2.482066e-54
```

---
class: inverse, center, middle

# Motifs

---

## Motifs

A common practice in transcription factor ChIPseq is to investigate the motifs enriched under peaks.

De novo motif enrichment can be performed in R/Bioconductor but this can be very time consuming. Here we will use the MEME-ChIP suite available online to identify de novo motifs.

MEME-ChIP requires a FASTA file of sequences under peaks as input so we extract this using the **BSgenome** package.

---

## Extracting sequences under regions

First we need to load the BSgenome object for the genome we are working on, UCSC's mm10 build for the mouse genome, **BSgenome.Mmusculus.UCSC.mm10**.

```r
library(BSgenome)
library(BSgenome.Mmusculus.UCSC.mm10)
BSgenome.Mmusculus.UCSC.mm10
```

```
## Mouse genome:
## # organism: Mus musculus (Mouse)
## # genome: mm10
## # provider: UCSC
## # release date: Dec. 2011
## # 66 sequences:
## #   chr1                 chr2                 chr3                
## #   chr4                 chr5                 chr6                
## #   chr7                 chr8                 chr9                
## #   chr10                chr11                chr12               
## #   chr13                chr14                chr15               
## #   ...                  ...                  ...                 
## #   chrUn_GL456372       chrUn_GL456378       chrUn_GL456379      
## #   chrUn_GL456381       chrUn_GL456382       chrUn_GL456383      
## #   chrUn_GL456385       chrUn_GL456387       chrUn_GL456389      
## #   chrUn_GL456390       chrUn_GL456392       chrUn_GL456393      
## #   chrUn_GL456394       chrUn_GL456396       chrUn_JH584304      
## # (use 'seqnames()' to see all the sequence names, use the '$' or '[[' operator
## # to access a given sequence)
```
---

## Extracting sequences under regions

The motif for the ChIP-ed transcription factor should near the summit of a peak. MEME-ChIP will trim our peaks to a common length internally if sequences are of different length.

It is best therefore to provide a peak set resized to a common length.

```r
macsSummits_GR <- GRanges(seqnames(macsPeaks_GR), IRanges(macsPeaks_GR$abs_summit,
    macsPeaks_GR$abs_summit), score = macsPeaks_GR$fold_enrichment)
macsSummits_GR <- resize(macsSummits_GR, 100, fix = "center")
```

---
## Extracting sequences under regions

We now have a GRanges, centred on the summit, highest point of signal for every peak.

```r
macsSummits_GR
```

```
## GRanges object with 13777 ranges and 1 metadata column:
##           seqnames              ranges strand |     score
##              <Rle>           <IRanges>  <Rle> | <numeric>
##       [1]     chr1     4785515-4785614      * |   5.23296
##       [2]     chr1     5083072-5083171      * |   4.26417
##       [3]     chr1     7397785-7397884      * |   9.84580
##       [4]     chr1     9545350-9545449      * |  11.85079
##       [5]     chr1     9700928-9701027      * |   5.14291
##       ...      ...                 ...    ... .       ...
##   [13773]     chrX 168674050-168674149      * |   5.21700
##   [13774]     chrX 169047705-169047804      * |   5.33257
##   [13775]     chrX 169320320-169320419      * |   4.67030
##   [13776]     chrY     1245745-1245844      * |   4.04550
##   [13777]     chrY     1286577-1286676      * |   5.00143
##   -------
##   seqinfo: 21 sequences from an unspecified genome; no seqlengths
```

---
## Extracting sequences under regions

Once we have recentered our peaks we can use the **getSeq** function with our GRanges of resized common peaks and the BSgenome object for mm10.

The **getSeq** function returns a *DNAStringSet* object containing sequences under peaks.

```r
peaksSequences <- getSeq(BSgenome.Mmusculus.UCSC.mm10, macsSummits_GR)
names(peaksSequences) <- paste0(seqnames(macsSummits_GR), ":", start(macsSummits_GR),
    "-", end(macsSummits_GR))

peaksSequences[1:2, ]
```

```
## DNAStringSet object of length 2:
##     width seq                                               names               
## [1]   100 CGTCCATCGCCAGAGTGACGCGG...CTGGGCTTCAGGTTGGCCAGGCT chr1:4785515-4785614
## [2]   100 CCGCCCTCAGCTGCGGTCACGTG...TGAGTGAGGCTCGCAACGTCTCC chr1:5083072-5083171
```
---
## Writing to FASTA file

The *writeXStringSet* function allows the user to write DNA/RNA/AA(aminoacid)StringSet objects out to file.

By default the *writeXStringSet* function writes the sequence information in FASTA format (as required for MEME-ChIP).

```r
writeXStringSet(peaksSequences, file = "mycMel_rep1.fa")
```

---
## MEME-ChIP

Now the file "mycMel_rep1.fa" contains sequences around the geometric center of peaks suitable for Motif analysis in MEME-ChIP.

In your own work you will typically run this from your own laptop with MEME installed locally but today we will upload our generated FASTA file to their [web portal](http://meme-suite.org/tools/meme-chip).

Follow instructions [here](http://meme-suite.org/doc/download.html) to install MEME localy.

Results files from MEME-ChIP can be found [here](http://rockefelleruniversity.github.io/myc_Meme_Example/meme-chip.html)

---
## Parsing back FIMO results

We can retrieve the locations of Myc motifs identified in MEME-ChIP from the FIMO output.

FIMO reports Myc motif locations as a GFF3 file which we should be able to vizualise in IGV. 
Sadly, this GFF file's naming conventions cause only a fraction of motifs to be reported.

---
## FIMO to R

Fortunately we can parse our motif's GFF file into R and address this using the **import** function in  the **rtracklayer** package.

```r
library(rtracklayer)
motifGFF <- import("~/Downloads/fimo.gff")
```

---
## FIMO to valid GFF3

We can give the sequences some more sensible names and export the GFF to file to visualise in IGV.

```r
motifGFF$Name <- paste0(seqnames(motifGFF), ":", start(motifGFF), "-", end(motifGFF))
motifGFF$ID <- paste0(seqnames(motifGFF), ":", start(motifGFF), "-", end(motifGFF))
export.gff3(motifGFF, con = "~/Downloads/fimoUpdated.gff")
```

---
## Scanning for known motifs

We saw previously we can scan sequences using some of the Biostrings functionality **matchPattern**.

Often with ChIPseq we may know the motif we are looking for or we can use a set of known motifs from a database such as a [JASPAR](http://jaspar.genereg.net).

We can access to JASPAR using the JASPAR2020 bioconductor library.

```r
library(JASPAR2020)
JASPAR2020
```

```
## An object of class "JASPAR2020"
## Slot "db":
## [1] "/usr/local/lib/R/host-site-library/JASPAR2020/extdata/JASPAR2020.sqlite"
```

---
## Get motifs from JASPAR with TFBStools

We can access the model for the our motif of interest using the **TFBSTools** package and its **getMatrixByName** function.

```r
library(TFBSTools)
pfm <- getMatrixByName(JASPAR2020, name = "MYC")
pfm
```

```
## An object of class PFMatrix
## ID: MA0147.3
## Name: MYC
## Matrix Class: Basic helix-loop-helix factors (bHLH)
## strand: +
## Tags: 
## $centrality_logp
## [1] "-296.164"
## 
## $family
## [1] "bHLH-ZIP factors"
## 
## $medline
## [1] "18555785"
## 
## $remap_tf_name
## [1] "MYC"
## 
## $source
## [1] "29126285"
## 
## $tax_group
## [1] "vertebrates"
## 
## $tfbs_shape_id
## [1] "140"
## 
## $type
## [1] "ChIP-seq"
## 
## $unibind
## [1] "1"
## 
## $collection
## [1] "CORE"
## 
## $species
##           9606 
## "Homo sapiens" 
## 
## $acc
## [1] "P01106"
## 
## Background: 
##    A    C    G    T 
## 0.25 0.25 0.25 0.25 
## Matrix: 
##   [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## A 1038 1019  438   44 4600   54  169   59   23   419   961   755
## C 1345 1287 3746 4840   59 4525  157  162   56  2801  1524  1785
## G 1649 1757  472   31  170   70 4578   63 4817  1198  1036  1348
## T  910  879  286   27  113  293   38 4658   46   524  1421  1054
```

---
## Motif scanning with motifmatchr

With this PWM we can use the **motifmatchr** package to scan our summits for the Myc motif and return the positions of the motifs.

We will need to provide our PWM, GRanges to scan within and BSGenome object to extract sequence from.

We also set the **out** paramter to positions for this instance.

```r
library(motifmatchr)
MycMotifs <- matchMotifs(pfm, macsSummits_GR, BSgenome.Mmusculus.UCSC.mm10, out = "positions")
MycMotifs
```

```
## GRangesList object of length 1:
## [[1]]
## GRanges object with 1600 ranges and 1 metadata column:
##          seqnames              ranges strand |     score
##             <Rle>           <IRanges>  <Rle> | <numeric>
##      [1]     chr1     9545433-9545444      - |   14.0312
##      [2]     chr1   13151280-13151291      + |   15.7259
##      [3]     chr1   13293444-13293455      + |   14.5201
##      [4]     chr1   16665901-16665912      - |   14.6571
##      [5]     chr1   43561269-43561280      + |   14.2468
##      ...      ...                 ...    ... .       ...
##   [1596]     chrX 155338504-155338515      + |   14.8300
##   [1597]     chrX 155338504-155338515      - |   14.8266
##   [1598]     chrX 160610006-160610017      + |   12.6818
##   [1599]     chrX 166454533-166454544      + |   13.6368
##   [1600]     chrX 166454533-166454544      - |   13.6368
##   -------
##   seqinfo: 21 sequences from an unspecified genome; no seqlengths
```

---
## Exporting motif matches

We can export the Myc motif positions within peaks for use later in IGV or for metaplot vizualisation.

```r
export.bed(MycMotifs[[1]], con = "MycMotifs.bed")
```

---
## Time for an exercise!

Exercise on ChIPseq data can be found [here](../../exercises/exercises/chipseq_part3_exercise.html)

---
## Answers to exercise

Answers can be found [here](../../exercises/answers/chipseq_part3_answers.html)