Mercurial > repos > artbio > small_rna_maps

comparison small_rna_maps.r @ 34:966bc5c46efd draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_maps commit b7c7c60d608694ca4f1638e4bb0d6db5b1afa944
author artbio
date Fri, 21 Oct 2022 23:19:55 +0000
parents f2e7ad3058e8
children
comparison
equal deleted inserted replaced
33:6fee74f6d2bd 34:966bc5c46efd
1 ## Setup R error handling to go to stderr 1 ## Setup R error handling to go to stderr
2 options(show.error.messages = F, 2 options(show.error.messages = FALSE,
3 error = function() { 3 error = function() {
4 cat(geterrmessage(), file = stderr()); q("no", 1, F) 4 cat(geterrmessage(), file = stderr())
5 q("no", 1, FALSE)
5 } 6 }
6 ) 7 )
7 options(warn = -1) 8 options(warn = -1)
8 library(RColorBrewer) 9 library(RColorBrewer)
9 library(lattice) 10 library(lattice)
29 args <- parse_args(parser) 30 args <- parse_args(parser)
30 31
31 # data frames implementation 32 # data frames implementation
32 33
33 ## first table 34 ## first table
34 table <- read.delim(args$first_dataframe, header = T, row.names = NULL) 35 table <- read.delim(args$first_dataframe, header = TRUE, row.names = NULL)
35 colnames(table)[1] <- "Dataset" 36 colnames(table)[1] <- "Dataset"
36 dropcol <- c("Strandness", "z.score") # not used by this Rscript and is dropped for backward compatibility 37 dropcol <- c("Strandness", "z.score") # not used by this Rscript and is dropped for backward compatibility
37 table <- table[, !(names(table) %in% dropcol)] 38 table <- table[, !(names(table) %in% dropcol)]
38 if (args$first_plot_method == "Counts" | args$first_plot_method == "Size") { 39 if (args$first_plot_method == "Counts" || args$first_plot_method == "Size") {
39 table <- within(table, Counts[Polarity == "R"] <- (Counts[Polarity == "R"] * - 1)) 40 table <- within(table, Counts[Polarity == "R"] <- (Counts[Polarity == "R"] * - 1)) # nolint
40 } 41 }
41 n_samples <- length(unique(table$Dataset)) 42 n_samples <- length(unique(table$Dataset))
42 samples <- unique(table$Dataset) 43 samples <- unique(table$Dataset)
43 if (args$normalization != "") { 44 if (args$normalization != "") {
44 norm_factors <- as.numeric(unlist(strsplit(args$normalization, " "))) 45 norm_factors <- as.numeric(unlist(strsplit(args$normalization, " ")))
45 } else { 46 } else {
46 norm_factors <- rep(1, n_samples) 47 norm_factors <- rep(1, n_samples)
47 } 48 }
48 if (args$first_plot_method == "Counts" | args$first_plot_method == "Size" | args$first_plot_method == "Coverage") { 49 if (args$first_plot_method == "Counts" || args$first_plot_method == "Size" || args$first_plot_method == "Coverage") {
49 i <- 1 50 i <- 1
50 for (sample in samples) { 51 for (sample in samples) {
51 # Warning Here the column is hard coded as the last column (dangerous) 52 # Warning Here the column is hard coded as the last column (dangerous)
52 # because its name changes with the method 53 # because its name changes with the method
53 table[, length(table)][table$Dataset == sample] <- table[, length(table)][table$Dataset == sample] * norm_factors[i] 54 table[, length(table)][table$Dataset == sample] <- table[, length(table)][table$Dataset == sample] * norm_factors[i]
59 per_gene_limit <- lapply(genes, function(x) c(1, unique(subset(table, Chromosome == x)$Chrom_length))) 60 per_gene_limit <- lapply(genes, function(x) c(1, unique(subset(table, Chromosome == x)$Chrom_length)))
60 n_genes <- length(per_gene_readmap) 61 n_genes <- length(per_gene_readmap)
61 62
62 # second table 63 # second table
63 if (args$extra_plot_method != "") { 64 if (args$extra_plot_method != "") {
64 extra_table <- read.delim(args$extra_dataframe, header = T, row.names = NULL) 65 extra_table <- read.delim(args$extra_dataframe, header = TRUE, row.names = NULL)
65 colnames(extra_table)[1] <- "Dataset" 66 colnames(extra_table)[1] <- "Dataset"
66 dropcol <- c("Strandness", "z.score") 67 dropcol <- c("Strandness", "z.score")
67 table <- table[, !(names(table) %in% dropcol)] 68 table <- table[, !(names(table) %in% dropcol)]
68 if (args$extra_plot_method == "Counts" | args$extra_plot_method == "Size") { 69 if (args$extra_plot_method == "Counts" || args$extra_plot_method == "Size") {
69 extra_table <- within(extra_table, Counts[Polarity == "R"] <- (Counts[Polarity == "R"] * -1)) 70 extra_table <- within(extra_table, Counts[Polarity == "R"] <- (Counts[Polarity == "R"] * -1)) # nolint
70 } 71 }
71 if (args$extra_plot_method == "Counts" | args$extra_plot_method == "Size" | args$extra_plot_method == "Coverage") { 72 if (args$extra_plot_method == "Counts" || args$extra_plot_method == "Size" || args$extra_plot_method == "Coverage") {
72 i <- 1 73 i <- 1
73 for (sample in samples) { 74 for (sample in samples) {
74 extra_table[, length(extra_table)][extra_table$Dataset == sample] <- extra_table[, length(extra_table)][extra_table$Dataset == sample] * norm_factors[i] 75 extra_table[, length(extra_table)][extra_table$Dataset == sample] <- extra_table[, length(extra_table)][extra_table$Dataset == sample] * norm_factors[i]
75 i <- i + 1 76 i <- i + 1
76 } 77 }
83 if (global == "yes") { 84 if (global == "yes") {
84 bc <- barchart(Counts ~ as.factor(Size) | factor(Dataset, levels = unique(Dataset)), 85 bc <- barchart(Counts ~ as.factor(Size) | factor(Dataset, levels = unique(Dataset)),
85 data = df, origin = 0, 86 data = df, origin = 0,
86 horizontal = FALSE, 87 horizontal = FALSE,
87 col = c("darkblue"), 88 col = c("darkblue"),
88 scales = list(y = list(tick.number = 4, rot = 90, relation = "same", cex = 0.5, alternating = T), x = list(rot = 0, cex = 0.6, tck = 0.5, alternating = c(3, 3))), 89 scales = list(y = list(tick.number = 4, rot = 90, relation = "same", cex = 0.5, alternating = TRUE), x = list(rot = 0, cex = 0.6, tck = 0.5, alternating = c(3, 3))),
89 xlab = list(label = bottom_first_method[[args$first_plot_method]], cex = .85), 90 xlab = list(label = bottom_first_method[[args$first_plot_method]], cex = .85),
90 ylab = list(label = legend_first_method[[args$first_plot_method]], cex = .85), 91 ylab = list(label = legend_first_method[[args$first_plot_method]], cex = .85),
91 main = title_first_method[[args$first_plot_method]], 92 main = title_first_method[[args$first_plot_method]],
92 layout = c(2, 6), newpage = T, 93 layout = c(2, 6), newpage = TRUE,
93 as.table = TRUE, 94 as.table = TRUE,
94 aspect = 0.5, 95 aspect = 0.5,
95 strip = strip.custom(par.strip.text = list(cex = 1), which.given = 1, bg = "lightblue"), 96 strip = strip.custom(par.strip.text = list(cex = 1), which.given = 1, bg = "lightblue"),
96 ... 97 ...
97 ) 98 )
100 data = df, origin = 0, 101 data = df, origin = 0,
101 horizontal = FALSE, 102 horizontal = FALSE,
102 group = Polarity, 103 group = Polarity,
103 stack = TRUE, 104 stack = TRUE,
104 col = c("red", "blue"), 105 col = c("red", "blue"),
105 scales = list(y = list(tick.number = 4, rot = 90, relation = "same", cex = 0.5, alternating = T), x = list(rot = 0, cex = 0.6, tck = 0.5, alternating = c(3, 3))), 106 scales = list(y = list(tick.number = 4, rot = 90, relation = "same", cex = 0.5, alternating = TRUE), x = list(rot = 0, cex = 0.6, tck = 0.5, alternating = c(3, 3))),
106 xlab = list(label = bottom_first_method[[args$first_plot_method]], cex = .85), 107 xlab = list(label = bottom_first_method[[args$first_plot_method]], cex = .85),
107 ylab = list(label = legend_first_method[[args$first_plot_method]], cex = .85), 108 ylab = list(label = legend_first_method[[args$first_plot_method]], cex = .85),
108 main = title_first_method[[args$first_plot_method]], 109 main = title_first_method[[args$first_plot_method]],
109 layout = c(2, 6), newpage = T, 110 layout = c(2, 6), newpage = TRUE,
110 as.table = TRUE, 111 as.table = TRUE,
111 aspect = 0.5, 112 aspect = 0.5,
112 strip = strip.custom(par.strip.text = list(cex = 1), which.given = 1, bg = "lightblue"), 113 strip = strip.custom(par.strip.text = list(cex = 1), which.given = 1, bg = "lightblue"),
113 ... 114 ...
114 ) 115 )
132 data = df, 133 data = df,
133 type = "h", 134 type = "h",
134 lwd = 1.5, 135 lwd = 1.5,
135 scales = list(relation = "free", x = list(rot = 0, cex = 0.7, axs = "i", tck = 0.5), y = list(tick.number = 4, rot = 90, cex = 0.7)), 136 scales = list(relation = "free", x = list(rot = 0, cex = 0.7, axs = "i", tck = 0.5), y = list(tick.number = 4, rot = 90, cex = 0.7)),
136 xlab = NULL, main = NULL, ylab = NULL, ylim = ylimits, 137 xlab = NULL, main = NULL, ylab = NULL, ylim = ylimits,
137 as.table = T, 138 as.table = TRUE,
138 origin = 0, 139 origin = 0,
139 horizontal = FALSE, 140 horizontal = FALSE,
140 group = Polarity, 141 group = Polarity,
141 col = c("red", "blue"), 142 col = c("red", "blue"),
142 par.strip.text = list(cex = 0.7), 143 par.strip.text = list(cex = 0.7),
148 type = ifelse(method == "Coverage", "l", "p"), 149 type = ifelse(method == "Coverage", "l", "p"),
149 pch = 19, 150 pch = 19,
150 cex = 0.35, 151 cex = 0.35,
151 scales = list(relation = "free", x = list(rot = 0, cex = 0.7, axs = "i", tck = 0.5), y = list(tick.number = 4, rot = 90, cex = 0.7)), 152 scales = list(relation = "free", x = list(rot = 0, cex = 0.7, axs = "i", tck = 0.5), y = list(tick.number = 4, rot = 90, cex = 0.7)),
152 xlab = NULL, main = NULL, ylab = NULL, ylim = ylimits, 153 xlab = NULL, main = NULL, ylab = NULL, ylim = ylimits,
153 as.table = T, 154 as.table = TRUE,
154 origin = 0, 155 origin = 0,
155 horizontal = FALSE, 156 horizontal = FALSE,
156 group = Polarity, 157 group = Polarity,
157 col = c("red", "blue"), 158 col = c("red", "blue"),
158 par.strip.text = list(cex = 0.7), 159 par.strip.text = list(cex = 0.7),
250 # main 251 # main
251 252
252 if (args$extra_plot_method != "") { 253 if (args$extra_plot_method != "") {
253 double_plot() 254 double_plot()
254 } 255 }
255 if (args$extra_plot_method == "" & !exists("global", where = args)) { 256 if (args$extra_plot_method == "" && !exists("global", where = args)) {
256 single_plot() 257 single_plot()
257 } 258 }
258 if (exists("global", where = args)) { 259 if (exists("global", where = args)) {
259 pdf(file = args$output, paper = "special", height = 11.69) 260 pdf(file = args$output, paper = "special", height = 11.69)
260 table <- within(table, Counts[Polarity == "R"] <- abs(Counts[Polarity == "R"])) 261 table <- within(table, Counts[Polarity == "R"] <- abs(Counts[Polarity == "R"])) # nolint
261 library(reshape2) 262 library(reshape2)
262 ml <- melt(table, id.vars = c("Dataset", "Chromosome", "Polarity", "Size")) 263 ml <- melt(table, id.vars = c("Dataset", "Chromosome", "Polarity", "Size"))
263 if (args$global == "nomerge") { 264 if (args$global == "nomerge") {
264 castml <- dcast(ml, Dataset + Polarity + Size ~ variable, function(x) sum(x)) 265 castml <- dcast(ml, Dataset + Polarity + Size ~ variable, function(x) sum(x))
265 castml <- within(castml, Counts[Polarity == "R"] <- (Counts[Polarity == "R"] * -1)) 266 castml <- within(castml, Counts[Polarity == "R"] <- (Counts[Polarity == "R"] * -1)) # nolint
266 bc <- globalbc(castml, global = "no") 267 bc <- globalbc(castml, global = "no")
267 } else { 268 } else {
268 castml <- dcast(ml, Dataset + Size ~ variable, function(x) sum(x)) 269 castml <- dcast(ml, Dataset + Size ~ variable, function(x) sum(x))
269 bc <- globalbc(castml, global = "yes") 270 bc <- globalbc(castml, global = "yes")
270 } 271 }