DESeq/EdgeR on counts of reads supporting variant

[brentp]

Is there precedent for and/or is it sensible to use DESeq or EdgeR to model counts of reads supporting variants in order to get significance?

For example, suppose I have 10 cases and 10 controls in a 20 sample VCF, can I extract the counts of alternate alleles (after filtering) into a counts matrix to these tools without violating their assumptions?

Some have recommended that I do a contingency table on number of cases / controls with / without the variant to get significance but that is discarding a lot of information.

thanks,
-Brent

EDIT: I suppose the difficulty is adjusting for depth at each site.

[dpryan]

I would try a beta-binomial test (you can do this with the vglm function in the VGAM package). The beta-binomial distribution is similar to a regular binomial distribution, except the rates are derived from beta distributions (i.e, they're not a single number but a distribution that describes a population).

The format would be something like:

library(VGAM)
groups <- factor(c(rep("case",10),rep("control", 10)), level=c("control","case"))
reference_count <- c(rbetabinom(10,seq(100,200,10),0.8), rbetabinom(10,seq(100,200,10),0.7)) #some example data
alt_count <- rep(seq(100,190,10),2) - reference_count
d <- data.frame(ref=reference_count, alt=alt_count, grp=groups)
fit <- vglm(cbind(ref,alt)~grp, betabinomial, d)
lrtest(fit)
#One can get the coefficients with
logit(coef(fit), inverse=T)
#But keep in mind only the intercepts will be correct, so
logit(sum(coef(fit)[c(1,3)]), inverse=T) #The actual "case" mu

One can also fit the shape parameters, if you prefer. BTW, this is more or less how BiSeq works (the problem ends up being similar), though it uses a different package.

[dpryan]

Keep in mind that the distribution can be a bit tricky to fit, so you sometimes have to play around with parameters (don't blindly accept a fit without looking at the resulting deviations. You can also use the "binomial.ab" to fit the normal paramaterization of the distribution, though I suspect you're interested in changes in mu.

[brentp]

Thanks Devon, I'll look into that more. I knew that BiSeq used the betaregression, but didn't think about the details.

Some others have also suggested that I look into plink/seq or pseq.

[dpryan]

It's just a question of whether the data follows a binomial distribution or is over-dispersed. In the former case logistic regression (as done by Plink or in R) would be appropriate, but in the latter case not.