#
Enrichment of MHBs in open chromatin regions

In [ ]:

    library("dplyr")
    library("ggplot2")

## An example using colon MHBs¶

In [ ]:

    # run computeCpgCov in mHapSuite
    system("java -jar mHapSuite-2.0-jar-with-dependencies.jar computeCpgCov \
            -bigwig colon_MM.bw -cpgPath hg19_CpG.gz -openChromatin colon_open.bed -tag colon_MHB \
            -bedPath colon_MHB.bed -missingDataAsZero")

In [11]:

    # the output file shows mean methylation of each CpG sites and whether it's covered by open chromatin regions
    result=read.table("colon_MHB.txt")
    names(result) <- c("chr", "beta", "Open")
    head(result)

A data.frame: 6 × 3
| chr | beta | Open
--- | --- | --- | ---
| <chr> | <dbl> | <int>
1 | chr1:540838-540839 | 0.6142857 | 0
2 | chr1:540841-540842 | 0.5692308 | 0
3 | chr1:540850-540851 | 0.5820895 | 0
4 | chr1:540854-540855 | 0.5909091 | 0
5 | chr1:540866-540867 | 0.5468750 | 0
6 | chr1:540868-540869 | 0.5781250 | 0

In [21]:

    # assign CpG into groups by mean methylation-levels
    result <- result[order(result$beta, decreasing = FALSE),]
    result$Cluster <- paste0("cluster_", rep(1:10, each = trunc(nrow(result) / 10) + 1))[1:nrow(result)]

    result$total=as.numeric("1")
    head(result);table(result$Cluster)
    res = as.data.frame(result[c("total","Open","Cluster")] %>% group_by(Cluster) %>% summarise_each(funs=sum))
    rownames(res) =res$Cluster;res$Cluster = NULL
    cluster_res = as.data.frame(result[c("beta","Cluster")] %>% group_by(Cluster) %>% summarise_each(funs=median))
    rownames(cluster_res) = cluster_res$Cluster;cluster_res$Cluster=NULL
    res = cbind(res,cluster_res)
    res$Freq = 100*round(res$Open/res$total,4)
    res$cluster=rownames(res)
    head(res);dim(res)
    colon_MHB=res[,c(3,4,5)]
    head(colon_MHB)

A data.frame: 6 × 5
| chr | beta | Open | Cluster | total
--- | --- | --- | --- | --- | ---
| <chr> | <dbl> | <int> | <chr> | <dbl>
30850 | chr16:103150-103151 | 0.00282486 | 1 | cluster_1 | 1
30851 | chr16:103156-103157 | 0.00284900 | 1 | cluster_1 | 1
30853 | chr16:103186-103187 | 0.00285714 | 1 | cluster_1 | 1
30855 | chr16:103203-103204 | 0.00289017 | 1 | cluster_1 | 1
30852 | chr16:103174-103175 | 0.00292398 | 1 | cluster_1 | 1
30854 | chr16:103195-103196 | 0.00292398 | 1 | cluster_1 | 1

     cluster_1 cluster_10  cluster_2  cluster_3  cluster_4  cluster_5  cluster_6
         12646      12641      12646      12646      12646      12646      12646
     cluster_7  cluster_8  cluster_9
         12646      12646      12646

A data.frame: 6 × 5
| total | Open | beta | Freq | cluster
--- | --- | --- | --- | --- | ---
| <dbl> | <int> | <dbl> | <dbl> | <chr>
cluster_1 | 12646 | 12428 | 0.0335383 | 98.28 | cluster_1
cluster_10 | 12641 | 9752 | 0.7291667 | 77.15 | cluster_10
cluster_2 | 12646 | 12351 | 0.1585033 | 97.67 | cluster_2
cluster_3 | 12646 | 12373 | 0.2166667 | 97.84 | cluster_3
cluster_4 | 12646 | 12243 | 0.2761194 | 96.81 | cluster_4
cluster_5 | 12646 | 12082 | 0.3562366 | 95.54 | cluster_5

1. 10
2. 5

A data.frame: 6 × 3
| beta | Freq | cluster
--- | --- | --- | ---
| <dbl> | <dbl> | <chr>
cluster_1 | 0.0335383 | 98.28 | cluster_1
cluster_10 | 0.7291667 | 77.15 | cluster_10
cluster_2 | 0.1585033 | 97.67 | cluster_2
cluster_3 | 0.2166667 | 97.84 | cluster_3
cluster_4 | 0.2761194 | 96.81 | cluster_4
cluster_5 | 0.3562366 | 95.54 | cluster_5

In [26]:

    # perform similar analysis for UMR, LMR, PMD and HMR
    input=read.table("colon_merge.txt",header = T);head(input);table(input$type)

A data.frame: 6 × 4
| beta | Freq | cluster | type
--- | --- | --- | --- | ---
| <dbl> | <dbl> | <chr> | <chr>
1 | 0.4705882 | 25.65 | cluster_1 | HMR
2 | 0.9802632 | 9.99 | cluster_10 | HMR
3 | 0.7368421 | 24.35 | cluster_2 | HMR
4 | 0.8248175 | 19.88 | cluster_3 | HMR
5 | 0.8715084 | 17.82 | cluster_4 | HMR
6 | 0.9010989 | 16.66 | cluster_5 | HMR

    HMR LMR MHB PMD UMR
     10  10  10  10  10

In [68]:

    options(repr.plot.width = 8, repr.plot.height=6)
    ggplot(input,aes(x=beta,y=Freq,color=type))+geom_smooth(se=F,size=1)+scale_y_continuous(expand=c(0,0))+
        scale_x_continuous(limits=c(0,1))+ggtitle("Colon")+
        scale_color_manual(values = c('#56A902','#2e409a','#fcbe32','#e97f02','#f15c5c'),name="")+
        labs(x="Mean methylation",y="The Percent of CpG sites \noverlapped by Open chromatin regions (%)")+theme_classic()+
        theme(plot.title = element_text(hjust=0.5,vjust=0.5),legend.text=element_text(size=13),axis.title=element_text(size=16),
            axis.text = element_text(size = 13))+coord_cartesian(ylim = c(0, 100))+
        geom_vline(aes(xintercept = 0.25),colour="black",linetype="dashed")

    `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

## Average profile across all tissue types¶

In [60]:

    total=read.table("merge.txt",header = T)
    total$tissue=(sapply(strsplit(total$tissue, "\\."), function(z) z[1]))
    total$type=sapply(strsplit(total$tissue, "\\_"), function(z) z[2])
    total$tissue=sapply(strsplit(total$tissue, "\\_"), function(z) z[1])
    head(total);table(total$tissue,total$type)

A data.frame: 6 × 5
| beta | Freq | cluster | tissue | type
--- | --- | --- | --- | --- | ---
| <dbl> | <dbl> | <chr> | <chr> | <chr>
1 | 0.3300000 | 9.25 | cluster_1 | Adipose | HMR
2 | 1.0000000 | 2.18 | cluster_10 | Adipose | HMR
3 | 0.7272727 | 9.72 | cluster_2 | Adipose | HMR
4 | 0.8333333 | 5.54 | cluster_3 | Adipose | HMR
5 | 0.8857143 | 3.99 | cluster_4 | Adipose | HMR
6 | 0.9180328 | 3.12 | cluster_5 | Adipose | HMR

                HMR LMR MHB PMD UMR
      Adipose    10  10  10  10  10
      Adrenal    10  10  10  10  10
      B-cell     10  10  10  10  10
      Breast     10  10  10  10  10
      colon      10  10  10  10  10
      esophagus  10  10  10  10  10
      heart      10  10  10  10  10
      liver      10  10  10  10  10
      lung       10  10  10  10  10
      ovary      10  10  10  10  10
      pancreas   10  10  10  10  10
      placenta   10  10  10  10  10
      spleen     10  10  10  10  10
      stomach    10  10  10  10  10
      T-cell     10  10  10  10  10
      thymus     10  10  10  10  10
      thyroid    10  10  10  10  10

In [64]:

    total_input=aggregate(total[1:2], by = list(total$cluster, total$type), mean);head(total_input)
    dim(total_input);table(total_input$`Group.2`)

A data.frame: 6 × 4
| Group.1 | Group.2 | beta | Freq
--- | --- | --- | --- | ---
| <chr> | <chr> | <dbl> | <dbl>
1 | cluster_1 | HMR | 0.4646884 | 16.720588
2 | cluster_10 | HMR | 0.9917105 | 5.696471
3 | cluster_2 | HMR | 0.7557492 | 13.751765
4 | cluster_3 | HMR | 0.8430695 | 10.236471
5 | cluster_4 | HMR | 0.8873174 | 8.662941
6 | cluster_5 | HMR | 0.9144253 | 7.845294

1. 50
2. 4

    HMR LMR MHB PMD UMR
     10  10  10  10  10

In [67]:

    options(repr.plot.width = 8, repr.plot.height=6)
    ggplot(total_input,aes(x=beta,y=Freq,color=Group.2))+geom_smooth(se=F,size=1)+scale_y_continuous(expand=c(0,0))+
        scale_x_continuous(limits=c(0,1))+ggtitle("17 tissues")+
        scale_color_manual(values = c('#56A902','#2e409a','#fcbe32','#e97f02','#f15c5c'),name="")+
        labs(x="Mean methylation",y="The Percent of CpG sites \noverlapped by Open chromatin regions (%)")+theme_classic()+
        theme(plot.title = element_text(hjust=0.5,vjust=0.5),legend.text=element_text(size=13),axis.title=element_text(size=16),
            axis.text = element_text(size = 13))+coord_cartesian(ylim = c(0, 100))+
        geom_vline(aes(xintercept = 0.25),colour="black",linetype="dashed")

    `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
    Warning message:
    “Removed 1 rows containing non-finite values (`stat_smooth()`).”

In [ ]: