#
Enrichment analysis of tissue-specific genes in MHBs

In [23]:

    library(rGREAT)
    library("org.Hs.eg.db")
    library(pheatmap)
    options(warn=-1)

## An example using colon MHBs¶

In [6]:

    gr=read.table("colon_MHB.bed")
    colnames(gr)=c("chr","start","end")
    gr=GRanges(gr)
    start(gr) = start(gr) + 1
    head(gr)

    GRanges object with 6 ranges and 0 metadata columns:
          seqnames        ranges strand
             <Rle>     <IRanges>  <Rle>
      [1]     chr1 540839-540886      *
      [2]     chr1 564471-564534      *
      [3]     chr1 566287-566317      *
      [4]     chr1 567312-567438      *
      [5]     chr1 567681-567786      *
      [6]     chr1 568529-568570      *
      -------
      seqinfo: 24 sequences from an unspecified genome; no seqlengths

In [7]:

    # loading tissue-specific genes
    TS = readRDS("tissueSpecificGenes.RDS")
    ls(TS)

1. 'Adipose Tissue'
2. 'Adrenal Gland'
3. 'Appendix'
4. 'Bone Marrow'
5. 'Breast'
6. 'Cerebral Cortex'
7. 'Cervix, uterine'
8. 'Colon'
9. 'Duodenum'
10. 'Endometrium'
11. 'Esophagus'
12. 'Fallopian Tube'
13. 'Gallbladder'
14. 'Heart Muscle'
15. 'Kidney'
16. 'Liver'
17. 'Lung'
18. 'Lymph Node'
19. 'Ovary'
20. 'Pancreas'
21. 'Placenta'
22. 'Prostate'
23. 'Rectum'
24. 'Salivary Gland'
25. 'Seminal Vesicle'
26. 'Skeletal Muscle'
27. 'Skin'
28. 'Small Intestine'
29. 'Smooth Muscle'
30. 'Spleen'
31. 'Stomach'
32. 'Testis'
33. 'Thyroid Gland'
34. 'Tonsil'
35. 'Urinary Bladder'

In [9]:

    # colon-specific genes
    example = as.character(TS[[8]][,1])
    head(example)

1. 'PIGY'
2. 'NOS2'
3. 'CEACAM7'
4. 'NOX1'
5. 'CLCA1'
6. 'CLCA4'

In [16]:

    SR = unique(mapIds(org.Hs.eg.db, example, 'ENTREZID', 'SYMBOL'))
    SR = SR[!is.na(SR)]
    gs = list(example = SR)

    res = great(gr, gs, "hg19", verbose = FALSE)
    tb = getEnrichmentTable(res)

    'select()' returned 1:1 mapping between keys and columns

    * extended_tss is already cached, directly use it.

In [17]:

    # enrichment of colon-specific genes in colon MHBs with high significance
    tb

A data.frame: 1 × 12
| id | genome_fraction | observed_region_hits | fold_enrichment | p_value | p_adjust | mean_tss_dist | observed_gene_hits | gene_set_size | fold_enrichment_hyper | p_value_hyper | p_adjust_hyper
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
| <chr> | <dbl> | <int> | <dbl> | <dbl> | <dbl> | <dbl> | <int> | <int> | <dbl> | <dbl> | <dbl>
1 | example | 0.009976902 | 208 | 1.438598 | 3.721069e-07 | 3.721069e-07 | 48792 | 70 | 143 | 1.109497 | 0.1394002 | 0.1394002

In [18]:

    test=tb[,c('id','p_adjust')]
    test$rnkPV=rownames(test)
    log=test$p_adjust>0.05
    test$rnkPV[log]<-NA
    test

A data.frame: 1 × 3
| id | p_adjust | rnkPV
--- | --- | --- | ---
| <chr> | <dbl> | <chr>
1 | example | 3.721069e-07 | 1

## Show enrichment by ranks of significance¶

In [54]:

    result=read.table("Exp.txt",sep="\t",header = T, row.names = 1)

In [55]:

    head(result)

A data.frame: 6 × 17
| Adipose | Adrenal | B.cell | Breast | colon | esophagus | heart | liver | lung | ovary | pancreas | placenta | spleen | stomach | T.cell | thymus | thyroid
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
| <int> | <int> | <lgl> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <int> | <lgl> | <int> | <int>
Adipose Tissue | 1 | NA | NA | 5 | NA | NA | 2 | NA | NA | NA | NA | 3 | NA | NA | NA | NA | NA
Adrenal Gland | NA | 1 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 2 | NA
Breast | NA | NA | NA | 1 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA
Colon | NA | NA | NA | NA | 1 | NA | NA | NA | NA | NA | 2 | NA | NA | NA | NA | NA | NA
Esophagus | NA | NA | NA | 3 | NA | 1 | NA | NA | NA | NA | NA | 2 | NA | 2 | NA | NA | NA
Heart Muscle | NA | NA | NA | NA | NA | NA | 1 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA

In [56]:

    # mask non-significant ones
    result_input=result
    result_input[is.na(result_input)] <- 14
    head(result_input)

A data.frame: 6 × 17
| Adipose | Adrenal | B.cell | Breast | colon | esophagus | heart | liver | lung | ovary | pancreas | placenta | spleen | stomach | T.cell | thymus | thyroid
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl>
Adipose Tissue | 1 | 14 | 14 | 5 | 14 | 14 | 2 | 14 | 14 | 14 | 14 | 3 | 14 | 14 | 14 | 14 | 14
Adrenal Gland | 14 | 1 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 2 | 14
Breast | 14 | 14 | 14 | 1 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14
Colon | 14 | 14 | 14 | 14 | 1 | 14 | 14 | 14 | 14 | 14 | 2 | 14 | 14 | 14 | 14 | 14 | 14
Esophagus | 14 | 14 | 14 | 3 | 14 | 1 | 14 | 14 | 14 | 14 | 14 | 2 | 14 | 2 | 14 | 14 | 14
Heart Muscle | 14 | 14 | 14 | 14 | 14 | 14 | 1 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14 | 14

## Heatmap¶

In [57]:

    options(repr.plot.width = 8, repr.plot.height=6, repr.plot.res = 150)
    pheatmap(result_input,
        show_rownames=T,main="Tissue-specific gene",
        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,14),legend_labels = c("1","14"),
        scale = "none",border_color = "grey90", display_numbers = matrix(ifelse(result_input > 1, "", "*"), nrow = nrow(result_input)),
        color = colorRampPalette(colors = c("#F71E35", "#FCD271","#fcebb6","white"))(50))