Question: Code for drawing a histogram with binned real gene transcription counts from an RNA-Seq experiment data and fit to negative binomial/normal/poission distribution curve
0
gravatar for pmanga
4.4 years ago by
pmanga20
United States
pmanga20 wrote:

Hi,

I have gene transcription counts from an RNA seq experiment and I wanted to draw a histogram with these  and fit it to negative binomial/normal/Poisson distribution curve. I want to get some comparisons like these http://www.plosone.org/article/slideshow.action?uri=info:doi/10.1371/journal.pone.0016327&imageURI=info:doi/10.1371/journal.pone.0016327.g001

Can someone help me with the code for this in R?

R • 3.2k views
ADD COMMENTlink modified 2.7 years ago by Biostar ♦♦ 20 • written 4.4 years ago by pmanga20

updated link: http://dx.doi.org/10.1371/journal.pone.0016327

ADD REPLYlink written 2.4 years ago by ek0
0
gravatar for Devon Ryan
4.4 years ago by
Devon Ryan89k
Freiburg, Germany
Devon Ryan89k wrote:

Use ggplot2 with the geom_histogram() function and just set binwidth to whatever you like (this is likely how that figure was made). For fitting, see the R functions glm.nb() and such.

ADD COMMENTlink written 4.4 years ago by Devon Ryan89k

Thanks..that helped!

ADD REPLYlink written 4.4 years ago by pmanga20
0
gravatar for dariober
4.4 years ago by
dariober10.0k
WCIP | Glasgow | UK
dariober10.0k wrote:

I think the general strategy should be that you find the set of parameters for you chosen model that best fits the data, for example via maximum likelihood. Then for plotting you plug these parameters in your model to get the corresponding values at a grid of points.

In a R these could be done like this:

require(MASS)
set.seed(2)
dat<- rnbinom(1000, mu = 5, size = 3) ## Obs data points
par.nb<- fitdistr(dat, "negative binomial")
par.pois<- fitdistr(dat, "Poisson")

at<- min(dat):max(dat)
p.pois<- dpois(at, par.pois$estimate)
p.nb<- dnbinom(at, size= par.nb$estimate[1], mu= par.nb$estimate[2])

hist(dat, breaks= 25, freq= FALSE, col= 'grey80',
    border= 'white', ylim= c(0, 0.2))
lines(names(table(dat)), table(dat)/length(dat), lwd= 2, col= 'grey30')
lines(at, p.pois, lwd= 2, col= 'blue')
lines(at, p.nb, lwd= 2, col= 'red')
legend('topright', lwd=2, cex= 0.8, col= c('grey30', 'red', 'blue'),
    legend= c('Obs.', 'Fitted neg. binom.', 'Fitted Poisson'))

ADD COMMENTlink written 4.4 years ago by dariober10.0k

Thanks a lot that is very helpful!!

ADD REPLYlink written 4.4 years ago by pmanga20
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1881 users visited in the last hour