U343 T24 Radiation GOSeq analysis

PWF for DOWN regulated genes

go_title.down <- 'GO enrichment for U343 T24_vs_T21_DOWN'
pwf.down  = nullp(u343_t25_down.genes, "hg38", "ensGene", bias.data = lengthData.down)

PWF for UP regulated genes

pwf.up  = nullp(u343_t25_up.genes, "hg38", "ensGene", bias.data = lengthData.up)

GO.wall.up <- goseq(pwf.up,"hg38","ensGene",test.cats=c("GO:BP"))
GO.wall.up$padj <- p.adjust(GO.wall.up$over_represented_pvalue, method="BH")
GO.wall.up

GO.wall.up <- subset(GO.wall.up, GO.wall.up$padj<.05)
GO.wall.up[["ratio"]] <- GO.wall.up[["numDEInCat"]] / GO.wall.up[["numInCat"]]

GO.wall.down <- goseq(pwf.down, "hg38","ensGene",test.cats=c("GO:BP"))
GO.wall.down$padj <- p.adjust(GO.wall.down$over_represented_pvalue, method="BH")
GO.wall.down <- subset(GO.wall.down, GO.wall.down$padj<.05)
GO.wall.down[["ratio"]] <- GO.wall.down[["numDEInCat"]] / GO.wall.down[["numInCat"]]

Go enrichment DOWN regulated

No enricment for DOWN in KEGG or UP in KEGG/GO.

barplot(GO.wall.down, showCategory = 15)

dotplot(GO.wall.down)

KEGG.up <- goseq(pwf.up, 'hg38', 'ensGene', test.cats="KEGG")
KEGG.up$padj <- p.adjust(KEGG.up$over_represented_pvalue, method="BH")
KEGG.up.sig <- subset(KEGG.up, KEGG.up$padj<.05)
KEGG.up.sig[["ratio"]] <- KEGG.up.sig[["numDEInCat"]] / KEGG.up[["numInCat"]]

KEGG.down <- goseq(pwf.down, 'hg38', 'ensGene', test.cats="KEGG")
KEGG.down$padj <- p.adjust(KEGG.down$over_represented_pvalue, method="BH")
KEGG.down.sig <- subset(KEGG.down, KEGG.down$padj<.05)
KEGG.down.sig[["ratio"]] <- KEGG.down.sig[["numDEInCat"]] / KEGG.down.sig[["numInCat"]]

human <- useMart('ensembl', dataset = 'hsapiens_gene_ensembl')
attributes <- c('ensembl_gene_id',
                'entrezgene')

ensembl.entrez <- getBM(attributes,
                                  values="*", 
                                  mart = human,
                                  bmHeader=FALSE)
ensembl2entrez <- as.character(ensembl.entrez$entrezgene)
names(ensembl2entrez) <- ensembl.entrez$ensembl_gene_id
ensembl2entrez <- as.list(ensembl2entrez)

get_kegg <- function(l)
{
  return( keggLink("pathway", l))
}

entrez.kegg <- paste('hsa', unique(ensembl2entrez), sep=":")
entrez.kegg.splitted <- split( entrez.kegg, ceiling(seq_along( entrez.kegg)/20))
entrez.kegg.pathways <- lapply(entrez.kegg.splitted, get_kegg) 
entrez.kegg.pathways.red <- Reduce(c, entrez.kegg.pathways)
grepKEGG <- function(geneid, pathway_map){
  if (!is.na(geneid))
  {
    return (unique(unlist(pathway_map[names(pathway_map) == paste('hsa', geneid, sep=':')])))
  }
}
gene2cat <- lapply(ensembl2entrez,grepKEGG, entrez.kegg.pathways.red)

getPathwayName <- function(x){
  y <- keggGet(x)
  return (as.character(unlist(strsplit(y[[1]]$NAME, '-'))[1]))
}

pwf.down  <- nullp(u343_t25_down.genes, "hg38", "ensGene", bias.data = lengthData.down)
KEGG.down.custom <- goseq(pwf.down, gene2cat = gene2cat)
KEGG.down.custom$term <- unlist(lapply(KEGG.down.custom$category, getPathwayName))
KEGG.down.custom[["ratio"]] <- KEGG.down.custom[["numDEInCat"]] / KEGG.down.custom[["numInCat"]]

KEGG.down.custom$padj <- p.adjust(KEGG.down.custom$over_represented_pvalue, method="BH")
KEGG.down.custom.sig <- subset(KEGG.down.custom, KEGG.down.custom$padj<.05)

dotplot(KEGG.down.custom)

KEGG.up.custom <- goseq(pwf.up, gene2cat=gene2cat)
KEGG.up.custom$term <- unlist(lapply(KEGG.up.custom$category, getPathwayName))
KEGG.up.custom[["ratio"]] <- KEGG.up.custom[["numDEInCat"]] / KEGG.up.custom[["numInCat"]]

KEGG.up.custom$padj <- p.adjust(KEGG.up.custom$over_represented_pvalue, method="BH")
KEGG.up.custom.sig <- subset(KEGG.up.custom, KEGG.up.custom$padj<.05)

dotplot(KEGG.up.custom)