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
5.2 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.7k views
ADD COMMENTlink
modified 3.5 years ago by Biostar ♦♦ 20 • written 5.2 years ago by pmanga20

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

ADD REPLYlink written 3.2 years ago by ek0
0
5.2 years ago by
Devon Ryan94k
Freiburg, Germany
Devon Ryan94k 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 5.2 years ago by Devon Ryan94k

Thanks..that helped!

ADD REPLYlink written 5.2 years ago by pmanga20
0
5.2 years ago by
dariober11k
WCIP | Glasgow | UK
dariober11k 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 5.2 years ago by dariober11k

Thanks a lot that is very helpful!!

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

Content
Help
Access

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