Mercurial > repos > galaxyp > psm_to_sam

comparison PSM2SAM.R @ 0:90ecb65017a0 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/psm2sam commit 9910ed076e4b8a3f083351b89fa861d0a4a93beb
author galaxyp
date Wed, 17 May 2017 20:23:27 -0400
parents
children 757ecf27b4a9
comparison
equal deleted inserted replaced
-1:000000000000 0:90ecb65017a0
1 #!/usr/bin/env Rscript
2
3 ## begin warning handler
4 withCallingHandlers({
5
6 library(methods) # Because Rscript does not always do this
7
8 options('useFancyQuotes' = FALSE)
9
10 suppressPackageStartupMessages(library("optparse"))
11 suppressPackageStartupMessages(library("RGalaxy"))
12
13
14 option_list <- list()
15
16 option_list$exon_anno <- make_option('--exon_anno', type='character')
17 option_list$proteinseq <- make_option('--proteinseq', type='character')
18 option_list$procodingseq <- make_option('--procodingseq', type='character')
19 option_list$bam_file <- make_option('--bam', type='character')
20 option_list$idpDB_file <- make_option('--idpDB', type='character')
21 option_list$pepXmlTab_file <- make_option('--pepXmlTab', type='character')
22 option_list$peptideShakerPsmReport_file <- make_option('--peptideShakerPsmReport', type='character')
23 option_list$variantAnnotation_file <- make_option('--variantAnnotation', type='character')
24 option_list$searchEngineScore <- make_option('--searchEngineScore', type='character')
25
26
27 opt <- parse_args(OptionParser(option_list=option_list))
28
29
30 psm2sam <- function(
31 exon_anno_file = GalaxyInputFile(required=TRUE),
32 proteinseq_file = GalaxyInputFile(required=TRUE),
33 procodingseq_file = GalaxyInputFile(required=TRUE),
34 bam_file = GalaxyInputFile(required=TRUE),
35 idpDB_file = GalaxyInputFile(required=FALSE),
36 pepXmlTab_file = GalaxyInputFile(required=FALSE),
37 peptideShakerPsmReport_file = GalaxyInputFile(required=FALSE),
38 variantAnnotation_file = GalaxyInputFile(required=FALSE),
39 searchEngineScore = GalaxyCharacterParam(required=FALSE)
40 )
41 {
42 options(stringsAsFactors = FALSE)
43
44 if (length(bam_file) == 0)
45 {
46 stop("BAM file must be specified to provide sequence headers")
47 }
48
49 outputHeader = grep("^@(?!PG)", readLines(bam_file, n=500, warn=FALSE), value=TRUE, perl=TRUE)
50 if (length(outputHeader) == 0)
51 {
52 stop("failed to read header lines from bam_file")
53 }
54
55 # load customProDB from GitHub (NOTE: downloading the zip is faster than cloning the repo with git2r or devtools::install_github)
56 download.file("https://github.com/chambm/customProDB/archive/master.zip", "customProDB.zip", quiet=TRUE)
57 unzip("customProDB.zip")
58 devtools::load_all("customProDB-master")
59
60 # load proBAMr from GitHub
61 download.file("https://github.com/chambm/proBAMr/archive/master.zip", "proBAMr.zip", quiet=TRUE)
62 unzip("proBAMr.zip")
63 devtools::load_all("proBAMr-master")
64
65 psmInputLength = length(idpDB_file)+length(pepXmlTab_file)+length(peptideShakerPsmReport_file)
66 if (psmInputLength == 0)
67 {
68 stop("one of the input PSM file parameters must be specified")
69 }
70 else if (psmInputLength > 1)
71 {
72 stop("only one of the input PSM file parameters can be specified")
73 }
74
75 if (length(idpDB_file) > 0)
76 {
77 if (length(searchEngineScore) == 0)
78 stop("searchEngineScore parameter must be specified when reading IDPicker PSMs, e.g. 'MyriMatch:MVH'")
79 passedPSM = readIdpDB(idpDB_file, searchEngineScore)
80 }
81 else if (length(pepXmlTab_file) > 0)
82 {
83 if (length(searchEngineScore) == 0)
84 stop("searchEngineScore parameter must be specified when reading pepXmlTab PSMs, e.g. 'mvh'")
85 passedPSM = readPepXmlTab(pepXmlTab_file, searchEngineScore)
86 }
87 else if (length(peptideShakerPsmReport_file) > 0)
88 {
89 if (length(searchEngineScore) > 0)
90 warning("searchEngineScore parameter is ignored when reading PeptideShaker PSM report")
91 passedPSM = readPeptideShakerPsmReport(peptideShakerPsmReport_file)
92 }
93
94 load(exon_anno_file)
95 load(proteinseq_file)
96 load(procodingseq_file)
97
98 if (length(variantAnnotation_file) > 0)
99 {
100 load(variantAnnotation_file) # variantAnnotation list, with members snvprocoding/snvproseq and indelprocoding/indelproseq
101
102 varprocoding = unique(rbind(variantAnnotation$snvprocoding, variantAnnotation$indelprocoding))
103 varproseq = unique(rbind(variantAnnotation$snvproseq, variantAnnotation$indelproseq))
104 }
105 else
106 {
107 varprocoding = NULL
108 varproseq = NULL
109 }
110
111 # add proBAMr program key
112 outputHeader = c(outputHeader, paste0("@PG\tID:proBAMr\tVN:", packageVersion("proBAMr")))
113
114 # first write header lines to the output SAM
115 writeLines(outputHeader, "output.sam")
116
117 # then write the PSM "reads"
118 PSMtab2SAM(passedPSM, exon,
119 proteinseq, procodingseq,
120 varproseq, varprocoding,
121 outfile = "output.sam",
122 show_progress = FALSE)
123
124 invisible(NULL)
125 }
126
127 params <- list()
128 for(param in names(opt))
129 {
130 if (!param == "help")
131 params[param] <- opt[param]
132 }
133
134 setClass("GalaxyRemoteError", contains="character")
135 wrappedFunction <- function(f)
136 {
137 tryCatch(do.call(f, params),
138 error=function(e) new("GalaxyRemoteError", conditionMessage(e)))
139 }
140
141
142 suppressPackageStartupMessages(library(RGalaxy))
143 do.call(psm2sam, params)
144
145 ## end warning handler
146 }, warning = function(w) {
147 cat(paste("Warning:", conditionMessage(w), "\n"))
148 invokeRestart("muffleWarning")
149 })