# A way to test for enrichment of differentially expressed genes in a genomic location

I have an experiment where I expect a certain genomic location to influence gene expression levels of nearby genes. I have data for expression levels (Agilent 4x44 microarrays, Drosophila) in two groups - one where I expect expression to be affected and the other wild-type and I would like to run a test for overrepresentation of differentially expressed genes in a genomic location.

My main problem is that I couldn't find a package (R/bioconductor) that would do it out of the box easily, so if you know about such a package, please let me know. In the meantime, this is what I figured out: I would run a sliding window over the whole genome and simply count number of differentially expressed genes in each window - this should tell me where I have the most differentially expressed genes in the genome. However, it will be dependent on gene density, so to obtain some sort of background distribution, I would run permutations of the samples (or p values), say, 1000 times, and check how often I am likely to find this number of windows with that number of differentially expressed genes compared to the observed numbers. Does this sound right?

I should add that while I know the location that would mess up things, I cannot exclude that any other genomic region would not be affected as well. So I have to test the whole genome.

• Why would you look at the entire genome when you have a control for the region that you know you are interested in? – learner Apr 15 '13 at 16:20
• The main reason I think is that my modification may also affect distant genomic regions (some trans effect, I imagine if I hit a transcription factor, many genes elsewhere may be messed up). So while I expect local effect, I cannot exclude distal ones... Also but less importantly, I think it would make my statement about local effect much stronger if I could show that only that region is affected in the whole genome. – yotiao Apr 15 '13 at 21:22
• Your idea is good. You might save a little effort and see what happens when you normalize your DE count by gene density; i.e. $\frac{DE_i}{GD_i}$ before attempting the permutation testing. – learner Apr 16 '13 at 12:20

This depends on what you mean by a genomic location. For each cytoband this would be rather straight forward to do. Roughly:

1) Get the cytoband locations for all genes. These are stored in the organism specific packages, e.g., org.Dm.eg.db, and are named as 'MAP' . You might need the chiptype specific annotation package to map between the probe identifiers and the genes first.

2) Once you have the cytoband annotations for the genes, you can then test each cytoband separately with the functionality offered by, e.g., the topGO package. There is a section with the heading 'Predefined list of interesting genes' in the vignette of the topGO package that shortly shows how to do this is a similar case.

For the smoothing approach you have thought of, it might be worth correcting the counts with the actual number of genes in any predefined window, taking into account that not all genes might be present on the chip. The exact gene locations are available in the organism specific annotation package (the same as above). Some difficulties might arise, since certain locations probably have a gene in both strands, so you just need to decide how to count them.

The cytoband based approach is available in, e.g., Chipster software (see the manual entry at http://chipster.csc.fi/manual/stat-hyperG-cytoband.html), and the source code for the analysis is available at https://github.com/chipster/chipster/blob/master/src/main/modules/microarray/R-2.12/stat-hyperG-cytoband.R, which might in some details, if you decide to use the cytobands.

• The only slight problem with this approach is that there are no chiptype specific annotation packages for Agilent drosophila chips. But it can be overcome, as the Agilent output chip itself maps the probes to genomic positions. Or I can do it the hard way and go through biomaRt and annotate the genes (not probes) to cytobands on my own. – yotiao Oct 14 '13 at 12:44
• You're right, of course, Bioconductor project does not distribute the annotations for the Drosophila array. Just in case you decide to try the cytoband approach, I have produced an annotation package using the refseq information from Agilent, and it is shared on my Google Drive at docs.google.com/file/d/0B5F_KFI2_sBKZzZXR3ZsV2pvQ1k/…. I just updated it today. – JTT Oct 14 '13 at 18:25

If you can get or create genesets by genomic location then geneset enrichment analysis (GSEA) would be one way to do this.

I don't know of any pre-defined genesets for Drosophila (not an organism I work with) but there is a lot of bioconductor/R support for both building them and carrying out the analysis.