R/path_vec.R

In github-gs/Quantitative-pathway-analysis: Quantitative Pathway Analysis

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)
}