Friday, August 19, 2016

GOSeq code for GO term enrichment of differentially expressed genes

If you have an output of the DESeq analysis of differentially expressed genes:

    X                id   baseMean  baseMeanA  baseMeanB foldChange
1   1 ENSG00000223972.4  0.3089417  0.5406481  0.0000000   0.000000
2   2 ENSG00000227232.4 51.6856965 46.1776533 59.0297541   1.278319
3   3 ENSG00000243485.2  0.0000000  0.0000000  0.0000000         NA
4   4 ENSG00000237613.2  0.0000000  0.0000000  0.0000000         NA
5   5 ENSG00000268020.2  0.0000000  0.0000000  0.0000000         NA
6   6 ENSG00000240361.1  0.0000000  0.0000000  0.0000000         NA
7   7 ENSG00000186092.4  0.0000000  0.0000000  0.0000000         NA
8   8 ENSG00000238009.2  0.5177663  0.2550477  0.8680578   3.403512
9   9 ENSG00000239945.1  0.0000000  0.0000000  0.0000000         NA
10 10 ENSG00000233750.3  0.4583062  0.8020358  0.0000000   0.000000
   log2FoldChange      pval      padj
1            -Inf 0.7521333 1.0000000
2       0.3542475 0.1023805 0.6247104
3              NA        NA        NA
4              NA        NA        NA
5              NA        NA        NA
6              NA        NA        NA
7              NA        NA        NA
8       1.7670242 0.6498254 1.0000000
9              NA        NA        NA
10           -Inf 0.5045118 1.0000000

and want to perform GOSeq analysis, that corrects for the length bias (shown in figure below - gene size effect shown in 400 gene bins as proportion of differentially expressed genes),



use the following code that will correct the ENSEMBL gene IDs (remove appended version number) and select genes with significant adjusted p-value:

library(goseq)
table<-read.delim("Result_Table.csv", header=T, sep=",")
table$id<-gsub("\\.[0-9].*", "",table$id)
genes=as.integer(table$padj<.05)
genes[is.na(genes)] <- 0
names(genes)=table$id
head(genes)
head(supportedOrganisms())
supportedOrganisms()[supportedOrganisms()$Genome=="hg19",]
pwf=nullp(genes,"hg19","ensGene")
GO.wall=goseq(pwf,"hg19","ensGene")
head(GO.wall)


Output in GO.wall table are significant GO terms:

> head(GO.wall)
        category over_represented_pvalue under_represented_pvalue numDEInCat
11673 GO:0044421            1.202999e-27                        1        610
2646  GO:0005576            1.593504e-26                        1        683
2677  GO:0005615            1.531784e-25                        1        261
11741 GO:0044699            1.158197e-22                        1       1764
11767 GO:0044763            3.101037e-22                        1       1645
5893  GO:0016477            1.439573e-21                        1        267
      numInCat                             term ontology
11673     3665        extracellular region part       CC
2646      4353             extracellular region       CC
2677      1302              extracellular space       CC
11741    13254          single-organism process       BP
11767    12106 single-organism cellular process       BP
5893      1176                   cell migration       BP

No comments:

Post a Comment