Testing for chromosome over-representation
Sometimes I am working on some data and notice certain biases, say differentially expressed genes appearing to originate more often from a chromosome. Or a factor binding more often to a class of transcripts. In these situations I tend to turn to Fisher’s exact test. Here I will put an example of what I do.
Get some data
For the sake of simplicity I will use data from the Pasilla1 dataset and run differential gene expression analysis with DESeq22 following the vignette’s instructions. The data.table package is used because I like it’s speed and syntax, specially for sub-setting. Since I am still learning all the ins and outs of it, I mix data.frame code with it. Whatever works.
library("data.table")
library("pasilla")
library("Biobase")
data("pasillaGenes")
library("DESeq2")
library("biomaRt")While loading the pasillaGenes R was throwing an error at me:
Warning: namespace ‘DESeq’ is not available and has been replaced
It turns out that DESeq needed to be installed to load the data. No idea why. So let’s get us some results:
countData <- counts(pasillaGenes)
colData <- pData(pasillaGenes)[,c("condition","type")]
dds <- DESeqDataSetFromMatrix(
countData = countData,
colData = colData,
design = ~ condition
)
dds$condition <- relevel(dds$condition, "untreated")
dds <- DESeq(dds)Add chromosome information
Now the fun starts. I created a data.table with the results plus gene IDs, And will now added some extra information with biomaRt. I have only added chromosome and biotype, but the amount of information one can add is large. I usually also add gene symbol and description. Very useful to have a quick idea of the function of a particular gene. Adding the full gene location might be useful to later on subset data and easily create a bed file.
res <- as.data.frame(results(dds))
res$ensembl_gene_id <- rownames(results(dds))
res <- data.table(res)
ensembl = useMart("ensembl", dataset = "dmelanogaster_gene_ensembl")
genemap <- getBM( attributes = c("ensembl_gene_id", "chromosome_name", "gene_biotype"),
filters = "ensembl_gene_id",
values = res$ensembl_gene_id,
mart = ensembl)
idx <- match( res$ensembl_gene_id, genemap$ensembl_gene_id )
res$chromosome <- genemap$chromosome_name[ idx ]
res$gene_biotype <- genemap$gene_biotype[ idx ]Now let’s have a look at down-regulated genes with padj < 0.1.
res[padj < 0.1 & log2FoldChange < 0]## baseMean log2FoldChange lfcSE stat pvalue
## 1: 85.00693 -0.7067901 0.2196386 -3.217969 1.291017e-03
## 2: 273.13010 -0.3878740 0.1360678 -2.850593 4.363776e-03
## 3: 1257.88508 -0.3896927 0.1321601 -2.948640 3.191754e-03
## 4: 280.32056 -0.4471012 0.1417352 -3.154483 1.607826e-03
## 5: 141.36961 -1.3409377 0.2076951 -6.456280 1.073079e-10
## ---
## 404: 998.56567 -0.4335869 0.1357833 -3.193228 1.406921e-03
## 405: 16.42232 -0.6973109 0.2276632 -3.062906 2.191992e-03
## 406: 373.24468 -0.4988888 0.1441100 -3.461860 5.364557e-04
## 407: 2891.60114 -1.6134844 0.1706094 -9.457183 3.163516e-21
## 408: 973.77088 -0.8325356 0.1332384 -6.248464 4.145082e-10
## padj ensembl_gene_id chromosome gene_biotype
## 1: 2.007727e-02 FBgn0000079 2R protein_coding
## 2: 5.230444e-02 FBgn0000244 3R protein_coding
## 3: 4.105429e-02 FBgn0000256 2L protein_coding
## 4: 2.389391e-02 FBgn0000286 2L protein_coding
## 5: 1.045641e-08 FBgn0000406 2L protein_coding
## ---
## 404: 2.158653e-02 FBgn0261545 2R protein_coding
## 405: 3.079189e-02 FBgn0261546 NA NA
## 406: 9.950931e-03 FBgn0261547 3L protein_coding
## 407: 1.106943e-18 FBgn0261552 3R protein_coding
## 408: 3.710330e-08 FBgn0261560 2L protein_codingMost of these appear to be coding genes. Not a great surprise. But… Oh! There are a lot of genes located in chromosome 2L. Is this a coincidence? Let us test it.
Test for over-representation
Chromosome 3L
chrom="3L"
all <- res$ensembl_gene_id
hits <- res[padj < 0.1 & log2FoldChange < 0]$ensembl_gene_id
hits_in_chr <- length(res[chromosome == chrom & padj < 0.1 & log2FoldChange < 0]$ensembl_gene_id)
genes_in_chr <- length(res[chromosome == chrom]$ensembl_gene_id)
hits_total <- length(hits)
genes_total <- length(all)Expected number of hits in chromosome 3L would be 71.8440912 and we have 61. Is this significant? I will use Fisher’s exact test, as applied here and here for very similar problems.
Firstly I construct a table/matrix with the events and from that calculated Fisher’s Exact test, or the probability of having more genes in chromosome 3L than expected by chance. There are more details here.
mat <- matrix(
c(
hits_in_chr,
genes_in_chr-hits_in_chr,
hits_total-hits_in_chr,
genes_total-hits_total-genes_in_chr+hits_in_chr
),
nrow=2,
ncol=2)
mat## [,1] [,2]
## [1,] 61 347
## [2,] 2487 11575fr <- fisher.test(mat, alternative="greater")
fr##
## Fisher's Exact Test for Count Data
##
## data: mat
## p-value = 0.935
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
## 0.6400274 Inf
## sample estimates:
## odds ratio
## 0.8181856So there is a p-value of 0.9349979 and we can reject the hypothesis that there is enrichment for or genes in 3L. My eyes were seeing patterns where this is none. Also, since I was testing only for over-representation, or enrichment, the option alternative="greater" was used in the test. Other options are available.
For all chromosomes
Is this this the case for any of the other chromosome? I will construct a function, and then loop over the chromosomes.
The function chrEnrichment is minimal and it only tests for down-regulated gens, but could be easily extended to add other arguments. It could also be used to test for biases in gene biotype.
chrEnrichment <- function(chr, df){
# Test for chromosome overpresentation in DESe2 results
all <- df$ensembl_gene_id
hits <- df[padj < 0.1 & log2FoldChange < 0]$ensembl_gene_id
hits_in_chr <- length(res[chromosome == chr & padj < 0.1 & log2FoldChange < 0]$ensembl_gene_id)
genes_in_chr <- length(res[chromosome == chr]$ensembl_gene_id)
hits_total <- length(hits)
genes_total <- length(all)
mat <- matrix(
c(
hits_in_chr,
genes_in_chr-hits_in_chr,
hits_total-hits_in_chr,
genes_total-hits_total-genes_in_chr+hits_in_chr
),
nrow=2,
ncol=2)
fr <- fisher.test(mat, alternative="greater")
df <- data.frame(
chromosome = chromosome,
observed = hits_in_chr,
expected = round(hits_total * genes_in_chr / genes_total, 1),
odds.ratio= fr$estimate[["odds ratio"]],
pvalue = fr$p.value
)
return(df)
}
chromosomes <- unique(na.omit(res[!like(chromosome,"Zv9")])$chromosome)
l_res <- list()
for (chromosome in chromosomes){
l_res[[chromosome]] <- chrEnrichment(chromosome, res)
}
table_fisher <- rbindlist(l_res)
table_fisher[,FDR:=p.adjust(pvalue, method='bonferroni'),] # add FDR## chromosome observed expected odds.ratio pvalue FDR
## 1: 2R 91 78.1 1.2197190 0.05874683 0.4112278
## 2: 3L 61 71.8 0.8181856 0.93499788 1.0000000
## 3: 2L 86 70.2 1.2958104 0.02276513 0.1593559
## 4: 3R 95 90.2 1.0719388 0.29753414 1.0000000
## 5: X 46 58.6 0.7523146 0.97261113 1.0000000
## 6: 4 0 2.3 0.0000000 1.00000000 1.0000000
## 7: YHet 0 0.4 0.0000000 1.00000000 1.0000000setkey(table_fisher, FDR) # will also sort by pvalueSince I am testing for all chromosomes, I have also calculated the adjusted p-value (bonferroni) to be on the safe side.
And what does Fisher say?
table_fisher## chromosome observed expected odds.ratio pvalue FDR
## 1: 2L 86 70.2 1.2958104 0.02276513 0.1593559
## 2: 2R 91 78.1 1.2197190 0.05874683 0.4112278
## 3: 3L 61 71.8 0.8181856 0.93499788 1.0000000
## 4: 3R 95 90.2 1.0719388 0.29753414 1.0000000
## 5: X 46 58.6 0.7523146 0.97261113 1.0000000
## 6: 4 0 2.3 0.0000000 1.00000000 1.0000000
## 7: YHet 0 0.4 0.0000000 1.00000000 1.0000000That we have a pretty much random distribution of down-regulated genes in the fly chromosomes.
References
sessionInfo()## R version 3.1.3 (2015-03-09)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.2 LTS
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=de_DE.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=de_DE.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel methods stats graphics grDevices utils
## [8] datasets base
##
## other attached packages:
## [1] biomaRt_2.22.0 DESeq2_1.6.3
## [3] RcppArmadillo_0.4.650.1.1 Rcpp_0.11.5
## [5] GenomicRanges_1.18.4 GenomeInfoDb_1.2.4
## [7] IRanges_2.0.1 S4Vectors_0.4.0
## [9] Biobase_2.26.0 BiocGenerics_0.12.1
## [11] pasilla_0.5.1 data.table_1.9.4
## [13] knitr_1.9
##
## loaded via a namespace (and not attached):
## [1] acepack_1.3-3.3 annotate_1.44.0 AnnotationDbi_1.28.1
## [4] base64enc_0.1-2 BatchJobs_1.5 BBmisc_1.9
## [7] BiocParallel_1.0.3 bitops_1.0-6 brew_1.0-6
## [10] checkmate_1.5.1 chron_2.3-45 cluster_2.0.1
## [13] codetools_0.2-11 colorspace_1.2-6 DBI_0.3.1
## [16] DESeq_1.18.0 digest_0.6.8 evaluate_0.5.5
## [19] fail_1.2 foreach_1.4.2 foreign_0.8-63
## [22] formatR_1.0 Formula_1.2-0 genefilter_1.48.1
## [25] geneplotter_1.44.0 ggplot2_1.0.1 grid_3.1.3
## [28] gtable_0.1.2 Hmisc_3.15-0 iterators_1.0.7
## [31] lattice_0.20-30 latticeExtra_0.6-26 locfit_1.5-9.1
## [34] MASS_7.3-39 munsell_0.4.2 nnet_7.3-9
## [37] plyr_1.8.1 proto_0.3-10 RColorBrewer_1.1-2
## [40] RCurl_1.95-4.5 reshape2_1.4.1 rpart_4.1-9
## [43] RSQLite_1.0.0 scales_0.2.4 sendmailR_1.2-1
## [46] splines_3.1.3 stringr_0.6.2 survival_2.38-1
## [49] tools_3.1.3 XML_3.98-1.1 xtable_1.7-4
## [52] XVector_0.6.0