Color intensity plot peak positions in R from BED file

Question

[j.mccarter1]

Hi Folks

I wanted to make a plot in R and Im struggling with it!

I was hoping to use a BED file to plot blocks of equal height but varying length (peak positions) along the x axis (genomic position) the reason I want equal height is so I could plot using color intensity based on the fold enrichment values from the MACs peak caller. Exactly like the UCSC genome browser "dense" view. Can somebody give me an example of this??

Many thanks in advance

Basically like this track on hg19

http://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr21%3A33031597-33041570&hgsid=202635987_ULlPiiRDsWNKFm0lTj8C7b9KG0IS

Answer 1

[dariober]

This is a pretty low level way of achieving it, if I correctly understand:

## Dummy data in ~bed format
n<- 10
chrom<- rep('chr1', n)
start<- seq(1, 1000, length.out= n)
end<- start + rpois(n, 50)
logfc<- rgamma(n, 2)

## Prepare colours based on logFC
pal<- colorRampPalette(c('blue', 'red'))
cols<- paste0(pal(n)[as.numeric(cut(logfc, breaks = n))], 99)

## Plot
plot(x= start, y= rep(0, n), type= 'n', bty= 'n', yaxt= 'n',
    ylab= '', xlab= 'Position', xlim= range(c(start, end)), ylim= c(0, 1))
rect(xleft= start, xright= end, ybottom= 0, ytop= 0.1, border= NA, col= cols)
text(x= start, y= 0.12, labels= sprintf("%.2f", logfc), adj= c(0,0))

## Legend
nlg<- 5
lgcols<- paste0(pal(nlg)[1:nlg], 99)
lg<- round(seq(min(logfc), max(logfc), length.out= nlg), 1)
legend('topleft', pch= 15, bty= 'n', col= lgcols, legend= lg, pt.cex= 2, cex= 1)

[j.mccarter1]

This is what I'm looking for - Ill have a go! Thanks

[j.mccarter1]

Ok so I got it to work with my data and it looks good - thank you - However I was hoping to get a color legend embedded also. Do you have any Idea how I would do that? I have looked at plot3D

[dariober]

See edit with legend. it's not very nice because the numbers are not nice & round, but maybe it gives you a hint...

Answer 2

[Alex Reynolds]

Maybe just export the UCSC browser shot directly from an existing session, or use a fetch tool, if automation is needed.

[j.mccarter1]

We are using a custom genome (not on UCSC) and as far as i know we cant add a custom genome to it. We like our previous outputs using R which is why we were hoping to use it again