Mercurial > repos > galaxyp > custom_pro_db_annotation_data_manager

comparison data_manager/customProDB_annotation.R @ 1:9b4ee836e35b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tools/bumbershoot/custom_pro_db commit 2174137cf8a15deefed5910ffa152c4ce9c81af6
author galaxyp
date Thu, 08 Jun 2017 10:55:08 -0400
parents 45755942ae7b
children 0a9ffebba65d
comparison
equal deleted inserted replaced
0:45755942ae7b 1:9b4ee836e35b
10 10
11 options('useFancyQuotes' = FALSE) 11 options('useFancyQuotes' = FALSE)
12 12
13 suppressPackageStartupMessages(library("optparse")) 13 suppressPackageStartupMessages(library("optparse"))
14 suppressPackageStartupMessages(library("RGalaxy")) 14 suppressPackageStartupMessages(library("RGalaxy"))
15 suppressPackageStartupMessages(library("GetoptLong"))
15 16
16 17
17 option_list <- list() 18 option_list <- list()
18 option_list$dbkey <- make_option('--dbkey', type='character') 19 option_list$dbkey <- make_option('--dbkey', type='character')
20 option_list$ensembl_host <- make_option('--ensembl_host', type='character')
21 option_list$ensembl_dataset <- make_option('--ensembl_dataset', type='character')
19 option_list$dbsnp <- make_option('--dbsnp', type='character') 22 option_list$dbsnp <- make_option('--dbsnp', type='character')
20 option_list$cosmic <- make_option('--cosmic', type='logical') 23 option_list$cosmic <- make_option('--cosmic', type='logical')
21 option_list$outputFile <- make_option('--outputFile', type='character') 24 option_list$outputFile <- make_option('--outputFile', type='character')
22 option_list$dbkey_description <- make_option('--dbkey_description', type='character') 25 option_list$dbkey_description <- make_option('--dbkey_description', type='character')
23 26
24 opt <- parse_args(OptionParser(option_list=option_list)) 27 opt <- parse_args(OptionParser(option_list=option_list))
25 28
26 29
27 customProDB_annotation <- function( 30 customProDB_annotation <- function(
28 dbkey = GalaxyCharacterParam(required=TRUE), 31 dbkey = GalaxyCharacterParam(required=FALSE),
32 ensembl_host = GalaxyCharacterParam(required=FALSE),
33 ensembl_dataset = GalaxyCharacterParam(required=FALSE),
29 dbsnp_str = GalaxyCharacterParam(required=FALSE), 34 dbsnp_str = GalaxyCharacterParam(required=FALSE),
30 cosmic = GalaxyLogicalParam(required=FALSE), 35 cosmic = GalaxyLogicalParam(required=FALSE),
31 dbkey_description = GalaxyCharacterParam(required=FALSE), 36 dbkey_description = GalaxyCharacterParam(required=FALSE),
32 outputFile = GalaxyOutput("output","json")) 37 outputFile = GalaxyOutput("output","json"))
33 { 38 {
39 options(stringsAsFactors = FALSE, gsubfn.engine = "R")
40
34 if (!file.exists(outputFile)) 41 if (!file.exists(outputFile))
35 { 42 {
36 gstop("json params file does not exist") 43 gstop("json params file does not exist")
37 } 44 }
38 45
39 if (length(dbkey_description) < 1) 46 if (length(dbkey)+length(ensembl_dataset)+length(ensembl_host) == 0)
40 { 47 {
41 dbkey_description = dbkey 48 gstop("one of the genome annotation sources must be specified; either dbkey or host and dataset")
49 }
50 else if (length(dbkey) > 0 &&
51 (length(ensembl_dataset) > 0 || length(ensembl_host) > 0))
52 {
53 gstop("only one genome annotation source can be specified; either dbkey or host and dataset")
42 } 54 }
43 55
44 if (length(dbsnp_str) > 0) 56 if (length(dbsnp_str) > 0)
45 { 57 {
46 dbsnp = dbsnp_str 58 dbsnp = dbsnp_str
51 } 63 }
52 64
53 use_cosmic = FALSE 65 use_cosmic = FALSE
54 if (length(cosmic) > 0) 66 if (length(cosmic) > 0)
55 { 67 {
56 if (grepl("^hg", dbkey)) 68 if (length(dbkey) > 0 && grepl("^hg", dbkey) ||
69 length(ensembl_dataset) > 0 && grepl("^hsapiens", ensembl_dataset))
57 { 70 {
58 use_cosmic = TRUE 71 use_cosmic = TRUE
59 } 72 }
60 else 73 else
61 { 74 {
74 }, error=function(err) 87 }, error=function(err)
75 { 88 {
76 gstop("failed to remove json params file after reading") 89 gstop("failed to remove json params file after reading")
77 }) 90 })
78 91
79 ucscTableCodingFastaURL = paste("http://genome.ucsc.edu/cgi-bin/hgTables?db=", dbkey, "&hgSeq.cdsExon=on&hgSeq.granularity=gene&hgSeq.casing=exon&hgSeq.repMasking=lower&hgta_doGenomicDna=get+sequence&hgta_group=genes&hgta_track=refGene&hgta_table=refGene&hgta_regionType=genome", sep="") 92 # load customProDB from GitHub (NOTE: downloading the zip is faster than cloning the repo with git2r or devtools::install_github)
80 ucscTableProteinFastaURL = paste("http://genome.ucsc.edu/cgi-bin/hgTables?db=", dbkey, "&hgta_geneSeqType=protein&hgta_doGenePredSequence=submit&hgta_track=refGene&hgta_table=refGene", sep="") 93 download.file("https://github.com/chambm/customProDB/archive/c57e5498392197bc598a18c26acb70d7530a921cc57e5498.zip", "customProDB.zip", quiet=TRUE)
81 codingFastaFilepath = paste(target_directory, "/", dbkey, ".cds.fa", sep="") 94 unzip("customProDB.zip")
82 proteinFastaFilepath = paste(target_directory, "/", dbkey, ".protein.fa", sep="") 95 devtools::load_all("customProDB-c57e5498392197bc598a18c26acb70d7530a921c")
83 96
84 suppressPackageStartupMessages(library(customProDB)) 97 #suppressPackageStartupMessages(library(customProDB))
85 options(timeout=3600) 98 options(timeout=3600)
86 99
87 cat(paste("Downloading coding FASTA from:", ucscTableCodingFastaURL, "\n")) 100 # download protein and coding sequences for UCSC annotation
88 download.file(ucscTableCodingFastaURL, codingFastaFilepath, quiet=T, mode='wb') 101 if (length(dbkey) > 0)
89 102 {
90 cat(paste("Downloading protein FASTA from:", ucscTableProteinFastaURL, "\n")) 103 proteinFastaFilepath = paste(dbkey, ".protein.fa", sep="")
91 download.file(ucscTableProteinFastaURL, proteinFastaFilepath, quiet=T, mode='wb') 104
92 105 cat(paste("Downloading protein FASTA from:", getProteinFastaUrlFromUCSC(dbkey), "\n"))
93 cat(paste("Preparing Refseq annotation files\n")) 106 download.file(getProteinFastaUrlFromUCSC(dbkey), proteinFastaFilepath, quiet=T, mode='wb')
94 customProDB::PrepareAnnotationRefseq(genome=dbkey, CDSfasta=codingFastaFilepath, pepfasta=proteinFastaFilepath, annotation_path=target_directory, dbsnp=dbsnp, COSMIC=use_cosmic) 107
95 108 local_cache_path = paste0("customProDB_annotation_", dbkey, "-", tools::md5sum(proteinFastaFilepath)[[1]])
96 outputPath = paste(dbkey, "/customProDB", sep="") 109 codingFastaFilepath = paste0(local_cache_path, "/", dbkey, ".cds.fa")
110 dir.create(local_cache_path, showWarnings=FALSE)
111
112 if (!file.exists(codingFastaFilepath)) {
113 cat(paste("Downloading coding FASTA from:", getCodingFastaUrlFromUCSC(dbkey), "\n"))
114 download.file(getCodingFastaUrlFromUCSC(dbkey), codingFastaFilepath, quiet=T, mode='wb')
115 }
116
117 cat(paste("Preparing Refseq annotation files\n"))
118 PrepareAnnotationRefseq(genome=dbkey, CDSfasta=codingFastaFilepath, pepfasta=proteinFastaFilepath, annotation_path=target_directory, dbsnp=dbsnp, COSMIC=use_cosmic, local_cache_path=local_cache_path)
119
120 if (length(dbkey_description) < 1)
121 {
122 dbkey_description = dbkey
123 }
124 }
125 else
126 {
127 local_cache_path = paste0("customProDB_annotation_", ensembl_dataset, "_", ensembl_host)
128
129 suppressPackageStartupMessages(library(biomaRt))
130 cat(paste("Preparing Ensembl annotation files\n"))
131 ensembl_mart = useMart("ENSEMBL_MART_ENSEMBL", dataset=ensembl_dataset, host=ensembl_host)
132 PrepareAnnotationEnsembl(mart=ensembl_mart, annotation_path=target_directory, dbsnp=dbsnp, COSMIC=use_cosmic, local_cache_path=local_cache_path)
133
134 metadata = sqldf::sqldf("SELECT value FROM metadata WHERE name='BioMart database version' OR name='BioMart dataset description' OR name='BioMart dataset version'",
135 dbname=file.path(target_directory, "txdb.sqlite"))
136 version = metadata$value[1] # Ensembl Genes 87
137 assembly = metadata$value[3]
138 dbkey = paste0(ensembl_dataset, "_", sub(".*?(\\d+)", "\\1", version, perl=TRUE))
139
140 # convert Ensembl chromosome names to UCSC for Galaxy compatibility
141 chromosomeMappingsBaseUrl = "https://raw.githubusercontent.com/dpryan79/ChromosomeMappings/master"
142 assemblyNoGrcPatch = sub("(\\S+?)(\\.p\\S+)?$", "\\1", assembly, perl=TRUE)
143 chromosomeMappingsUrl = qq("@{chromosomeMappingsBaseUrl}/@{assemblyNoGrcPatch}_ensembl2UCSC.txt")
144 if (RCurl::url.exists(chromosomeMappingsUrl))
145 {
146 cat(qq("Converting Ensembl chromosome names from: @{chromosomeMappingsUrl}\n"))
147 e2u = read.delim(chromosomeMappingsUrl, header=FALSE, col.names=c("ensembl", "ucsc"))
148 e2u = setNames(as.list(e2u$ucsc), e2u$ensembl)
149 load(file.path(target_directory, "exon_anno.RData"))
150 exon$chromosome_name = sapply(exon$chromosome_name, function(x) e2u[[as.character(x)]])
151 exon = exon[nzchar(exon$chromosome_name), ] # omit genome patches with no mapping
152 save(exon, file=file.path(target_directory, "exon_anno.RData"))
153 }
154 else
155 {
156 gwarning(qq("unable to convert Ensembl chromosome names to UCSC; mapping file @{assemblyNoGrcPatch}_ensembl2UCSC.txt does not exist"))
157 }
158
159 if (length(dbkey_description) < 1)
160 {
161 dbkey_description = qq("@{ensembl_dataset} (@{version}) (@{assembly})")
162 }
163 }
164
165 qualified_dbkey = dbkey
166
167 if (length(dbsnp_str) > 0 && nzchar(dbsnp_str))
168 {
169 qualified_dbkey = qq("@{qualified_dbkey}_db@{dbsnp_str}")
170 dbkey_description = qq("@{dbkey_description} (db@{dbsnp_str})")
171 }
172
173 if (length(cosmic) > 0)
174 {
175 qualified_dbkey = qq("@{qualified_dbkey}_cosmic")
176 dbkey_description = qq("@{dbkey_description} (COSMIC)")
177 }
178
179 outputPath = paste0(qualified_dbkey, "/customProDB")
97 output = list(data_tables = list()) 180 output = list(data_tables = list())
98 output[["data_tables"]][["customProDB"]]=c(path=outputPath, name=dbkey_description, dbkey=dbkey, value=dbkey) 181 output[["data_tables"]][["customProDB"]]=c(path=outputPath, name=dbkey_description, dbkey=qualified_dbkey, value=qualified_dbkey)
99 write(toJSON(output), file=outputFile) 182 write(toJSON(output), file=outputFile)
100 } 183 }
101 184
102 185
103 params <- list() 186 params <- list()