== Introduction ==

Given a ChIP-seq experiment for a transcription factor (TF) of interest and a microarray or RNA-Seq experiment comparing WT samples to KO or RNAi for the same TF, the goal of this analysis is to find out if the genes potentially regulated by the TF are changing in expression and if there is an enrichment of genes differentially expressed (DE) within the genes bound by the TF and vice versa.

We will be using the term “genes bound” or “genes potentially regulated by the TF” to refer to genes that are within a certain distance from a ChIP-Seq peak. This is defined by the user. Some examples of criteria to define “genes bound” by the TF would be: 1) The region from 2kb upstream to 500 pb downstream of the TSS overlaps with a peak. 2) The gene body or the 10Kb upstream or 10 Kb downstream of the gene overlap with a peak.


== Step 1: Define peaks ("bound regions") from a ChIP-Seq experiment == 

[wiki:SOPs/chip_seq_peaks Follow the standard operating procedure for identifying ChIP-Seq peaks]

== Step 2: Find genes differentially expressed from an expression microarray experiment == 

Follow standard operating procedures for [wiki:SOPs/normalize_ma normalizing microarray experiments] and [wiki:SOPs/diff_microarry identifying differentially expressed genes] to get a list of genes differentially expressed or a list of all the genes and their fold change on the microarray experiment.

For RNA-seq data, you can follow the steps in [wiki:SOPs/rna-seq-diff-expressions] to get a list of genes differentially expressed or a list of all the genes and their fold change.

== Step 3: Find out if the genes potentially regulated by the TF are changing in expression ==

'''Method 1'''. Use R to draw cumulative distribution function (cdf) plots that show the distribution of expression changes for all genes and for genes bound. If binding results in activation or repression of the genes bound we should see a shift in the distributions. We can test if the difference is significant with a Kolmogorov–Smirnov (KS) test.

Naming: FC stands for fold change in expression between control cells and cells where the TF of interest has been knock down. GenesBound is the list of genes potentially regulated by the TF according to the ChIP-seq experiment (see introduction above). AllTranscripts.FC.txt and GenesBound.FC.txt files have the FC for all genes or for the subset of genes bound. We plot both distributions and look for differences.

The input files used for the R code below are like:

AllTranscripts.FC.txt

{{{
 NM_000022        0.13    -0.02   -0.53     0.29   -0.06   -0.22   -0.18   -0.07
 NM_000073        0.03    -0.39   -0.42     0.29    0.01   -0.23   -0.02   -0.14
 NM_000188       -0.18    -0.43   -0.41    -0.72   -0.04   -0.29   -0.38   -0.47
 NM_000270        0.31     0.54    0.50     0.67    0.32    0.18    0.31    0.54
 NM_000314       -0.12     0.13    0.04    -0.23    0.03    0.14    0.11   -0.06
 NM_000401        0.00    -0.40    0.13    -0.06   -0.32    0.10    0.15    0.18
 NM_000404       -0.08     0.17    0.30    -0.20   -0.03    0.22   -0.08    0.11
    .....
}}}

GenesBound.FC.txt
{{{
 NM_000073       0.03    -0.39  -0.42   0.29    0.01  -0.23   -0.02   -0.14
 NM_000401       0.00    -0.40   0.13  -0.06   -0.32   0.10    0.15    0.18
 ...
}}}

These input files have the logFC of 8 know-down experiments. The expression changes we are interested in plotting are on the 6th expression column.
{{{
 control <- read.delim("AllTranscripts.FC.txt")
 GenesBound <- read.delim("GenesBound.FC.txt")
 png(file="cdfPlot.png")
 plot(ecdf(control[,6]), xlab="logFC_hsTF", ylab="Cumulative Fraction")
 lines(ecdf(GenesBound[,6]), verticals=T,do.points=FALSE, col="purple", lty=1, xlim=c(-6,6))
 dev.off()
}}}
{{{
 # Perform a KS test
 ks.test(control[,6], GenesBound[,6])
}}}

'''Method 2'''. Pick a cut-off for genes that are differentially expressed and test if the bound genes are enriched for differentially expressed genes using Fisher's exact test.

For example: Of the 20,000 transcripts assayed by the microarray, 700 were differentially expressed and 400 of the DE expressed genes were also found to be bound by the TF. The total number of genes bound by the TF was 3000.

{{{
              DE  non-DE  Totals
   Bound      400   2600  3000
   Not-Bound  300  16700  17000
   Totals     700  19300  20000
}}}

You can perform a Chi-squared test or a Fisher's exact test in R:
{{{
   data <- matrix(c(400, 2600, 300, 16700), ncol=2, byrow=T)
   chisq.test(data)
   fisher.test(data)
}}}

'''Method 3'''. Use the GSEA tool from the Broad Institute. Use the raw array data and genes bound as the gene data set to look for enrichment [http://www.broadinstitute.org/gsea/index.jsp GSEA]

'''Method 4'''. For visualization, draw a scatterplot on log2 gene expression levels with WT and KO on X and Y axis, highlight the bound genes on top of all genes. This can be achieved with simple R codes. In the example below, the two text files each has three columns, with geneID, WT and KO as header.  

{{{
 control=read.delim("all_genes_exp.txt")
 GenesBound = read.delim("GenesBound_exp.txt")
 par(pty="s")
 plot(log2(control$WT), log2(control$KO))
 abline(0,1)
 points (log2(GenesBound$WT, log2(GenesBound$KO),  col="red")
}}}