Codes in paper

Bar graph of the relationship between the amount of DNA methylation valleys and H3K4me3 and H3K27me3


DNA methylation valleys (DNA methylation valleys, DMV) are a wide range of hypomethylated regions, which are related to the activation of histone marker H3K4me3 and the inhibition of histone marker H3K27me3. In metastatic castration-resistant prostate cancer (mCRPC, metastatic castration-resistant prostate cancer), as shown in the figure, the greater the number of DMVs in the sample, the higher the signal ratio of H3K4me3 to H3K27me3.

Sample-level log2(OR) calculated from the number of DMVs that overlap H3K4me3 vs. H3K27me3 sites. Lower values favor H3K27me3 and higher values favor H3K4me3. Lower: Sample-level count of DMVs, order matching upper panel.

Code explanation

Read data
#Loop through HMRs per sample and assign them to groups
output <- data.frame()
samples2use <- sample_ids_wgbs
valleydf <- read.delim(file=fn_valley,sep='\t',header=F,stringsAsFactors=F)
colnames(valleydf) <- c('chr','start','end','sample')
valley <- makeGRangesFromDataFrame(valleydf,keep.extra.columns=T)
Calculate the number of H3K4me3 and H3K27me3 peaks in the methylation valley, and do fisher test to calculate the OR of their relationship with the methylation valley
for(i in 1:length(samples2use)) {
    sample <- samples2use[i]
    valleyi <- valley[valley$sample==sample]

    valleywidth <- sum(width(valleyi))
    output[sample,'Sample'] <- sample  #Sample name 
    output[sample,'width'] <- valleywidth  # Total width of methylated valley
    output[sample,'count'] <- length(valleyi)  # Number of methylated valleys

    print(paste(i, output[sample,'count']))

    in_h3k4me3 <- countOverlaps(tracks[['pca100_H3K4me3']],valleyi,ignore.strand=T)>0
    in_h3k27me3 <- countOverlaps(tracks[['pca100_H3K27me3']],valleyi,ignore.strand=T)>0
    in_cpg_island <- countOverlaps(tracks[['cpg_island']],valleyi,ignore.strand=T)>0

    count_h3k4me3_valley <- sum(in_h3k4me3)
    count_h3k27me3_valley <- sum(in_h3k27me3)
    count_h3k4me3_novalley <- sum(!in_h3k4me3)
    count_h3k27me3_novalley <- sum(!in_h3k27me3)

    tab <- rbind(c(count_h3k4me3_valley,count_h3k4me3_novalley),
                 c(count_h3k27me3_valley,count_h3k27me3_novalley)) # Generate quadruple table
    ft <- fisher.test(tab)  # Fisher test, whether the signal values of h3k4me3 and h3k27me3 are related to the methylation valley

    output[sample,'OR'] <- log2(ft$estimate)
    output[sample,'p'] <- ft$p.value

output <- output[rownames(output) %in% sample_ids_wgbs,]
output$fdr <- p.adjust(output$p,method='fdr')
output$logfdr <- -log10(output$fdr)

lvl <- output[order(output$count),'Sample']
#lvl <- output[order(output$width),'Sample']
output$Sample <- factor(output$Sample,levels=lvl)
Drawing according to the output matrix, ggarrange layout drawing
#Set up data frames for plotting
df1 <-$Sample,count=output$OR,SS=output$fdr<=0.05)
df3 <-$Sample,count=output$width/HG38_SIZE*100)
df4 <-$Sample,count=output$count)

cols <- brewer.pal(n=4, name='Dark2')
plist <- list()
plist[[1]] <- ggplot(df1,aes(x=sample,y=count))+theme_classic()+geom_bar(stat='identity')+
    theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5))+ylab('Log2 Odds Ratio')+
plist[[3]] <- ggplot(df3,aes(x=sample,y=count))+theme_classic()+geom_bar(stat='identity',fill=cols[2])+
    ylab('DMV %Genome')  # This picture is not shown in the article in the end
plist[[2]] <- ggplot(df4,aes(x=sample,y=count))+theme_classic()+geom_bar(stat='identity',fill=cols[3])+
    ylab('DMV Count')



Zhao S G , Chen W S , Li H , et al. The DNA methylation landscape of advanced prostate cancer[J]. Nature Genetics, 2020.

Original code

Leave a Reply

Your email address will not be published. Required fields are marked *