Very low probability densities

Question

I'm generating a density plot using the ggplot2 library:

ggplot(df_chr20, aes(x=start1)) + geom_density() + theme_bw() 
ggplot(df_chr21, aes(x=start1)) + geom_density() + theme_bw() 

The data frame for df_chr20 has 1573 rows, and I wanted to plot the probability density of genomic regions as x = start1. The same for df_chr21, with 1271 rows. What I see in the plot is what I expect, but is there a reason why the y-axis probablilty densities are so low? The area under the curve is surely not equal to 1.

I expected the high localised high density in chr 21, and both the plots have similar scales in terms of y-axis values. But since they are so low, how can I interpret this plot?

I found some information here, but my values are much lower than the examples.

Any help would be very kind!

Answer 1

It is not like your conventional frequency which adds up to one. The density is low because the width of your bins is huge, and the number of observations you have is low. From this by Wickham, the basic kernel is:

where K is the kernel and h is the bandwidth (Scott, 1992b)

If you don't specify it, by default h will revert to nrd0:

library(rtracklayer)
gr_narrowPeak <- import("https://www.encodeproject.org/files/ENCFF072FBH/@@download/ENCFF072FBH.bed.gz",format = "BED",extraCols = extraCols_narrowPeak)
gr_chr21 = gr_narrowPeak[seqnames(gr_narrowPeak)=="chr21"]

da = as.data.frame(gr_chr21[rep(1:length(gr_chr21),gr_chr21$signalValue)])
dens = density(da$start)
dens$bw
[1] 1133889

plot(dens)

ggplot(da,aes(x=start))+geom_density()

You can decrease h such that you get something smaller (but at the cost of looking coarse). Important thing is it doesn't sum up to one.

If you would like to explain it easier like a histogram, consider converting it into counts, like explained in this post:

binW = 50000

ggplot(da,aes(x=start))+
geom_density(aes(y = ..density..*(nrow(da)*binW)),bw=50000) +
geom_histogram(binwidth=binW,fill="steelblue",alpha=0.5) + 
theme_bw()

ggplot(da,aes(x=start))+
geom_histogram(aes(y=..count../nrow(da)),binwidth=50000) + 
geom_density(aes(y = ..density..*binW),bw=binW,fill="steelblue",alpha=0.5) +
theme_bw()