Data analysis step 8: Pathway analysis with GSEA
We use Gene Set Enrichment Analysis (GSEA) because it can detect pathway changes more sensitively and robustly than some methods. A 2013 paper compared a bunch of gene set analyses software with microarrays and is worth a look.
Generate a rank file
The rank file is a list of detected genes and a rank metric score. At the top of the list are genes with the "strongest" up-regulation, at the bottom of the list are the genes with the "strongest" down-regulation and the genes not changing are in the middle. The metric score I like to use is the sign of the fold change multiplied by the inverse of the p-value, although there may be better methods out there (link).Here is what the first 10 lines in the DESeq output look like:
 | |
|
#!/bin/bash
DGE=$1
RNK=`echo $DGE | sed 's/.xls/.rnk/'`
sed 1d $DGE \
| sort -k7g \
| cut -d '_' -f2- \
| awk '!arr[$1]++' \
| awk '{OFS="\t"}
{ if ($6>0) printf "%s\t%4.3e\n", $1, 1/$7 ;
else printf "%s\t%4.3e\n", $1, -1/$7 }' \
| sort -k2gr > $RNK
And here's the top and bottom 10 lines of the rank file:
$cat DEB_DESeq.rnk | (head;tail)
AC011899.9 3.352e+64
CAMK1D 1.471e+51
HSD11B1 1.919e+48
GNG4 2.279e+43
CPLX1 3.873e+40
CGNL1 1.200e+39
APCDD1 1.237e+35
MANEAL 9.303e+32
MKX 4.441e+28
HOXA13 1.007e+26
HLA-B -2.908e+132
PRAME -3.063e+165
LYZ -9.701e+185
CPXM1 -1.169e+186
HSH2D -1.658e+186
TGFBI -8.856e+210
APOC2 -9.633e+243
APOC4-APOC2 -1.984e+245
COL1A2 -1.298e+274
SLC12A3 -5.422e+304
GNG4 2.279e+43
CPLX1 3.873e+40
CGNL1 1.200e+39
APCDD1 1.237e+35
MANEAL 9.303e+32
MKX 4.441e+28
HOXA13 1.007e+26
HLA-B -2.908e+132
PRAME -3.063e+165
LYZ -9.701e+185
CPXM1 -1.169e+186
HSH2D -1.658e+186
TGFBI -8.856e+210
APOC2 -9.633e+243
APOC4-APOC2 -1.984e+245
COL1A2 -1.298e+274
SLC12A3 -5.422e+304
Run GSEA using the Java GUI
Download the javaGSEA Desktop Application from the Broad website (you will need to register). Open a terminal and cd into the directory containing the jar file and write (substitute the XX for appropriate version)
java -jar gsea2-XX.jar
Hit the "Load data" tab and browse for the rnk file we just generated.
Next, select "GSEA Preranked" from the "Tools" pull-down menu.
Here is the list of parameters I used:
Gene sets database: Initially, I like to use REACTOME to see that GSEA is working as well as quickly give a broad idea of the biology going on. Later on we will run with everything in the MSigDB database.
Number of permutations: 1000 is normally OK.
Ranked List: Ensure that GSEA is looking at the newly generated rnk file.
Collapse dataset to gene symbols: False
Chip platform: We're using is the gene symbols, you can leave this blank and GSEA will download the gene symbol chip file from the Broad website.
Enrichment statistic: Classic, this behaves better when we have extreme p-values such as this case
Max gene set size: 5000
Min gene set size: 15
Collapsing mode for probe sets>1 gene: Max_probe
Normalisation mode: meandiv
Omit features with no symbol: True
REACTOME result
From the initial REACTOME data base of 674 sets, 242 were removed because there either were less than 15 members detected or more than 5000. This left us with 432 sets to analyse. There were 220 sets with FDR p-value (called q-value here) ≤0.05. From these, 125 were up-regulated and 95 were down-regulated.
This result shows up-regulation of cell cycle and DNA repair pathways and down-regulation of innate immune response.
 | |
|
 | |
|
MSigDB Result
In the case of MSigDB, we began with 10295 gene sets and 2148 were filtered out (too few or too many), leaving 8147 sets. From these, 1613 were up-regulated (FDR p-value≤0.05) and 3692 were down-regulated.
 | |
|
 | |
|
