Question: Color intensity plot peak positions in R from BED file
0
j.mccarter10 wrote:

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??

Basically like this track on hg19

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

ucsc bed chip-seq ngs R • 2.1k views
modified 5.7 years ago by dariober11k • written 5.7 years ago by j.mccarter10
2
dariober11k wrote:

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)
``` This is what I'm looking for - Ill have a go! Thanks

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

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

0
Alex Reynolds31k wrote:

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