#
Region-set enrichment analysis using LOLA

In [ ]:

    library("LOLA")
    library("GenomicRanges")
    library("Hmisc")
    library("pheatmap")

## An example using MHBs from adrenal gland¶

In [ ]:

    MHB="BED/Adrenal_MHB.bed"

    # MHB regions
    background="BED/MHB_Cat.bed"
    mhb=readBed(MHB)

    # Universe_Sets Background
    universe_set=readBed(background)

    # Overlaped Regions
    states=loadRegionDB("BED/hg19_test")

    # Runing LOLA
    locResults = runLOLA(mhb, universe_set, states, cores=5)

In [58]:

    constant = 17
    locResults$qValue = (10^(-locResults$pValueLog)) * (constant / locResults$rnkPV)
    head(locResults);dim(locResults)

A data.table: 6 × 23
userSet | dbSet | collection | pValueLog | oddsRatio | support | rnkPV | rnkOR | rnkSup | maxRnk | ⋯ | d | description | cellType | tissue | antibody | treatment | dataSource | filename | size | qValue
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
<int> | <int> | <chr> | <dbl> | <dbl> | <int> | <int> | <int> | <int> | <int> | ⋯ | <int> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <dbl> | <dbl>
1 | 2 | ucsc_example | 210.191253 | 2.185228 | 3529 | 1 | 1 | 1 | 1 | ⋯ | 69405 | ucsc_example | NA | NA | NA | NA | NA | Adrenal_H3K27ac.bed | 158566 | 1.094450e-209
1 | 10 | ucsc_example | 58.564747 | 1.537888 | 2500 | 2 | 2 | 11 | 11 | ⋯ | 75085 | ucsc_example | NA | NA | NA | NA | NA | Ovary_H3K27ac.bed | 194882 | 2.315642e-58
1 | 1 | ucsc_example | 10.963068 | 1.188053 | 2707 | 3 | 3 | 10 | 10 | ⋯ | 66551 | ucsc_example | NA | NA | NA | NA | NA | Adipose_H3K27ac.bed | 123493 | 6.169644e-11
1 | 13 | ucsc_example | 10.092572 | 1.174367 | 3167 | 4 | 4 | 4 | 4 | ⋯ | 59554 | ucsc_example | NA | NA | NA | NA | NA | Spleen_H3K27ac.bed | 130465 | 3.434132e-10
1 | 14 | ucsc_example | 8.211183 | 1.153098 | 3454 | 5 | 5 | 2 | 5 | ⋯ | 54817 | ucsc_example | NA | NA | NA | NA | NA | Stomach_H3K27ac.bed | 135084 | 2.090721e-08
1 | 17 | ucsc_example | 6.673084 | 1.135736 | 3124 | 6 | 6 | 5 | 6 | ⋯ | 59348 | ucsc_example | NA | NA | NA | NA | NA | Thyroid_H3K27ac.bed | 142299 | 6.014693e-07

1. 17
2. 23

In [56]:

    # mask non-significant results
    log=locResults$qValue>0.05
    locResults$rnkPV[log]<-NA
    result=locResults[,c('rnkPV','filename')];result

A data.table: 17 × 2
rnkPV | filename
--- | ---
<int> | <chr>
1 | Adrenal_H3K27ac.bed
2 | Ovary_H3K27ac.bed
3 | Adipose_H3K27ac.bed
4 | Spleen_H3K27ac.bed
5 | Stomach_H3K27ac.bed
6 | Thyroid_H3K27ac.bed
7 | Lung_H3K27ac.bed
8 | Heart_H3K27ac.bed
9 | Liver_H3K27ac.bed
NA | Esophagus_H3K27ac.bed
NA | Pancreas_H3K27ac.bed
NA | Colon_H3K27ac.bed
NA | Thymus_H3K27ac.bed
NA | Breast_H3K27ac.bed
NA | Placenta_H3K27ac.bed
NA | Bcell_H3K27ac.bed
NA | Tcell_H3K27ac.bed

## Enrichment of MHBs in H3K27ac ChIP-seq peaks¶

In [9]:

    Result=read.table("BED/Result_H3K27ac.txt",sep="\t",header = T)
    colnames(Result)=(sapply(strsplit(colnames(Result), "\\."), function(z) z[1]))
    colnames(Result)=capitalize(colnames(Result))
    rownames(Result)=Result$Filename
    Result=Result[,-1]
    rownames(Result)=capitalize(rownames(Result))
    colnames(Result)=(sapply(strsplit(colnames(Result), "\\_"), function(z) z[1]))
    head(Result)

A data.frame: 6 × 17
| Adipose | Adrenal | B | Breast | Colon | Esophagus | Heart | Liver | Lung | Ovary | Pancreas | Placenta | Spleen | Stomach | T | Thymus | Thyroid
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
| <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int>
Adipose | 2 | 3 | 9 | 4 | 17 | 9 | 1 | 17 | 6 | 5 | 17 | 17 | 8 | 17 | 7 | 17 | 17
Adrenal gland | 5 | 1 | 14 | 17 | 17 | 17 | 2 | 17 | 13 | 3 | 17 | 17 | 6 | 17 | 15 | 17 | 17
B-cell | 13 | 17 | 2 | 17 | 2 | 17 | 17 | 17 | 8 | 17 | 17 | 17 | 2 | 17 | 2 | 2 | 17
Breast | 17 | 17 | 4 | 1 | 9 | 6 | 17 | 17 | 11 | 17 | 17 | 17 | 17 | 17 | 5 | 17 | 17
Colon | 6 | 17 | 11 | 17 | 1 | 10 | 17 | 17 | 3 | 17 | 17 | 17 | 7 | 3 | 6 | 17 | 17
Esophagus | 9 | 17 | 10 | 2 | 5 | 1 | 17 | 17 | 7 | 17 | 17 | 17 | 17 | 5 | 11 | 17 | 17

## Heatmap¶

In [34]:

    options(repr.plot.width = 8, repr.plot.height=6)
    pheatmap(Result,
        show_rownames=T,main="H3K27ac Enrichment rnkPV",
        cluster_cols= F,angle_col=45,fontsize_row=11,cellwidth=21,
        cluster_rows=F,cellheight=18,fontsize_col=12,
        how_colnames = T,
        legend_breaks = c(1,17),legend_labels = c("1","17"),
        scale = "none",border_color = "grey90", display_numbers = matrix(ifelse(Result > 1, "", "*"),
                                                                         nrow = nrow(Result)),
        color = colorRampPalette(colors = c("#F71E35", "#FCD271","#fcebb6","white"))(50))

In [ ]: