Upload New File

scripts/KEGG1_enrichment.R

+.libPaths(c('/local/data/course/project/araformatics/pipeline/test_R',.libPaths()))
+#!/usr/bin/env Rscript
+#BIF30806 Project group: Araformatic
+#(2021)
+#KEGG pathway enrichment
+#Scrip to run enrichKEGG, from clusterProfile library to get enriched 
+#KEGG pathways and their adjusted p values.
+#Usage: Rscript KEGG3.R input_file.txt
+#    input_file: txt file, genes as rows, leafs as columns and values
+#    are adjusted p value from DGE
+#Output: txt file for each leaf(sample) in the input file
+
+#source("https://bioconductor.org/biocLite.R")
+#biocLite(c("KEGGREST", "org.At.tair.db", "Rgraphviz","clusterProfile"))
+#if (!requireNamespace("BiocManager", quietly = TRUE))
+  install.packages("BiocManager")
+#BiocManager::install("org.At.tair.db")
+#BiocManager::install("signatureSearch")
+#BiocManager::install("clusterProfiler")
+library("org.At.tair.db")
+library(KEGGREST)
+library(clusterProfiler)
+
+#Take arguments from command line
+args = commandArgs(trailingOnly = TRUE)
+
+
+#Get A.thaliana pathways from KEGG 
+pathway_list <- keggList("pathway", "ath")
+
+#genes per pathway
+pathway_codes <- sub("path:", "", names(pathway_list)) #remove path: string
+
+genes_pathway <- sapply(pathway_codes,
+                        function(pwid){
+                          path <- keggGet(pwid)
+                          if (is.null(path[[1]]$GENE)) return(NA)
+                          path2 <- path[[1]]$GENE[c(TRUE,FALSE)]
+                          path2 <- unlist(lapply(strsplit(path2,split=";",fixed =T),
+                                                 function(x)x[1])) #atomic vector, apply function to all [1] in path2
+                          return(path2)
+                        })
+all_geneIDs <- keys(org.At.tair.db)
+
+#Read DE file
+#MULTIPLE SAMPLES IN ONE FILE
+#read file with p adj value data
+de_multiple<- read.delim("multiple_input.txt",stringsAsFactors = F,row.names = 1) 
+de_multiple<- read.delim(args[1],stringsAsFactors = F,row.names = 1)
+sample_names <-colnames(de_multiple)
+
+for(leaf in sample_names){
+  current_leaf = leaf
+  print(current_leaf)
+  
+  #geneList
+  geneList <- rownames(de_multiple[which(de_multiple[which(colnames(de_multiple) == current_leaf)]<0.05),])
+  print(geneList)
+  #enrichment 
+  kk <- enrichKEGG(gene=geneList,organism="ath",keyType = "kegg",
+                   pAdjustMethod= "BY",universe=names(all_geneIDs),qvalueCutoff = 0.05)
+  
+  
+  #write table with enriched pathways
+  if(is.null(kk)){
+    next
+    }
+  else{
+    #enrichKEGG2 gives a A feaResult instance turn it into data frame
+    enrich_k = data.frame(kk@result)
+    
+    table = paste(current_leaf,"KEGG1",sep="_")
+    table = paste(table,"txt",sep=".")
+    write.table(enrich_k, col.names=T, row.names=F,
+                file = paste("/enrichment/",table,sep=""), sep ="\t",quote = F)
+    }
+
+
+  }
\ No newline at end of file