In CharlesJB/metageneVignettes: Vignette for the metagene paper

Enrichment of the H3K27ac histone mark in active and bivalent transcription start sites (TSS)

Introduction

The H3K27ac histone mark is normally enriched at active transcriptional start sites (TSS) compared to poised TSS (Roadmap Epigenomics Consortium et al.,2015). metagene has been used to robustly test for the difference in enrichment of H3K27ac between these two groups of promoter regions in A549 cells.

1. Load metagene results data

The results have been saved in a list that can be loaded and used by similaRpeak. The list contains 2 elements: a data.frame and a Friedman test results. Only the data.frame is needed.

suppressMessages(library(packrat))
packrat::on()
load("results_S4.RData")
# Extract data.frame 
DF<-results[["DF"]]

The data.frame content can be used to recreate the metagene graph.

suppressMessages(library(ggplot2))
p <- ggplot(DF, aes(x=position, y=value, ymin=qinf, ymax=qsup))
p + geom_ribbon(aes(fill=group), alpha=0.3) + 
    geom_line(aes(color=group), size=1, litype=1, bg="transparent") +
    ylab("Reads per million (RPM)")

Figure S6. H3K27ac TssA versus H3K27ac TssBiv

2. Load similaRpeak package

The similaRpeak package must be loaded.

library(similaRpeak)

3. Calculate metrics using similaRpeak

The metrics are calculated using the similarity function which takes as arguments the two profiles vectors and the threshold values.

# H3K27ac_TssA vector
H3K27ac_TssA <-DF$value[DF$group == "H3K27ac_TssA"]
# H3K27ac_TssBiv vector
H3K27ac_TssBiv<-DF$value[DF$group == "H3K27ac_TssBiv"]
# Calculate metrics
metrics <- similarity(H3K27ac_TssA,
                            H3K27ac_TssBiv,
                            ratioAreaThreshold=0.1,
                            ratioMaxMaxThreshold=0.01,
                            ratioIntersectThreshold=0.1,
                            ratioNormalizedIntersectThreshold=0.1,
                            diffPosMaxThresholdMinValue=0.01,
                            diffPosMaxThresholdMaxDiff=100,
                            diffPosMaxTolerance=0.01)
metrics

The RATIO_AREA and RATIO_MAX_MAX reflect the difference of amplitude between both profiles. The small RATIO_INTERSECT is caused by the small values of the H3K27ac_TssBiv profile.

The RATIO_NORMALIZED_INTERSECT is quite high. However, the graph of the normalized histone profiles shows that this value is caused by the relatively flat profiles.

plot(H3K27ac_TssA*length(H3K27ac_TssA)/sum(H3K27ac_TssA, na.rm=TRUE), xaxt = "n", 
        type="l", col="blue", xlab="", ylab="", ylim=c(0, 9))
par(new=TRUE)
plot(H3K27ac_TssBiv*length(H3K27ac_TssBiv)/sum(H3K27ac_TssBiv, na.rm=TRUE), 
        type="l", col="darkgreen", xlab="position", 
        ylab="Normalized Coverage (Coverage/Mean Coverage)", 
        ylim=c(0, 9), xaxt = "n")
legend("topright", c("H3K27ac_TssA","H3K27ac_TssBiv"), cex=1.2, 
        col=c("blue","darkgreen"), lty=1)

Figure S7. Normalized H3K27ac TssA versus normalized H3K27ac TssBiv

In conclusion, the profile of H3K27ac histone mark in active transcriptional start sites (TSS) compared to poised TSS are not sharing any of the tested patterns.