I am using edgeR for analyzing shRNA data. So far I got the DE pvalue per shRNA. Now I want to summarize the results per gene rather than per shRNA. What is the "official" way to do this? 

There are ad hoc things such as for example: pick significant shRNAs and then summarize per gene, taking care of conflicting differentials from probes.

I was hoping to do something more elegant, that might involve a linear model which per genes is something like ~probe+condition, this way I can read off the DE with the probe effect removed from the condition coefficient in a linear model fit per gene. edgeR includes a GLM feature which I was hoping will help with this, but I could not figure out how, because different genes have different number of shRNAs so a table with rows = "genes" and columns = "conditions repeated for every possible probe (up to the max number of shRNA probes in a gene)" ends up with a lot of NAs and edgeR dose not seem to be capable of dealing with this. This is also funny because edgeR thinks of columns as samples and ends normalizing the table... It is a mess.

Any ideas? 

There are several options here. The simplest is to take the shRNA with the highest average abundance, and to use that as a representative of the entire gene. However, this doesn't really take advantage of the additional evidence provided by the other shRNAs of that gene.

The second option is to use Simes' method (check out the combineTests function in the csaw package). The null hypothesis here is that none of the shRNAs within a gene are DE. Simes' method will compute a combined p-value for each gene, based on the p-values for all shRNAs within that gene. This combined p-value represents the evidence against the aforementioned null, and incorporates information from all shRNAs in the gene. Genes with low combined p-values are of interest as they will contain at least one DE shRNA.

However, it may not be appropriate in your case, as strong DE for one shRNA can result in a low combined p-value for the entire gene. From what I understand, you are looking for genes that have consistent DE across several shRNAs. That brings me to the third option, which is to treat all shRNAs in a gene as a "gene set" (i.e., shRNAs are the "genes" in the set, and the gene containing the shRNAs is the "gene set"). You can then use methods like ROAST to do a self-contained gene set test. The set will be statistically significant only when the majority of the shRNAs are DE, which might be what you want. Check out ?roast and ?roast.DGEList for more information.

edgeR and other GLM- and LM-based differential expression packages are designed to fit the same design model to every gene, so you will not be able to use them directly in this way.

However, edgeR has a function called spliceVariants that fits a separate model for each gene to test the interaction between exon and group while estimating dispersions at the gene level. You could treat your shRNAs like "exons", but you don't want an interaction term, you want an additive mode. So you could write a similar function that fits your desired model to each gene to test for condition effects using shRNA ID as a blocking factor.

Aaron has described several sensible options for summarizing shRNAs per gene. Generally however we advise people not to average shRNAs if at all possible. This is because some shRNAs are much better than others, and averaging a good one with several bad ones isn't a gain. I would rather that you got a list of DE shRNAs, and then examine the genes they belong to. If you want the "official" line, this is it.

The edgeR glm features are not designed to fit models across rows, so they can't be used for this purpose. 

Summarizing shRNA data per gene is a difficult problem with no good generic solution to my knowledge. There are two papers that made important contributions towards solving it:

• Use very many shRNAs per gene ("Ultracomplex Libraries"), trust in the statistics of large numbers, and learn from the data which shRNAs are effective. See e.g. "A Systematic Mammalian Genetic Interaction Map Reveals Pathways Underlying Ricin Susceptibility" by  Mike Bassik et al.  http://dx.doi.org/10.1016/j.cell.2013.01.030
• Tune your algorithm using a 'ground truth', see e.g. "Measuring error rates in genomic perturbation screens: gold standards for human functional genomics" by Traver Hart and Jason Moffat's lab, DOI: 10.15252/msb.20145216