annotate id_converter_UniProt.R @ 14:659f1248f535 draft

planemo upload commit 688c456ca57914a63c20eba942ec5fe81e896099-dirty

# Read file and return file content as data.frame
readfile = function(filename, header) {
  if (header == "true") {
    # Read only first line of the file as header:
    headers <- read.table(filename, nrows = 1, header = FALSE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE, na.strings=c("", "NA"), blank.lines.skip = TRUE, quote = "")
    #Read the data of the files (skipping the first row)
    file <- read.table(filename, skip = 1, header = FALSE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE, na.strings=c("", "NA"), blank.lines.skip = TRUE, quote = "")
    # Remove empty rows
    file <- file[!apply(is.na(file) | file == "", 1, all), , drop=FALSE]
    #And assign the header to the data
    names(file) <- headers
  }
  else {
    file <- read.table(filename, header = FALSE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE, na.strings=c("", "NA"), blank.lines.skip = TRUE, quote = "")
    # Remove empty rows
    file <- file[!apply(is.na(file) | file == "", 1, all), , drop=FALSE]
  }
  return(file)
}

# Mapping IDs using file built from
#   - HUMAN_9606_idmapping_selected.tab
#     Tarball downloaded from ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/idmapping/by_organism/
#   - nextprot_ac_list_all.txt 
#     Downloaded from ftp://ftp.nextprot.org/pub/current_release/ac_lists/
# Available databases: 
#   UNIPROT_AC: Uniprot accession number (e.g. P31946)
#   UNIPROT_ID: Uniprot identifiers (e.g 1433B_HUMAN)
#   GeneID_EntrezGene: Entrez gene ID (serie of digit) (e.g. 7529)
#   RefSeq: RefSeq (NCBI) protein (e.g.  NP_003395.1; NP_647539.1; XP_016883528.1)
#   GI_number: GI (NCBI GI number) ID (serie of digits) assigned to each sequence record processed by NCBI (e.g; 21328448; 377656701; 67464627; 78101741) 
#   PDB: Protein DataBank Identifiers (e.g. 2BR9:A; 3UAL:A;   3UBW:A) 
#   GO_ID: GOterms (Gene Ontology) Identifiers (e.g. GO:0070062; GO:0005925; GO:0042470; GO:0016020; GO:0005739; GO:0005634)
#   PIR: Protein Information Resource ID (e.g. S34755)	
#   OMIM: OMIM (Online Mendelian Inheritance in Man database) ID (serie of digits) (e.g: 601289)	
#   UniGene: Unigene Identifier (e.g. Hs.643544)
#   Ensembl_ENSG: Ensembl gene identifiers (e.g. ENSG00000166913) 
#   Ensembl_ENST: Ensembl transcript identifiers (e.g. ENST00000353703; ENST00000372839)
#   Ensembl_ENSP: Ensembl protein identifiers (e.g. ENSP00000300161; ENSP00000361930)

mapping = function() {
  args <- commandArgs(TRUE)
  if(length(args)<1) {
    args <- c("--help")
  }
  
  # Help section
  if("--help" %in% args) {
    cat("Selection and Annotation HPA
    Arguments:
        --ref_file: path to reference file (id_mapping_file.txt)
        --input_type: type of input (list of id or filename)
        --id_type: type of input IDs
        --input: list of IDs (text or filename)
        --column_number: the column number which contains list of input IDs
        --header: true/false if your file contains a header
        --target_ids: target IDs to map to 
        --output: output filename \n")
    q(save="no")
  }
  
  # Parse arguments
  parseArgs <- function(x) strsplit(sub("^--", "", x), "=")
  argsDF <- as.data.frame(do.call("rbind", parseArgs(args)))
  args <- as.list(as.character(argsDF$V2))
  names(args) <- argsDF$V1
  
  input_id_type = args$id_type # Uniprot, ENSG....
  list_id_input_type = args$input_type # list or file
  options = strsplit(args$target_ids, ",")[[1]]
  output = args$output
  id_mapping_file = args$ref_file
    
  # Extract input IDs
  if (list_id_input_type == "list") {
    print(args$input)
    list_id = trimws(strsplit(args$input, ",")[[1]])
    list_id = list_id[list_id != ""]    #remove empty entry
    # Remove isoform accession number (e.g. "-2")
    list_id = gsub("-.+", "", list_id)
  } else if (list_id_input_type == "file") {
    filename = args$input
    column_number = as.numeric(gsub("c", "" ,args$column_number))
    header = args$header
    file_all = readfile(filename, header)
    list_id = trimws(gsub("[$,\xc2\xa0]","",sapply(strsplit(file_all[,column_number], ";"), "[", 1)))
    # Remove isoform accession number (e.g. "-2")
    list_id = gsub("-.+", "", list_id)
  }

  # Extract ID maps
  id_map = read.table(id_mapping_file, header = TRUE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE, na.strings = "", quote = "")
  
  names = c()
    
  # Map IDs
  res = matrix(nrow=length(list_id), ncol=0)

  for (opt in options) {
    names = c(names, opt)
    mapped = id_map[match(list_id, id_map[input_id_type][,]),][opt][,]
    res = cbind(res, matrix(mapped))
  }
     
  # Write output
  if (list_id_input_type == "list") {
    res = cbind(as.matrix(list_id), res)
    names = c(input_id_type, names)
    colnames(res) = names
    write.table(res, output, row.names = FALSE, sep = "\t", quote = FALSE)
  }
  else if (list_id_input_type == "file") {
    names(res) = options
    if (all(names(file_all) == file_all[1,1:length(names(file_all))])){ #if header of file is the same as the first line of file
      names(file_all)[column_number] = input_id_type
    }
    names = c(names(file_all), names)
    output_content = cbind(file_all, res)
    colnames(output_content) = names
    write.table(output_content, output, row.names = FALSE, sep = "\t", quote = FALSE)
  }
}

mapping()

#Rscript id_converter_UniProt.R "UniProt.AC" "test-data/UnipIDs.txt,c1,false" "file" "Ensembl_PRO,Ensembl,neXtProt_ID" "test-data/output.txt" ../../utils/mapping_file.txt