annotate get_expression_profiles.R @ 13:2173192d4824 draft default tip

"planemo upload commit 59b014e9f6e2d668cbd7c4844b10db3d59baefd8"

1 # Read file and return file content as data.frame
2 read_file <- function(path,header){
3   file <- try(read.csv(path,header=header, sep="\t",stringsAsFactors = FALSE, quote="\"", check.names = F),silent=TRUE)
4   if (inherits(file,"try-error")){
5     stop("File not found !")
6   }else{
7     return(file)
8   }
9 }
10 
11 str2bool <- function(x){
12   if (any(is.element(c("t","true"),tolower(x)))){
13     return (TRUE)
14   }else if (any(is.element(c("f","false"),tolower(x)))){
15     return (FALSE)
16   }else{
17     return(NULL)
18   }
19 }
20 
21 # input has to be a list of IDs in ENSG format
22 # tissue is one of unique(HPA.normal.tissue$Tissue)
23 # level is one, or several, or 0 (=ALL) of "Not detected", "Medium", "High", "Low"
24 # reliability is one, or several, or 0 (=ALL) of "Approved", "Supported", "Uncertain"
25 annot.HPAnorm<-function(input, HPA_normal_tissue, tissue, level, reliability, not_mapped_option) {
26   dat <- subset(HPA_normal_tissue, Gene %in% input)
27   res.Tissue<-subset(dat, Tissue %in% tissue) 
28   res.Level<-subset(res.Tissue, Level %in% level) 
29   res.Rel<-subset(res.Level, Reliability %in% reliability) 
30   
31   if (not_mapped_option) {
32     if (length(setdiff(intersect(input, unique(dat$Gene)), unique(res.Rel$Gene)))>0) {
33       not_match_IDs <- matrix(setdiff(intersect(input, unique(dat$Gene)), unique(res.Rel$Gene)), ncol = 1, nrow = length(setdiff(intersect(input, unique(dat$Gene)), unique(res.Rel$Gene))))
34       not.match <- matrix("no match", ncol = ncol(HPA_normal_tissue) - 1, nrow = length(not_match_IDs))
35       not.match <- cbind(not_match_IDs, unname(not.match))
36       colnames(not.match) <- colnames(HPA_normal_tissue)
37       res <- rbind(res.Rel, not.match)
38     } else {
39       res <- res.Rel
40     } 
41     
42     if (length(setdiff(input, unique(dat$Gene)))>0) {
43       not.mapped <- matrix(ncol = ncol(HPA_normal_tissue) - 1, nrow = length(setdiff(input, unique(dat$Gene))))
44       not.mapped <- cbind(matrix(setdiff(input, unique(dat$Gene)), ncol = 1, nrow = length(setdiff(input, unique(dat$Gene)))), unname(not.mapped))
45       colnames(not.mapped) <- colnames(HPA_normal_tissue)
46       res <- rbind(res, not.mapped)
47     }
48     
49   } else {
50     res <- res.Rel
51   }
52   
53   return(res)
54   
55 }
56 
57 annot.HPAcancer<-function(input, HPA_cancer_tissue, cancer, not_mapped_option) {
58   dat <- subset(HPA_cancer_tissue, Gene %in% input)
59   res.Cancer<-subset(dat, Cancer %in% cancer) 
60 
61   if (not_mapped_option) {
62     not.mapped <- matrix(ncol=ncol(HPA_cancer_tissue)-1, nrow=length(setdiff(input, unique(dat$Gene))))
63     not.mapped <- cbind(matrix(setdiff(input, unique(dat$Gene)), ncol = 1, nrow = length(setdiff(input, unique(dat$Gene)))), unname(not.mapped))
64     colnames(not.mapped) <- colnames(HPA_cancer_tissue)
65     res <- rbind(res.Cancer, not.mapped)
66   } else {
67     res <- res.Cancer
68   }
69   return(res)
70 }
71 
72 clean_ids <- function(ids){
73   
74   ids = gsub(" ","",ids)
75   ids = ids[which(ids!="")]
76   ids = ids[which(ids!="NA")]
77   ids = ids[!is.na(ids)]
78  
79   return(ids) 
80 }
81 
82 main <- function() {
83   args <- commandArgs(TRUE)
84   if(length(args)<1) {
85     args <- c("--help")
86   }
87   
88   # Help section
89   if("--help" %in% args) {
90     cat("Selection and Annotation HPA
91     Arguments:
92         --ref_file: HPA normal/cancer tissue file path
93         --input_type: type of input (list of id or filename)
94         --input: list of IDs in ENSG format
95         --column_number: the column number which you would like to apply...
96         --header: true/false if your file contains a header
97         --atlas: normal/cancer
98           if normal:
99             --tissue: list of tissues
00             --level: Not detected, Low, Medium, High
01             --reliability: Supportive, Uncertain
02           if cancer:
03             --cancer: Cancer tissues
04         --not_mapped: true/false if your output file should contain not-mapped and not-match IDs 
05         --output: output filename \n")
06     q(save="no")
07   }
08   
09   # Parse arguments
10   parseArgs <- function(x) strsplit(sub("^--", "", x), "=")
11   argsDF <- as.data.frame(do.call("rbind", parseArgs(args)))
12   args <- as.list(as.character(argsDF$V2))
13   names(args) <- argsDF$V1
14   
15   #save(args,file = "/home/dchristiany/proteore_project/ProteoRE/tools/Get_expression_profiles/args.rda")
16   #load("/home/dchristiany/proteore_project/ProteoRE/tools/Get_expression_profiles/args.rda")
17   
18   # Extract input
19   input_type = args$input_type
20   if (input_type == "list") {
21     ids = unlist(strsplit(strsplit(args$input, "[ \t\n]+")[[1]],";"))
22   } else if (input_type == "file") {
23     filename = args$input
24     column_number = as.numeric(gsub("c", "" ,args$column_number))
25     header = str2bool(args$header)
26     file = read_file(filename, header)
27     ids = unlist(strsplit(file[,column_number], ";"))
28   }
29   #filter ids
30   ids = clean_ids(ids)
31   
32   # Read reference file
33   reference_file = read_file(args$ref_file, TRUE)
34 
35   # Extract other options
36   atlas = args$atlas
37   not_mapped_option = str2bool(args$not_mapped)
38   if (atlas=="normal") {
39     tissue = strsplit(args$tissue, ",")[[1]]
40     level = strsplit(args$level, ",")[[1]]
41     reliability = strsplit(args$reliability, ",")[[1]]
42     # Calculation
43     res = annot.HPAnorm(ids, reference_file, tissue, level, reliability, not_mapped_option)
44   } else if (atlas=="cancer") {
45     cancer = strsplit(args$cancer, ",")[[1]]
46     # Calculation
47     res = annot.HPAcancer(ids, reference_file, cancer, not_mapped_option)
48   }
49   
50   # Write output
51   output = args$output
52   res <- apply(res, c(1,2), function(x) gsub("^$|^ $", NA, x))
53   write.table(res, output, sep = "\t", quote = FALSE, row.names = FALSE)
54 }
55 
56 main()