Monday, April 4, 2016

Calculate Fisher test p-values for overlaps of two bed files in a given genomic background

This is a code that will generate a Fisher test p-value that shows significance for the overlap of two sets of genomic regions (for example from ChIP-Seq experiments). To calculate Fisher p-value you would need three BED files:
1. First bed file to overlap
2. Second bed file for overlap
3. Bed file that will serve as genomic background for overlaps.

As genomic background you can use e.g. combined ENCODE set of open chromatin regions or a similar data set. Genomic background is necessary to calculate the constituents of the Fisher test contingency matrix: overlap of A and B in the background BG set, A but not B in BG, B but not A in BG, BG without A and B.

Code:

bedtools intersect -a A.bed -b background.bed > 1_BG
bedtools intersect -a B.bed -b background.bed > 2_BG

bedtools intersect -wa -a 1_BG -b 2_BG | uniq > 1_2_BG
bedtools intersect -wa -a 2_BG -b 1_BG | uniq > 2_1_BG
bedtools intersect -v -a 1_BG -b 1_2_BG | uniq > 1_BG_no2
bedtools intersect -v -a 2_BG -b 2_1_BG | uniq > 2_BG_no1
bedtools intersect -v -a background.bed -b 1_2_BG | uniq > BG_no1_no2

wc -l 1_2_BG
wc -l 2_1_BG
wc -l 1_BG_no2
wc -l 2_BG_no1
wc -l BG_no1_no2
Next, create a contingency matrix in R using the outputs of wc -l command, e.g.

365 1_2_BG
360 2_1_BG
45854 1_BG_no2
3670 2_BG_no1
2233184 BG_no1_no2
Create R code:

AHR_TCF21  <-
matrix(c(365, 3670, 45854, 2233184),
       nrow = 2,
       dimnames = list(Guess = c("b+", "b-"),
                       Truth = c("a+", "a-")))

AHR_TCF21.fisher <- fisher.test(AHR_TCF21)
str(AHR_TCF21.fisher)
Output in R:

> AHR_TCF21  <-
+ matrix(c(365, 3670, 45854, 2233184),
+        nrow = 2,
+        dimnames = list(Guess = c("b+", "b-"),
+                        Tru .... [TRUNCATED]

> AHR_TCF21.fisher <- fisher.test(AHR_TCF21)

> str(AHR_TCF21.fisher)
List of 7
 $ p.value    : num 1.88e-121
 $ conf.int   : atomic [1:2] 4.34 5.4
  ..- attr(*, "conf.level")= num 0.95
 $ estimate   : Named num 4.84
  ..- attr(*, "names")= chr "odds ratio"
 $ null.value : Named num 1
  ..- attr(*, "names")= chr "odds ratio"
 $ alternative: chr "two.sided"
 $ method     : chr "Fisher's Exact Test for Count Data"
 $ data.name  : chr "AHR_TCF21"
 - attr(*, "class")= chr "htest"
In this example, we can see that the overlap is highly significant, with p-value of 1.88e-121 and CI 4.34-5.4.


No comments:

Post a Comment