annotate goprofiles.R @ 0:6eeb2fb0c4bd draft default tip

planemo upload

1 # Load necessary libraries
2 library("org.Hs.eg.db", quietly=TRUE)
3 library("goProfiles", quietly=TRUE)
4 
5 # Read file and return file content as data.frame?
6 readfile = function(filename, header) {
7   if (header == "true") {
8     # Read only the first two lines of the files as data (without headers):
9     headers <- read.table(filename, nrows = 1, header = FALSE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE)
10     #print("header")
11     #print(headers)
12     # Create the headers names with the two (or more) first rows, sappy allows to make operations over the columns (in this case paste) - read more about sapply here :
13     #headers_names <- sapply(headers, paste, collapse = "_")
14     #print(headers_names)
15     #Read the data of the files (skipping the first 2 rows):
16     file <- read.table(filename, skip = 1, header = FALSE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE)
17     #print(file[1,])
18     #And assign the headers of step two to the data:
19     names(file) <- headers
20   }
21   else {
22     file <- read.table(filename, header = FALSE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE)
23   }
24   return(file)
25 }
26 
27 #filename = "/Users/LinCun/Documents/ProteoRE/usecase1/Check/HPA.Selection.134.txt"
28 #test = readfile(filename)
29 #str(test)
30 #str(test$Gene.names)
31 getprofile = function(ids, id_type, level) {
32   
33   # Check if level is number
34   if (! as.numeric(level) %% 1 == 0) {
35     stop("Please enter an integer for level")
36   }
37   else {
38     level = as.numeric(level)
39   }
40   #genes = as.vector(file[,ncol])
41   
42   # Extract Gene Entrez ID
43   if (id_type == "Entrez") {
44     id = select(org.Hs.eg.db, ids, "ENTREZID", multiVals = "first")
45     genes_ids = id$ENTREZID[which( ! is.na(id$ENTREZID))]
46   }
47   else {
48     genes_ids = c()
49     id = select(org.Hs.eg.db, ids, "ENTREZID", "UNIPROT", multiVals = "first")
50     #print(id[[1]][1])
51     genes_ids = id$ENTREZID[which( ! is.na(id$ENTREZID))]
52     # IDs that have NA ENTREZID
53     NAs = id$UNIPROT[which(is.na(id$ENTREZID))]
54     print("IDs unable to convert to ENTREZID: ")
55     print(NAs)
56   }
57   #print(genes_ids)
58   # Convert Protein IDs into entrez ids
59   
60   # for (i in 1:length(id$UNIPROT)) {
61   #   print(i)
62   #   if (is.na(id[[2]][i])) {
63   #     print(id[[2]][i])
64   #   }
65   # }
66   # a = id[which(id$ENTREZID == "NA"),]
67   # print(a)
68   # print(a$UNIPROT)
69   #print(id[[1]][which(is.na(id$ENTREZID))])
70   #print(genes_ids)
71   # for (gene in genes) {
72   #   #id = as.character(mget(gene, org.Hs.egALIAS2EG, ifnotfound = NA))
73   #   id = select(org.Hs.eg.db, genes, "ENTREZID", "UNIPROT")
74   #   print(id)
75   #   genes_ids = append(genes_ids, id$ENTREZID)
76   # }
77   #print(genes_ids)
78   
79   # Create basic profiles
80   profile.CC = basicProfile(genes_ids, onto='CC', level=level, orgPackage="org.Hs.eg.db", empty.cats=F, ord=T, na.rm=T)
81   profile.BP = basicProfile(genes_ids, onto='BP', level=level, orgPackage="org.Hs.eg.db", empty.cats=F, ord=T, na.rm=T)
82   profile.MF = basicProfile(genes_ids, onto='MF', level=level, orgPackage="org.Hs.eg.db", empty.cats=F, ord=T, na.rm=T)
83   profile.ALL = basicProfile(genes_ids, onto='ANY', level=level, orgPackage="org.Hs.eg.db", empty.cats=F, ord=T, na.rm=T)
84   
85   # Print profile
86   # printProfiles(profile)
87   
88   return(c(profile.CC, profile.MF, profile.BP, profile.ALL))
89 }
90 
91 # Plot profiles to PNG
92 plotPNG = function(profile.CC = NULL, profile.BP = NULL, profile.MF = NULL, profile.ALL = NULL, per = TRUE, title = TRUE) {
93   if (!is.null(profile.CC)) {
94     png("profile.CC.png")
95     plotProfiles(profile.CC, percentage=per, multiplePlots=FALSE, aTitle=title)
96     dev.off()
97   }
98   if (!is.null(profile.BP)) {
99     png("profile.BP.png")
00     plotProfiles(profile.BP, percentage=per, multiplePlots=FALSE, aTitle=title)
01     dev.off()
02   }
03   if (!is.null(profile.MF)) {
04     png("profile.MF.png")
05     plotProfiles(profile.MF, percentage=per, multiplePlots=FALSE, aTitle=title)
06     dev.off()
07   }
08   if (!is.null(profile.ALL)) {
09     png("profile.ALL.png")
10     plotProfiles(profile.ALL, percentage=per, multiplePlots=T, aTitle=title)
11     dev.off()
12   }
13 }
14 
15 # Plot profiles to JPEG
16 plotJPEG = function(profile.CC = NULL, profile.BP = NULL, profile.MF = NULL, profile.ALL = NULL, per = TRUE, title = TRUE) {
17   if (!is.null(profile.CC)) {
18     jpeg("profile.CC.jpeg")
19     plotProfiles(profile.CC, percentage=per, multiplePlots=FALSE, aTitle=title)
20     dev.off()
21   }
22   if (!is.null(profile.BP)) {
23     jpeg("profile.BP.jpeg")
24     plotProfiles(profile.BP, percentage=per, multiplePlots=FALSE, aTitle=title)
25     dev.off()
26   }
27   if (!is.null(profile.MF)) {
28     jpeg("profile.MF.jpeg")
29     plotProfiles(profile.MF, percentage=per, multiplePlots=FALSE, aTitle=title)
30     dev.off()
31   }
32   if (!is.null(profile.ALL)) {
33     jpeg("profile.ALL.jpeg")
34     plotProfiles(profile.ALL, percentage=per, multiplePlots=FALSE, aTitle=title)
35     dev.off()
36   }
37 }
38 
39 # Plot profiles to PDF
40 plotPDF = function(profile.CC = NULL, profile.BP = NULL, profile.MF = NULL, profile.ALL = NULL, per = TRUE, title = TRUE) {
41   if (!is.null(profile.CC)) {
42     pdf("profile.CC.pdf")
43     plotProfiles(profile.CC, percentage=per, multiplePlots=FALSE, aTitle=title)
44     dev.off()
45   }
46   if (!is.null(profile.BP)) {
47     pdf("profile.BP.pdf")
48     plotProfiles(profile.BP, percentage=per, multiplePlots=FALSE, aTitle=title)
49     dev.off()
50   }
51   if (!is.null(profile.MF)) {
52     pdf("profile.MF.pdf")
53     plotProfiles(profile.MF, percentage=per, multiplePlots=FALSE, aTitle=title)
54     dev.off()
55   }
56   if (!is.null(profile.ALL)) {
57     #print("all")
58     pdf("profile.ALL.pdf")
59     plotProfiles(profile.ALL, percentage=per, multiplePlots=FALSE, aTitle=title)
60     dev.off()
61   }
62 }
63 
64 goprofiles = function() {
65   args = commandArgs(trailingOnly = TRUE)
66   #print(args)
67   # arguments: filename.R inputfile ncol "CC,MF,BP,ALL" "PNG,JPEG,PDF" level "TRUE"(percentage) "Title"
68   if (length(args) != 8) {
69     stop("Not enough/Too many arguments", call. = FALSE)
70   }
71   else {
72     input_type = args[2]
73     if (input_type == "text") {
74       input = strsplit(args[1], "\\s+")[[1]]
75     }
76     else if (input_type == "file") {
77       filename = strsplit(args[1], ",")[[1]][1]
78       ncol = strsplit(args[1], ",")[[1]][2]
79       # Check ncol
80       if (! as.numeric(gsub("c", "", ncol)) %% 1 == 0) {
81         stop("Please enter an integer for level")
82       }
83       else {
84         ncol = as.numeric(gsub("c", "", ncol))
85       }
86       header = strsplit(args[1], ",")[[1]][3]
87       # Get file content
88       file = readfile(filename, header)
89       # Extract Protein IDs list
90       input = c()
91       for (row in as.character(file[,ncol])) {
92         input = c(input, strsplit(row, ";")[[1]][1])
93       }
94     }
95     id_type = args[3]
96     ontoopt = strsplit(args[4], ",")[[1]]
97     #print(ontoopt)
98     #plotopt = strsplit(args[3], ",")
99     plotopt = args[5]
00     level = args[6]
01     per = as.logical(args[7])
02     title = args[8]
03 
04     profiles = getprofile(input, id_type, level)
05     profile.CC = profiles[1]
06     #print(profile.CC)
07     profile.MF = profiles[2]
08     #print(profile.MF)
09     profile.BP = profiles[3]
10     #print(profile.BP)
11     profile.ALL = profiles[-3:-1]
12     #print(profile.ALL)
13     #c(profile.ALL, profile.CC, profile.MF, profile.BP)
14     if ("CC" %in% ontoopt) {
15       if (grepl("PNG", plotopt)) {
16         plotPNG(profile.CC=profile.CC, per=per, title=title)
17       }
18       if (grepl("JPEG", plotopt)) {
19         plotJPEG(profile.CC = profile.CC, per=per, title=title)
20       }
21       if (grepl("PDF", plotopt)) {
22         plotPDF(profile.CC = profile.CC, per=per, title=title)
23       }
24     }
25     if ("MF" %in% ontoopt) {
26       if (grepl("PNG", plotopt)) {
27         plotPNG(profile.MF = profile.MF, per=per, title=title)
28       }
29       if (grepl("JPEG", plotopt)) {
30         plotJPEG(profile.MF = profile.MF, per=per, title=title)
31       }
32       if (grepl("PDF", plotopt)) {
33         plotPDF(profile.MF = profile.MF, per=per, title=title)
34       }
35     }
36     if ("BP" %in% ontoopt) {
37       if (grepl("PNG", plotopt)) {
38         plotPNG(profile.BP = profile.BP, per=per, title=title)
39       }
40       if (grepl("JPEG", plotopt)) {
41         plotJPEG(profile.BP = profile.BP, per=per, title=title)
42       }
43       if (grepl("PDF", plotopt)) {
44         plotPDF(profile.BP = profile.BP, per=per, title=title)
45       }
46     }
47     
48     #if (grepl("PNG", plotopt)) {
49     # plotPNG(profile.ALL = profile.ALL, per=per, title=title)
50     #}
51     #if (grepl("JPEG", plotopt)) {
52     # plotJPEG(profile.ALL = profile.ALL, per=per, title=title)
53     #}
54     #if (grepl("PDF", plotopt)) {
55     # plotPDF(profile.ALL = profile.ALL, per=per, title=title)
56     #}
57   }
58   
59 }
60 
61 goprofiles()
62 
63 #Rscript go.R ../proteinGroups_Maud.txt "1" "CC" "PDF" 2 "TRUE" "Title"