Replace GOenrichment_proj_test.R

scripts/GOenrichment_proj.R

@@ -49,12 +49,12 @@ library(multtest, verbose = FALSE)
 
 #list with samples names
 multiple_dexp<- read.delim(args[1],stringsAsFactors = F) 
-#multiple_dexp<- read.delim("multiple_input.txt",stringsAsFactors = F) 
+#multiple_dexp<- read.delim("Padj_21genes.txt",stringsAsFactors = F) 
 sample_names <-colnames(multiple_dexp)# all names
-sample_names[2:length(sample_names)]
 pcutoff <- 0.05
 #gene names from DGE analysis output
 genes <- multiple_dexp$Gene
+genes <- genes[2:length(genes)]
 
 
 #useful info: #multiple_dexp[,c(1,which(colnames(multiple_dexp) == current_leaf))]
@@ -64,17 +64,19 @@ for(leaf in sample_names[2:length(sample_names)]){
   
   #filtering genes.replace to 1 if above threshold, otherwise 0 (significant)
   filter_multiple <- ifelse(multiple_dexp[,which(colnames(multiple_dexp) == current_leaf)] < pcutoff,1,0) # adj.P.Val name of the column 
-  names(filter_multiple) <- genes
-  geneList<- filter_multiple
+  geneList <- filter_multiple
+  names(geneList) <- unlist(multiple_dexp$Gene)
+  #print(geneList)
   
   # create topGOdata object
   GOdata_sample <- new("topGOdata",
                 ontology = "BP",
+                nodeSize = 5,
                 allGenes = geneList,
                 geneSelectionFun = function(x)(x == 1),
                 annot = annFUN.org, mapping = "org.At.tair.db")
   
-  #Run the Fishers exact Test (raw p-vals)
+  #Run the Fisher´s exact Test (raw p-vals)
   fisher_sample <- runTest(GOdata_sample, algorithm = "elim", statistic = "fisher")
   
   sample_result <- GenTable(GOdata_sample, raw.p.value = fisher_sample, 
@@ -86,7 +88,7 @@ for(leaf in sample_names[2:length(sample_names)]){
   #http://www.bioconductor.org/packages/release/bioc/html/multtest.html
   raw_pvalues <- as.double(sample_result$raw.p.value)
   cor_pval <- mt.rawp2adjp(raw_pvalues, "BY")
-  adjp <-cor_pval$adjp #if want order_adjp <- res_FDR$adjp[order(res_FDR$index), ]
+  adjp <-cor_pval$adjp[order(cor_pval$index), ]
   
   #Add FDR p value to topGO output
   FDR_by_sample <- sample_result
@@ -100,7 +102,8 @@ for(leaf in sample_names[2:length(sample_names)]){
   
   #write table with information per leaf
   table = paste(current_leaf,"txt",sep=".")
-  write.table(FDR_by_sample, col.names=T, row.names=F, file = paste("/enrichment/",table,sep=""), sep ="\t",
+  write.table(FDR_by_sample, col.names=T, row.names=F, file = paste("enrichment/",table,sep=""), sep ="\t",
               quote = F)
   
 }
+