# HG changeset patch # User lnguyen # Date 1505567827 14400 # Node ID 6eeb2fb0c4bd14a398af4a639e839229211371ca planemo upload diff -r 000000000000 -r 6eeb2fb0c4bd goprofiles.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/goprofiles.R Sat Sep 16 09:17:07 2017 -0400 @@ -0,0 +1,263 @@ +# Load necessary libraries +library("org.Hs.eg.db", quietly=TRUE) +library("goProfiles", quietly=TRUE) + +# Read file and return file content as data.frame? +readfile = function(filename, header) { + if (header == "true") { + # Read only the first two lines of the files as data (without headers): + headers <- read.table(filename, nrows = 1, header = FALSE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE) + #print("header") + #print(headers) + # 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 : + #headers_names <- sapply(headers, paste, collapse = "_") + #print(headers_names) + #Read the data of the files (skipping the first 2 rows): + file <- read.table(filename, skip = 1, header = FALSE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE) + #print(file[1,]) + #And assign the headers of step two to the data: + names(file) <- headers + } + else { + file <- read.table(filename, header = FALSE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE) + } + return(file) +} + +#filename = "/Users/LinCun/Documents/ProteoRE/usecase1/Check/HPA.Selection.134.txt" +#test = readfile(filename) +#str(test) +#str(test$Gene.names) +getprofile = function(ids, id_type, level) { + + # Check if level is number + if (! as.numeric(level) %% 1 == 0) { + stop("Please enter an integer for level") + } + else { + level = as.numeric(level) + } + #genes = as.vector(file[,ncol]) + + # Extract Gene Entrez ID + if (id_type == "Entrez") { + id = select(org.Hs.eg.db, ids, "ENTREZID", multiVals = "first") + genes_ids = id$ENTREZID[which( ! is.na(id$ENTREZID))] + } + else { + genes_ids = c() + id = select(org.Hs.eg.db, ids, "ENTREZID", "UNIPROT", multiVals = "first") + #print(id[[1]][1]) + genes_ids = id$ENTREZID[which( ! is.na(id$ENTREZID))] + # IDs that have NA ENTREZID + NAs = id$UNIPROT[which(is.na(id$ENTREZID))] + print("IDs unable to convert to ENTREZID: ") + print(NAs) + } + #print(genes_ids) + # Convert Protein IDs into entrez ids + + # for (i in 1:length(id$UNIPROT)) { + # print(i) + # if (is.na(id[[2]][i])) { + # print(id[[2]][i]) + # } + # } + # a = id[which(id$ENTREZID == "NA"),] + # print(a) + # print(a$UNIPROT) + #print(id[[1]][which(is.na(id$ENTREZID))]) + #print(genes_ids) + # for (gene in genes) { + # #id = as.character(mget(gene, org.Hs.egALIAS2EG, ifnotfound = NA)) + # id = select(org.Hs.eg.db, genes, "ENTREZID", "UNIPROT") + # print(id) + # genes_ids = append(genes_ids, id$ENTREZID) + # } + #print(genes_ids) + + # Create basic profiles + profile.CC = basicProfile(genes_ids, onto='CC', level=level, orgPackage="org.Hs.eg.db", empty.cats=F, ord=T, na.rm=T) + profile.BP = basicProfile(genes_ids, onto='BP', level=level, orgPackage="org.Hs.eg.db", empty.cats=F, ord=T, na.rm=T) + profile.MF = basicProfile(genes_ids, onto='MF', level=level, orgPackage="org.Hs.eg.db", empty.cats=F, ord=T, na.rm=T) + profile.ALL = basicProfile(genes_ids, onto='ANY', level=level, orgPackage="org.Hs.eg.db", empty.cats=F, ord=T, na.rm=T) + + # Print profile + # printProfiles(profile) + + return(c(profile.CC, profile.MF, profile.BP, profile.ALL)) +} + +# Plot profiles to PNG +plotPNG = function(profile.CC = NULL, profile.BP = NULL, profile.MF = NULL, profile.ALL = NULL, per = TRUE, title = TRUE) { + if (!is.null(profile.CC)) { + png("profile.CC.png") + plotProfiles(profile.CC, percentage=per, multiplePlots=FALSE, aTitle=title) + dev.off() + } + if (!is.null(profile.BP)) { + png("profile.BP.png") + plotProfiles(profile.BP, percentage=per, multiplePlots=FALSE, aTitle=title) + dev.off() + } + if (!is.null(profile.MF)) { + png("profile.MF.png") + plotProfiles(profile.MF, percentage=per, multiplePlots=FALSE, aTitle=title) + dev.off() + } + if (!is.null(profile.ALL)) { + png("profile.ALL.png") + plotProfiles(profile.ALL, percentage=per, multiplePlots=T, aTitle=title) + dev.off() + } +} + +# Plot profiles to JPEG +plotJPEG = function(profile.CC = NULL, profile.BP = NULL, profile.MF = NULL, profile.ALL = NULL, per = TRUE, title = TRUE) { + if (!is.null(profile.CC)) { + jpeg("profile.CC.jpeg") + plotProfiles(profile.CC, percentage=per, multiplePlots=FALSE, aTitle=title) + dev.off() + } + if (!is.null(profile.BP)) { + jpeg("profile.BP.jpeg") + plotProfiles(profile.BP, percentage=per, multiplePlots=FALSE, aTitle=title) + dev.off() + } + if (!is.null(profile.MF)) { + jpeg("profile.MF.jpeg") + plotProfiles(profile.MF, percentage=per, multiplePlots=FALSE, aTitle=title) + dev.off() + } + if (!is.null(profile.ALL)) { + jpeg("profile.ALL.jpeg") + plotProfiles(profile.ALL, percentage=per, multiplePlots=FALSE, aTitle=title) + dev.off() + } +} + +# Plot profiles to PDF +plotPDF = function(profile.CC = NULL, profile.BP = NULL, profile.MF = NULL, profile.ALL = NULL, per = TRUE, title = TRUE) { + if (!is.null(profile.CC)) { + pdf("profile.CC.pdf") + plotProfiles(profile.CC, percentage=per, multiplePlots=FALSE, aTitle=title) + dev.off() + } + if (!is.null(profile.BP)) { + pdf("profile.BP.pdf") + plotProfiles(profile.BP, percentage=per, multiplePlots=FALSE, aTitle=title) + dev.off() + } + if (!is.null(profile.MF)) { + pdf("profile.MF.pdf") + plotProfiles(profile.MF, percentage=per, multiplePlots=FALSE, aTitle=title) + dev.off() + } + if (!is.null(profile.ALL)) { + #print("all") + pdf("profile.ALL.pdf") + plotProfiles(profile.ALL, percentage=per, multiplePlots=FALSE, aTitle=title) + dev.off() + } +} + +goprofiles = function() { + args = commandArgs(trailingOnly = TRUE) + #print(args) + # arguments: filename.R inputfile ncol "CC,MF,BP,ALL" "PNG,JPEG,PDF" level "TRUE"(percentage) "Title" + if (length(args) != 8) { + stop("Not enough/Too many arguments", call. = FALSE) + } + else { + input_type = args[2] + if (input_type == "text") { + input = strsplit(args[1], "\\s+")[[1]] + } + else if (input_type == "file") { + filename = strsplit(args[1], ",")[[1]][1] + ncol = strsplit(args[1], ",")[[1]][2] + # Check ncol + if (! as.numeric(gsub("c", "", ncol)) %% 1 == 0) { + stop("Please enter an integer for level") + } + else { + ncol = as.numeric(gsub("c", "", ncol)) + } + header = strsplit(args[1], ",")[[1]][3] + # Get file content + file = readfile(filename, header) + # Extract Protein IDs list + input = c() + for (row in as.character(file[,ncol])) { + input = c(input, strsplit(row, ";")[[1]][1]) + } + } + id_type = args[3] + ontoopt = strsplit(args[4], ",")[[1]] + #print(ontoopt) + #plotopt = strsplit(args[3], ",") + plotopt = args[5] + level = args[6] + per = as.logical(args[7]) + title = args[8] + + profiles = getprofile(input, id_type, level) + profile.CC = profiles[1] + #print(profile.CC) + profile.MF = profiles[2] + #print(profile.MF) + profile.BP = profiles[3] + #print(profile.BP) + profile.ALL = profiles[-3:-1] + #print(profile.ALL) + #c(profile.ALL, profile.CC, profile.MF, profile.BP) + if ("CC" %in% ontoopt) { + if (grepl("PNG", plotopt)) { + plotPNG(profile.CC=profile.CC, per=per, title=title) + } + if (grepl("JPEG", plotopt)) { + plotJPEG(profile.CC = profile.CC, per=per, title=title) + } + if (grepl("PDF", plotopt)) { + plotPDF(profile.CC = profile.CC, per=per, title=title) + } + } + if ("MF" %in% ontoopt) { + if (grepl("PNG", plotopt)) { + plotPNG(profile.MF = profile.MF, per=per, title=title) + } + if (grepl("JPEG", plotopt)) { + plotJPEG(profile.MF = profile.MF, per=per, title=title) + } + if (grepl("PDF", plotopt)) { + plotPDF(profile.MF = profile.MF, per=per, title=title) + } + } + if ("BP" %in% ontoopt) { + if (grepl("PNG", plotopt)) { + plotPNG(profile.BP = profile.BP, per=per, title=title) + } + if (grepl("JPEG", plotopt)) { + plotJPEG(profile.BP = profile.BP, per=per, title=title) + } + if (grepl("PDF", plotopt)) { + plotPDF(profile.BP = profile.BP, per=per, title=title) + } + } + + #if (grepl("PNG", plotopt)) { + # plotPNG(profile.ALL = profile.ALL, per=per, title=title) + #} + #if (grepl("JPEG", plotopt)) { + # plotJPEG(profile.ALL = profile.ALL, per=per, title=title) + #} + #if (grepl("PDF", plotopt)) { + # plotPDF(profile.ALL = profile.ALL, per=per, title=title) + #} + } + +} + +goprofiles() + +#Rscript go.R ../proteinGroups_Maud.txt "1" "CC" "PDF" 2 "TRUE" "Title" diff -r 000000000000 -r 6eeb2fb0c4bd goprofiles.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/goprofiles.xml Sat Sep 16 09:17:07 2017 -0400 @@ -0,0 +1,149 @@ + + Identify enriched biological themes, GO terms from your protein list. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + +
+
+ + + + + + + + + + + + + + + + + + + +
+ +
+ + + + + +
+
+ + + +
diff -r 000000000000 -r 6eeb2fb0c4bd test-data/UnipIDs.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/UnipIDs.txt Sat Sep 16 09:17:07 2017 -0400 @@ -0,0 +1,25 @@ +P04637 +P08246 +P63244 +P10275 +P00533 +Q14524 +P05067 +P35555 +P35222 +O95273 +P00451 +P38398 +Q05086 +Q12802 +P68871 +P04585 +Q96EB6 +Q9NYL2 +P31749 +P01137 +Q5S007 +Q08379 +P02649 +P35498 +P12931 diff -r 000000000000 -r 6eeb2fb0c4bd test-data/goprofile.BP.pdf Binary file test-data/goprofile.BP.pdf has changed diff -r 000000000000 -r 6eeb2fb0c4bd test-data/goprofile.CC.pdf Binary file test-data/goprofile.CC.pdf has changed diff -r 000000000000 -r 6eeb2fb0c4bd test-data/goprofile.MF.pdf Binary file test-data/goprofile.MF.pdf has changed