Mercurial > repos > iuc > ribowaltz_process

comparison ribowaltz.R @ 0:6d4c94373bba draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ribowaltz commit ff002df702f544829d1b500ac4b517c1e70ad14d
author iuc
date Thu, 22 Sep 2022 20:30:54 +0000
parents
children 042cab870a39
comparison
equal deleted inserted replaced
-1:000000000000 0:6d4c94373bba
1 # setup R error handling to go to stderr
2 options(show.error.messages = FALSE, error = function() {
3 cat(geterrmessage(), file = stderr())
4 q("no", 1, FALSE)
5 })
6
7 # we need that to not crash galaxy with an UTF8 error on German LC settings.
8 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
9 library("getopt")
10 options(stringAsFactors = FALSE, useFancyQuotes = FALSE)
11 args <- commandArgs(trailingOnly = TRUE)
12
13 # get options, using the spec as defined by the enclosed list.
14 # we read the options from the default: commandArgs(TRUE).
15 spec <- matrix(c(
16 "quiet", "q", 0, "logical",
17 "help", "h", 0, "logical",
18 "bamdir", "b", 1, "character",
19 "gtffile", "g", 1, "character",
20 "codon_coverage_info", "Y", 1, "character",
21 "cds_coverage_info", "Z", 1, "character",
22 "psite_info_rdata", "O", 0, "character",
23 "refseq_sep", "s", 0, "character",
24 "params_duplicate_filterting", "d", 0, "character",
25 "params_peridiocity_filterting", "l", 0, "character",
26 "params_custom_filterting", "c", 0, "character",
27 "params_psite_additional", "p", 0, "character",
28 "params_coverage_additional", "C", 0, "character"
29 ), byrow = TRUE, ncol = 4)
30 opt <- getopt(spec)
31
32 # if help was asked for print a friendly message
33 # and exit with a non-zero error code
34 if (!is.null(opt$help)) {
35 cat(getopt(spec, usage = TRUE))
36 q(status = 1)
37 }
38
39 verbose <- is.null(opt$quiet)
40
41 library("riboWaltz")
42
43 # create annotation data table
44 annotation_dt <- create_annotation(opt$gtffile)
45
46 sep <- opt$refseq_sep
47 if (opt$refseq_sep == "") {
48 sep <- NULL
49 }
50 # convert alignments in BAM files into list of data tables
51 reads_list <- bamtolist(bamfolder = opt$bamdir, annotation = annotation_dt, refseq_sep = sep)
52
53 library("jsonlite")
54 # remove duplicate reads
55 if (!is.null(opt$params_duplicate_filterting)) {
56 json_duplicate_filterting <- fromJSON(opt$params_duplicate_filterting)
57 reads_list <- duplicates_filter(
58 data = reads_list,
59 extremity = json_duplicate_filterting$extremity,
60 keep = json_duplicate_filterting$keep
61 )
62 }
63
64 # selection of read lengths - periodicity filtering
65 if (!is.null(opt$params_peridiocity_filterting)) {
66 json_peridiocity_filterting <- fromJSON(opt$params_peridiocity_filterting)
67 reads_list <- length_filter(
68 data = reads_list,
69 length_filter_mode = "periodicity",
70 periodicity_threshold = json_peridiocity_filterting$periodicity_threshold
71 )
72 }
73 # selection of read lengths - length range filtering
74 if (!is.null(opt$params_custom_filterting)) {
75 json_custom_filterting <- fromJSON(opt$params_custom_filterting)
76 reads_list <- length_filter(
77 data = reads_list,
78 length_filter_mode = "custom",
79 length_range = json_custom_filterting$length_range
80 )
81 }
82
83 # compute P-site offset
84 json_psite_additional <- fromJSON(opt$params_psite_additional)
85 psite_offset <- psite(
86 reads_list,
87 start = json_psite_additional$use_start,
88 flanking = json_psite_additional$flanking,
89 extremity = json_psite_additional$psite_extrimity,
90 plot = TRUE,
91 cl = json_psite_additional$cl,
92 plot_format = "pdf",
93 plot_dir = "plots"
94 )
95 psite_offset
96
97 reads_psite_list <- psite_info(reads_list, psite_offset)
98 reads_psite_list
99 # write a separate P-site offset info table for each sample
100 for (sample in names(reads_psite_list)) {
101 write.table(
102 reads_psite_list[[sample]],
103 file = paste(sample, "psite_info.tsv", sep = "_"),
104 sep = "\t", row.names = FALSE, quote = FALSE
105 )
106 print(paste(sample, "psite_info.tsv", sep = "_"))
107 }
108
109 # write R object to a file
110 if (!is.null(opt$psite_info_rdata)) {
111 save(reads_psite_list, annotation_dt, file = opt$psite_info_rdata)
112 }
113
114 json_coverage_additional <- fromJSON(opt$params_coverage_additional)
115 # codon coverage
116 codon_coverage_out <- codon_coverage(
117 reads_psite_list,
118 annotation_dt,
119 psite = json_coverage_additional$psites_per_region,
120 min_overlap = json_coverage_additional$min_overlap
121 )
122 write.table(codon_coverage_out, file = opt$codon_coverage_info, sep = "\t", row.names = FALSE, quote = FALSE)
123
124 # CDS coverage
125 cds_coverage_out <- cds_coverage(
126 reads_psite_list,
127 annotation_dt,
128 start_nts = json_coverage_additional$start_nts,
129 stop_nts = json_coverage_additional$stop_nts
130 )
131 write.table(cds_coverage_out, file = opt$cds_coverage_info, sep = "\t", row.names = FALSE, quote = FALSE)