path_vec<-function(nmerge_tab,condition,group1,group2,kg.sets){
index1<-is.na(factor(condition,levels=group1))==FALSE
edata1=nmerge_tab[,index1]
edatamean1=apply(edata1,1,function(x) mean(x) )
ex1=data.frame(mrna=rownames(nmerge_tab),edatamean1)
diff1=edata1
index2<-is.na(factor(condition,levels=group2))==FALSE
edata2=nmerge_tab[,index2]
edatamean2=apply(edata2,1,function(x) mean(x) )
ex2=data.frame(mrna=rownames(nmerge_tab),edatamean2)
diff2=edata2
#####把大鼠所有通路信息保存到kg.sets???
# path.list<-keggLink("pathway",'rno')
# genes=gsub(paste('rno',":",sep=""),"",names(path.list))
# paths=gsub(paste('path',":",sep=""),"",path.list)
# kg.sets=split(genes,paths)
####### 通路富集函数
enrichkegg<-function(ensemblid){
geneid<-ensemblid
keggpath<-enrichKEGG(as.character(geneid),organism = "rno", pvalueCutoff = 0.05,pAdjustMethod = "none", qvalueCutoff = Inf)
pathway<-summary(keggpath)
return(pathway)
}
pathway1=enrichkegg(rownames(diff1))[,1]
pathway2=enrichkegg(rownames(diff2))[,1]
####### ENSEMBL ID to ENTREZ ID
pathway<-intersect(pathway1,pathway2)
# transID<-function(id){
# geneid<-bitr(id,fromType = 'ENSEMBL',toType = 'ENTREZID',OrgDb = 'org.Rn.eg.db', drop = FALSE)
# return(geneid)
# }
###### 建立通路定量表达谱
EstablishProfile<-function(id,expression,pathway){
geneid<-id
# na<-geneid[is.na(geneid[,2])==TRUE,][,1]
profile<-data.frame()
for(i in 1:length(id)){
exvector<-c()
proid<-geneid
# if(id[i]%in%na){proid<-'NA'}
# else{proid<-bitr(id[i],fromType = "ENSEMBL",toType = "ENTREZID",OrgDb = 'org.Rn.eg.db')[1,2]}
for(j in 1:length(pathway)){
if(proid%in%kg.sets[[pathway[j]]]){exvector<-c(exvector,expression[expression[,1]==id[i],2])}
else{exvector<-c(exvector,0)}
}
if(sum(exvector)==0){next}
profile<-rbind(profile,exvector)
}
colnames(profile)<-pathway
return(profile)
}
profile1=EstablishProfile(rownames(diff1),ex1,pathway)
profile2=EstablishProfile(rownames(diff2),ex2,pathway)
list=list(profile1,profile2)
return(list)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.