Entering edit mode
10.4 years ago
Ming Tommy Tang
★
4.5k
Hi biostarers,
I was following the GAGE tutorial for pathway enrichment analysis after DESeq. 1724 differentially expressed genes were used for pathway analysis.
res <- nbinomTest( cds, 'control, 'treat' )
resSig <- res[ res$padj < 0.01 & (res$log2FoldChange >1| res$log2FoldChange < -1), ]
resSig <- na.omit(resSig)
require(gage)
data(kegg.gs)
deseq.fc<- resSig$log2FoldChange
names(deseq.fc)<- resSig$id
sum(is.infinite(deseq.fc)) # there are some infinite numbers, if use DESeq2, no such problem.
deseq.fc[deseq.fc>10]=10
deseq.fc[deseq.fc<-10]=-10
exp.fc<- deseq.fc
#kegg.gsets works with 3000 KEGG speicies
data(korg)
head(korg[,1:3], n=20)
#let's get the annotation files for mouse and convert the gene set to gene symbol format
kg.mouse<- kegg.gsets("mouse")
kegg.gs<- kg.mouse$kg.sets[kg.mouse$sigmet.idx]
lapply(kegg.gs[1:3],head)
# egSymb is only for human data, so eg2sym and sym2eg functions are only for human data.
#data(egSymb)
#kegg.gs.sym<- lapply(kegg.gs, eg2sym)
#lapply(kegg.gs.sym[1:3],head)
# to convert IDs among gene/transcript ID to Entrez GeneID or reverse, use eg2id and id2eg in the pathview package #written by the same person.
library(pathview)
data(bods)
bods
gene.symbol.eg<- id2eg(ids=names(exp.fc), category='SYMBOL', org='Mm') # convert the gene symbol to Entrez Gene ID
head(gene.symbol.eg, n=100)
head(gene.symbol.eg[,2], n=10)
names(exp.fc)<- gene.symbol.eg[,2]
fc.kegg.p<- gage(exp.fc, gsets= kegg.gs, ref=NULL, samp=NULL)
sel<- fc.kegg.p$greater[,"q.val"] < 0.1 & !is.na(fc.kegg.p$greater[,"q.val"])
table(sel)
sel.l<- fc.kegg.p$less[,"q.val"] < 0.1 & !is.na(fc.kegg.p$greater[,"q.val"])
table(sel.l)
> table(sel.l)
sel.l
FALSE
202
but no pathways are significantly changed table(sel)
...anything wrong with my code? Thank you!