Mercurial > repos > iuc > limma_voom

comparison limma_voom.R @ 21:58c35179ebf0 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/limma_voom commit 127882bd6729d92500ce2a7a51eb5f8949a4c2b5"
author iuc
date Fri, 04 Jun 2021 20:37:04 +0000
parents 0921444c832d
children 708348a17fa1
comparison
equal deleted inserted replaced
20:0921444c832d 21:58c35179ebf0
44 # 44 #
45 # Author: Shian Su - registertonysu@gmail.com - Jan 2014 45 # Author: Shian Su - registertonysu@gmail.com - Jan 2014
46 # Modified by: Maria Doyle - Jun 2017, Jan 2018, May 2018 46 # Modified by: Maria Doyle - Jun 2017, Jan 2018, May 2018
47 47
48 # Record starting time 48 # Record starting time
49 timeStart <- as.character(Sys.time()) 49 time_start <- as.character(Sys.time())
50 50
51 # Load all required libraries 51 # Load all required libraries
52 library(methods, quietly=TRUE, warn.conflicts=FALSE) 52 library(methods, quietly = TRUE, warn.conflicts = FALSE)
53 library(statmod, quietly=TRUE, warn.conflicts=FALSE) 53 library(statmod, quietly = TRUE, warn.conflicts = FALSE)
54 library(splines, quietly=TRUE, warn.conflicts=FALSE) 54 library(splines, quietly = TRUE, warn.conflicts = FALSE)
55 library(edgeR, quietly=TRUE, warn.conflicts=FALSE) 55 library(edgeR, quietly = TRUE, warn.conflicts = FALSE)
56 library(limma, quietly=TRUE, warn.conflicts=FALSE) 56 library(limma, quietly = TRUE, warn.conflicts = FALSE)
57 library(scales, quietly=TRUE, warn.conflicts=FALSE) 57 library(scales, quietly = TRUE, warn.conflicts = FALSE)
58 library(getopt, quietly=TRUE, warn.conflicts=FALSE) 58 library(getopt, quietly = TRUE, warn.conflicts = FALSE)
59 library(gplots, quietly=TRUE, warn.conflicts=FALSE) 59 library(gplots, quietly = TRUE, warn.conflicts = FALSE)
60 60
61 ################################################################################ 61 ################################################################################
62 ### Function Declaration 62 ### Function Declaration
63 ################################################################################ 63 ################################################################################
64 # Function to sanitise contrast equations so there are no whitespaces 64 # Function to sanitise contrast equations so there are no whitespaces
65 # surrounding the arithmetic operators, leading or trailing whitespace 65 # surrounding the arithmetic operators, leading or trailing whitespace
66 sanitiseEquation <- function(equation) { 66 sanitise_equation <- function(equation) {
67 equation <- gsub(" *[+] *", "+", equation) 67 equation <- gsub(" *[+] *", "+", equation)
68 equation <- gsub(" *[-] *", "-", equation) 68 equation <- gsub(" *[-] *", "-", equation)
69 equation <- gsub(" *[/] *", "/", equation) 69 equation <- gsub(" *[/] *", "/", equation)
70 equation <- gsub(" *[*] *", "*", equation) 70 equation <- gsub(" *[*] *", "*", equation)
71 equation <- gsub("^\\s+|\\s+$", "", equation) 71 equation <- gsub("^\\s+|\\s+$", "", equation)
72 return(equation) 72 return(equation)
73 } 73 }
74 74
75 # Function to sanitise group information 75 # Function to sanitise group information
76 sanitiseGroups <- function(string) { 76 sanitise_groups <- function(string) {
77 string <- gsub(" *[,] *", ",", string) 77 string <- gsub(" *[,] *", ",", string)
78 string <- gsub("^\\s+|\\s+$", "", string) 78 string <- gsub("^\\s+|\\s+$", "", string)
79 return(string) 79 return(string)
80 } 80 }
81 81
82 # Function to make contrast contain valid R names 82 # Function to make contrast contain valid R names
83 sanitiseContrast <- function(string) { 83 sanitise_contrast <- function(string) {
84 string <- strsplit(string, split="-") 84 string <- strsplit(string, split = "-")
85 string <- lapply(string, make.names) 85 string <- lapply(string, make.names)
86 string <- lapply(string, paste, collapse="-") 86 string <- lapply(string, paste, collapse = "-")
87 return(string) 87 return(string)
88 } 88 }
89 89
90 # Function to change periods to whitespace in a string 90 # Function to change periods to whitespace in a string
91 unmake.names <- function(string) { 91 unmake_names <- function(string) {
92 string <- gsub(".", " ", string, fixed=TRUE) 92 string <- gsub(".", " ", string, fixed = TRUE)
93 return(string) 93 return(string)
94 } 94 }
95 95
96 # Generate output folder and paths 96 # Generate output folder and paths
97 makeOut <- function(filename) { 97 make_out <- function(filename) {
98 return(paste0(opt$outPath, "/", filename)) 98 return(paste0(opt$outPath, "/", filename))
99 } 99 }
100 100
101 # Generating design information 101 # Generating design information
102 pasteListName <- function(string) { 102 paste_listname <- function(string) {
103 return(paste0("factors$", string)) 103 return(paste0("factors$", string))
104 } 104 }
105 105
106 # Create cata function: default path set, default seperator empty and appending 106 # Create cata function: default path set, default seperator empty and appending
107 # true by default (Ripped straight from the cat function with altered argument 107 # true by default (Ripped straight from the cat function with altered argument
121 } 121 }
122 .Internal(cat(list(...), file, sep, fill, labels, append)) 122 .Internal(cat(list(...), file, sep, fill, labels, append))
123 } 123 }
124 124
125 # Function to write code for html head and title 125 # Function to write code for html head and title
126 HtmlHead <- function(title) { 126 html_head <- function(title) {
127 cata("<head>\n") 127 cata("<head>\n")
128 cata("<title>", title, "</title>\n") 128 cata("<title>", title, "</title>\n")
129 cata("</head>\n") 129 cata("</head>\n")
130 } 130 }
131 131
132 # Function to write code for html links 132 # Function to write code for html links
133 HtmlLink <- function(address, label=address) { 133 html_link <- function(address, label = address) {
134 cata("<a href=\"", address, "\" target=\"_blank\">", label, "</a><br />\n") 134 cata("<a href=\"", address, "\" target=\"_blank\">", label, "</a><br />\n")
135 } 135 }
136 136
137 # Function to write code for html images 137 # Function to write code for html images
138 HtmlImage <- function(source, label=source, height=500, width=500) { 138 html_image <- function(source, label = source, height = 500, width = 500) {
139 cata("<img src=\"", source, "\" alt=\"", label, "\" height=\"", height) 139 cata("<img src=\"", source, "\" alt=\"", label, "\" height=\"", height)
140 cata("\" width=\"", width, "\"/>\n") 140 cata("\" width=\"", width, "\"/>\n")
141 } 141 }
142 142
143 # Function to write code for html list items 143 # Function to write code for html list items
144 ListItem <- function(...) { 144 list_item <- function(...) {
145 cata("<li>", ..., "</li>\n") 145 cata("<li>", ..., "</li>\n")
146 } 146 }
147 147
148 TableItem <- function(...) { 148 table_item <- function(...) {
149 cata("<td>", ..., "</td>\n") 149 cata("<td>", ..., "</td>\n")
150 } 150 }
151 151
152 TableHeadItem <- function(...) { 152 table_head_item <- function(...) {
153 cata("<th>", ..., "</th>\n") 153 cata("<th>", ..., "</th>\n")
154 } 154 }
155 155
156 ################################################################################ 156 ################################################################################
157 ### Input Processing 157 ### Input Processing
158 ################################################################################ 158 ################################################################################
159 159
160 # Collect arguments from command line 160 # Collect arguments from command line
161 args <- commandArgs(trailingOnly=TRUE) 161 args <- commandArgs(trailingOnly = TRUE)
162 162
163 # Get options, using the spec as defined by the enclosed list. 163 # Get options, using the spec as defined by the enclosed list.
164 # Read the options from the default: commandArgs(TRUE). 164 # Read the options from the default: commandArgs(TRUE).
165 spec <- matrix(c( 165 spec <- matrix(c(
166 "htmlPath", "R", 1, "character", 166 "htmlPath", "R", 1, "character",
188 "weightOpt", "w", 0, "logical", 188 "weightOpt", "w", 0, "logical",
189 "topgenes", "G", 1, "integer", 189 "topgenes", "G", 1, "integer",
190 "treatOpt", "T", 0, "logical", 190 "treatOpt", "T", 0, "logical",
191 "plots", "P", 1, "character", 191 "plots", "P", 1, "character",
192 "libinfoOpt", "L", 0, "logical"), 192 "libinfoOpt", "L", 0, "logical"),
193 byrow=TRUE, ncol=4) 193 byrow = TRUE, ncol = 4)
194 opt <- getopt(spec) 194 opt <- getopt(spec)
195 195
196 196
197 if (is.null(opt$matrixPath) & is.null(opt$filesPath)) { 197 if (is.null(opt$matrixPath) & is.null(opt$filesPath)) {
198 cat("A counts matrix (or a set of counts files) is required.\n") 198 cat("A counts matrix (or a set of counts files) is required.\n")
199 q(status=1) 199 q(status = 1)
200 } 200 }
201 201
202 if (is.null(opt$cpmReq)) { 202 if (is.null(opt$cpmReq)) {
203 filtCPM <- FALSE 203 filt_cpm <- FALSE
204 } else { 204 } else {
205 filtCPM <- TRUE 205 filt_cpm <- TRUE
206 } 206 }
207 207
208 if (is.null(opt$cntReq) || is.null(opt$sampleReq)) { 208 if (is.null(opt$cntReq) || is.null(opt$sampleReq)) {
209 filtSmpCount <- FALSE 209 filt_smpcount <- FALSE
210 } else { 210 } else {
211 filtSmpCount <- TRUE 211 filt_smpcount <- TRUE
212 } 212 }
213 213
214 if (is.null(opt$totReq)) { 214 if (is.null(opt$totReq)) {
215 filtTotCount <- FALSE 215 filt_totcount <- FALSE
216 } else { 216 } else {
217 filtTotCount <- TRUE 217 filt_totcount <- TRUE
218 } 218 }
219 219
220 if (is.null(opt$rdaOpt)) { 220 if (is.null(opt$rdaOpt)) {
221 wantRda <- FALSE 221 want_rda <- FALSE
222 } else { 222 } else {
223 wantRda <- TRUE 223 want_rda <- TRUE
224 } 224 }
225 225
226 if (is.null(opt$annoPath)) { 226 if (is.null(opt$annoPath)) {
227 haveAnno <- FALSE 227 have_anno <- FALSE
228 } else { 228 } else {
229 haveAnno <- TRUE 229 have_anno <- TRUE
230 } 230 }
231 231
232 if (is.null(opt$filtCounts)) { 232 if (is.null(opt$filtCounts)) {
233 wantFilt <- FALSE 233 want_filt <- FALSE
234 } else { 234 } else {
235 wantFilt <- TRUE 235 want_filt <- TRUE
236 } 236 }
237 237
238 if (is.null(opt$normCounts)) { 238 if (is.null(opt$normCounts)) {
239 wantNorm <- FALSE 239 want_norm <- FALSE
240 } else { 240 } else {
241 wantNorm <- TRUE 241 want_norm <- TRUE
242 } 242 }
243 243
244 if (is.null(opt$robOpt)) { 244 if (is.null(opt$robOpt)) {
245 wantRobust <- FALSE 245 want_robust <- FALSE
246 } else { 246 } else {
247 wantRobust <- TRUE 247 want_robust <- TRUE
248 } 248 }
249 249
250 if (is.null(opt$weightOpt)) { 250 if (is.null(opt$weightOpt)) {
251 wantWeight <- FALSE 251 want_weight <- FALSE
252 } else { 252 } else {
253 wantWeight <- TRUE 253 want_weight <- TRUE
254 } 254 }
255 255
256 if (is.null(opt$trend)) { 256 if (is.null(opt$trend)) {
257 wantTrend <- FALSE 257 want_trend <- FALSE
258 deMethod <- "limma-voom" 258 de_method <- "limma-voom"
259 } else { 259 } else {
260 wantTrend <- TRUE 260 want_trend <- TRUE
261 deMethod <- "limma-trend" 261 de_method <- "limma-trend"
262 priorCount <- opt$trend 262 prior_count <- opt$trend
263 } 263 }
264 264
265 if (is.null(opt$treatOpt)) { 265 if (is.null(opt$treatOpt)) {
266 wantTreat <- FALSE 266 want_treat <- FALSE
267 } else { 267 } else {
268 wantTreat <- TRUE 268 want_treat <- TRUE
269 } 269 }
270 270
271 if (is.null(opt$libinfoOpt)) { 271 if (is.null(opt$libinfoOpt)) {
272 wantLibinfo <- FALSE 272 want_libinfo <- FALSE
273 } else { 273 } else {
274 wantLibinfo <- TRUE 274 want_libinfo <- TRUE
275 } 275 }
276 276
277 277
278 if (!is.null(opt$filesPath)) { 278 if (!is.null(opt$filesPath)) {
279 # Process the separate count files (adapted from DESeq2 wrapper) 279 # Process the separate count files (adapted from DESeq2 wrapper)
280 library("rjson") 280 library("rjson")
281 parser <- newJSONParser() 281 parser <- newJSONParser()
282 parser$addData(opt$filesPath) 282 parser$addData(opt$filesPath)
283 factorList <- parser$getObject() 283 factor_list <- parser$getObject()
284 factors <- sapply(factorList, function(x) x[[1]]) 284 factors <- sapply(factor_list, function(x) x[[1]])
285 filenamesIn <- unname(unlist(factorList[[1]][[2]])) 285 filenames_in <- unname(unlist(factor_list[[1]][[2]]))
286 sampleTable <- data.frame(sample=basename(filenamesIn), 286 sampletable <- data.frame(sample = basename(filenames_in),
287 filename=filenamesIn, 287 filename = filenames_in,
288 row.names=filenamesIn, 288 row.names = filenames_in,
289 stringsAsFactors=FALSE) 289 stringsAsFactors = FALSE)
290 for (factor in factorList) { 290 for (factor in factor_list) {
291 factorName <- factor[[1]] 291 factorname <- factor[[1]]
292 sampleTable[[factorName]] <- character(nrow(sampleTable)) 292 sampletable[[factorname]] <- character(nrow(sampletable))
293 lvls <- sapply(factor[[2]], function(x) names(x)) 293 lvls <- sapply(factor[[2]], function(x) names(x))
294 for (i in seq_along(factor[[2]])) { 294 for (i in seq_along(factor[[2]])) {
295 files <- factor[[2]][[i]][[1]] 295 files <- factor[[2]][[i]][[1]]
296 sampleTable[files,factorName] <- lvls[i] 296 sampletable[files, factorname] <- lvls[i]
297 } 297 }
298 sampleTable[[factorName]] <- factor(sampleTable[[factorName]], levels=lvls) 298 sampletable[[factorname]] <- factor(sampletable[[factorname]], levels = lvls)
299 } 299 }
300 rownames(sampleTable) <- sampleTable$sample 300 rownames(sampletable) <- sampletable$sample
301 rem <- c("sample","filename") 301 rem <- c("sample", "filename")
302 factors <- sampleTable[, !(names(sampleTable) %in% rem), drop=FALSE] 302 factors <- sampletable[, !(names(sampletable) %in% rem), drop = FALSE]
303 303
304 #read in count files and create single table 304 #read in count files and create single table
305 countfiles <- lapply(sampleTable$filename, function(x){read.delim(x, row.names=1)}) 305 countfiles <- lapply(sampletable$filename, function(x) {
306 read.delim(x, row.names = 1)
307 })
306 counts <- do.call("cbind", countfiles) 308 counts <- do.call("cbind", countfiles)
307 309
308 } else { 310 } else {
309 # Process the single count matrix 311 # Process the single count matrix
310 counts <- read.table(opt$matrixPath, header=TRUE, sep="\t", strip.white=TRUE, stringsAsFactors=FALSE, check.names=FALSE) 312 counts <- read.table(opt$matrixPath, header = TRUE, sep = "\t", strip.white = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
311 row.names(counts) <- counts[, 1] 313 row.names(counts) <- counts[, 1]
312 counts <- counts[ , -1] 314 counts <- counts[, -1]
313 countsRows <- nrow(counts) 315 countsrows <- nrow(counts)
314 316
315 # Process factors 317 # Process factors
316 if (is.null(opt$factInput)) { 318 if (is.null(opt$factInput)) {
317 factorData <- read.table(opt$factFile, header=TRUE, sep="\t", strip.white=TRUE) 319 factordata <- read.table(opt$factFile, header = TRUE, sep = "\t", strip.white = TRUE)
318 if(!setequal(factorData[, 1], colnames(counts))) 320 if (!setequal(factordata[, 1], colnames(counts)))
319 stop("Sample IDs in counts and factors files don't match") 321 stop("Sample IDs in counts and factors files don't match")
320 # order samples as in counts matrix 322 # order samples as in counts matrix
321 factorData <- factorData[match(colnames(counts), factorData[, 1]), ] 323 factordata <- factordata[match(colnames(counts), factordata[, 1]), ]
322 factors <- factorData[, -1, drop=FALSE] 324 factors <- factordata[, -1, drop = FALSE]
323 } else { 325 } else {
324 factors <- unlist(strsplit(opt$factInput, "|", fixed=TRUE)) 326 factors <- unlist(strsplit(opt$factInput, "|", fixed = TRUE))
325 factorData <- list() 327 factordata <- list()
326 for (fact in factors) { 328 for (fact in factors) {
327 newFact <- unlist(strsplit(fact, split="::")) 329 newfact <- unlist(strsplit(fact, split = "::"))
328 factorData <- rbind(factorData, newFact) 330 factordata <- rbind(factordata, newfact)
329 } # Factors have the form: FACT_NAME::LEVEL,LEVEL,LEVEL,LEVEL,... The first factor is the Primary Factor. 331 } # Factors have the form: FACT_NAME::LEVEL,LEVEL,LEVEL,LEVEL,... The first factor is the Primary Factor.
330 332
331 # Set the row names to be the name of the factor and delete first row 333 # Set the row names to be the name of the factor and delete first row
332 row.names(factorData) <- factorData[, 1] 334 row.names(factordata) <- factordata[, 1]
333 factorData <- factorData[, -1] 335 factordata <- factordata[, -1]
334 factorData <- sapply(factorData, sanitiseGroups) 336 factordata <- sapply(factordata, sanitise_groups)
335 factorData <- sapply(factorData, strsplit, split=",") 337 factordata <- sapply(factordata, strsplit, split = ",")
336 # Transform factor data into data frame of R factor objects 338 # Transform factor data into data frame of R factor objects
337 factors <- data.frame(factorData) 339 factors <- data.frame(factordata)
338 } 340 }
339 } 341 }
340 # check there are the same number of samples in counts and factors 342 # check there are the same number of samples in counts and factors
341 if(nrow(factors) != ncol(counts)) { 343 if (nrow(factors) != ncol(counts)) {
342 stop("There are a different number of samples in the counts files and factors") 344 stop("There are a different number of samples in the counts files and factors")
343 } 345 }
344 # make groups valid R names, required for makeContrasts 346 # make groups valid R names, required for makeContrasts
345 factors <- sapply(factors, make.names) 347 factors <- sapply(factors, make.names)
346 factors <- data.frame(factors) 348 factors <- data.frame(factors)
347 349
348 # if annotation file provided 350 # if annotation file provided
349 if (haveAnno) { 351 if (have_anno) {
350 geneanno <- read.table(opt$annoPath, header=TRUE, sep="\t", quote= "", strip.white=TRUE, stringsAsFactors=FALSE) 352 geneanno <- read.table(opt$annoPath, header = TRUE, sep = "\t", quote = "", strip.white = TRUE, stringsAsFactors = FALSE)
351 } 353 }
352 354
353 #Create output directory 355 #Create output directory
354 dir.create(opt$outPath, showWarnings=FALSE) 356 dir.create(opt$outPath, showWarnings = FALSE)
355 357
356 # Process contrasts 358 # Process contrasts
357 if (is.null(opt$contrastInput)) { 359 if (is.null(opt$contrastInput)) {
358 contrastData <- read.table(opt$contrastFile, header=TRUE, sep="\t", quote= "", strip.white=TRUE, stringsAsFactors=FALSE) 360 contrast_data <- read.table(opt$contrastFile, header = TRUE, sep = "\t", quote = "", strip.white = TRUE, stringsAsFactors = FALSE)
359 contrastData <- contrastData[, 1, drop=TRUE] 361 contrast_data <- contrast_data[, 1, drop = TRUE]
360 } else { 362 } else {
361 # Split up contrasts seperated by comma into a vector then sanitise 363 # Split up contrasts seperated by comma into a vector then sanitise
362 contrastData <- unlist(strsplit(opt$contrastInput, split=",")) 364 contrast_data <- unlist(strsplit(opt$contrastInput, split = ","))
363 } 365 }
364 contrastData <- sanitiseEquation(contrastData) 366 contrast_data <- sanitise_equation(contrast_data)
365 contrastData <- gsub(" ", ".", contrastData, fixed=TRUE) 367 contrast_data <- gsub(" ", ".", contrast_data, fixed = TRUE)
366 368
367 # in case input groups start with numbers make the names valid R names, required for makeContrasts 369 # in case input groups start with numbers make the names valid R names, required for makeContrasts
368 cons <- NULL 370 cons <- NULL
369 cons_d <- NULL 371 cons_d <- NULL
370 for (i in contrastData) { 372 for (i in contrast_data) {
371 373
372 # if the contrast is a difference of differences e.g. (A-B)-(X-Y) 374 # if the contrast is a difference of differences e.g. (A-B)-(X-Y)
373 if (grepl("\\)-\\(", i)) { 375 if (grepl("\\)-\\(", i)) {
374 i <- unlist(strsplit(i, split="\\)-\\(")) 376 i <- unlist(strsplit(i, split = "\\)-\\("))
375 i <- gsub("\\(|\\)","", i) 377 i <- gsub("\\(|\\)", "", i)
376 for (j in i) { 378 for (j in i) {
377 j <- sanitiseContrast(j) 379 j <- sanitise_contrast(j)
378 j <- paste0("(", j, ")") 380 j <- paste0("(", j, ")")
379 cons_d <- append(cons_d, unlist(j)) 381 cons_d <- append(cons_d, unlist(j))
380 } 382 }
381 cons_d <- paste(cons_d, collapse = '-') 383 cons_d <- paste(cons_d, collapse = "-")
382 cons <- append(cons, unlist(cons_d)) 384 cons <- append(cons, unlist(cons_d))
383 } else { 385 } else {
384 i <- sanitiseContrast(i) 386 i <- sanitise_contrast(i)
385 cons <- append(cons, unlist(i)) 387 cons <- append(cons, unlist(i))
386 } 388 }
387 } 389 }
388 390
389 plots <- character() 391 plots <- character()
390 if (!is.null(opt$plots)) { 392 if (!is.null(opt$plots)) {
391 plots <- unlist(strsplit(opt$plots, split=",")) 393 plots <- unlist(strsplit(opt$plots, split = ","))
392 } 394 }
393 395
394 denOutPng <- makeOut("densityplots.png") 396 den_png <- make_out("densityplots.png")
395 denOutPdf <- makeOut("densityplots.pdf") 397 den_pdf <- make_out("densityplots.pdf")
396 cpmOutPdf <- makeOut("cpmplots.pdf") 398 cpm_pdf <- make_out("cpmplots.pdf")
397 boxOutPng <- makeOut("boxplots.png") 399 box_png <- make_out("boxplots.png")
398 boxOutPdf <- makeOut("boxplots.pdf") 400 box_pdf <- make_out("boxplots.pdf")
399 mdsscreeOutPng <- makeOut("mdsscree.png") 401 mdsscree_png <- make_out("mdsscree.png")
400 mdsscreeOutPdf <- makeOut("mdsscree.pdf") 402 mdsscree_pdf <- make_out("mdsscree.pdf")
401 mdsxOutPdf <- makeOut("mdsplot_extra.pdf") 403 mdsx_pdf <- make_out("mdsplot_extra.pdf")
402 mdsxOutPng <- makeOut("mdsplot_extra.png") 404 mdsx_png <- make_out("mdsplot_extra.png")
403 mdsamOutPdf <- makeOut("mdplots_samples.pdf") 405 mdsam_pdf <- make_out("mdplots_samples.pdf")
404 mdOutPdf <- character() # Initialise character vector 406 md_pdf <- character() # Initialise character vector
405 volOutPdf <- character() 407 vol_pdf <- character()
406 heatOutPdf <- character() 408 heat_pdf <- character()
407 stripOutPdf <- character() 409 strip_pdf <- character()
408 mdvolOutPng <- character() 410 mdvol_png <- character()
409 topOut <- character() 411 top_out <- character()
410 glimmaOut <- character() 412 glimma_out <- character()
411 for (i in 1:length(cons)) { 413 for (i in seq_along(cons)) {
412 con <- cons[i] 414 con <- cons[i]
413 con <- gsub("\\(|\\)", "", con) 415 con <- gsub("\\(|\\)", "", con)
414 mdOutPdf[i] <- makeOut(paste0("mdplot_", con, ".pdf")) 416 md_pdf[i] <- make_out(paste0("mdplot_", con, ".pdf"))
415 volOutPdf[i] <- makeOut(paste0("volplot_", con, ".pdf")) 417 vol_pdf[i] <- make_out(paste0("volplot_", con, ".pdf"))
416 heatOutPdf[i] <- makeOut(paste0("heatmap_", con, ".pdf")) 418 heat_pdf[i] <- make_out(paste0("heatmap_", con, ".pdf"))
417 stripOutPdf[i] <- makeOut(paste0("stripcharts_", con, ".pdf")) 419 strip_pdf[i] <- make_out(paste0("stripcharts_", con, ".pdf"))
418 mdvolOutPng[i] <- makeOut(paste0("mdvolplot_", con, ".png")) 420 mdvol_png[i] <- make_out(paste0("mdvolplot_", con, ".png"))
419 topOut[i] <- makeOut(paste0(deMethod, "_", con, ".tsv")) 421 top_out[i] <- make_out(paste0(de_method, "_", con, ".tsv"))
420 glimmaOut[i] <- makeOut(paste0("glimma_", con, "/MD-Plot.html")) 422 glimma_out[i] <- make_out(paste0("glimma_", con, "/MD-Plot.html"))
421 } 423 }
422 filtOut <- makeOut(paste0(deMethod, "_", "filtcounts")) 424 filt_out <- make_out(paste0(de_method, "_", "filtcounts"))
423 normOut <- makeOut(paste0(deMethod, "_", "normcounts")) 425 norm_out <- make_out(paste0(de_method, "_", "normcounts"))
424 rdaOut <- makeOut(paste0(deMethod, "_analysis.RData")) 426 rda_out <- make_out(paste0(de_method, "_analysis.RData"))
425 sessionOut <- makeOut("session_info.txt") 427 session_out <- make_out("session_info.txt")
426 428
427 # Initialise data for html links and images, data frame with columns Label and 429 # Initialise data for html links and images, data frame with columns Label and
428 # Link 430 # Link
429 linkData <- data.frame(Label=character(), Link=character(), 431 link_data <- data.frame(Label = character(), Link = character(),
430 stringsAsFactors=FALSE) 432 stringsAsFactors = FALSE)
431 imageData <- data.frame(Label=character(), Link=character(), 433 image_data <- data.frame(Label = character(), Link = character(),
432 stringsAsFactors=FALSE) 434 stringsAsFactors = FALSE)
433 435
434 # Initialise vectors for storage of up/down/neutral regulated counts 436 # Initialise vectors for storage of up/down/neutral regulated counts
435 upCount <- numeric() 437 up_count <- numeric()
436 downCount <- numeric() 438 down_count <- numeric()
437 flatCount <- numeric() 439 flat_count <- numeric()
438 440
439 ################################################################################ 441 ################################################################################
440 ### Data Processing 442 ### Data Processing
441 ################################################################################ 443 ################################################################################
442 444
443 # Extract counts and annotation data 445 # Extract counts and annotation data
444 print("Extracting counts") 446 print("Extracting counts")
445 data <- list() 447 data <- list()
446 data$counts <- counts 448 data$counts <- counts
447 if (haveAnno) { 449 if (have_anno) {
448 # order annotation by genes in counts (assumes gene ids are in 1st column of geneanno) 450 # order annotation by genes in counts (assumes gene ids are in 1st column of geneanno)
449 annoord <- geneanno[match(row.names(counts), geneanno[,1]), ] 451 annoord <- geneanno[match(row.names(counts), geneanno[, 1]), ]
450 data$genes <- annoord 452 data$genes <- annoord
451 } else { 453 } else {
452 data$genes <- data.frame(GeneID=row.names(counts)) 454 data$genes <- data.frame(GeneID = row.names(counts))
453 } 455 }
454 456
455 # Creating naming data 457 # Creating naming data
456 samplenames <- colnames(data$counts) 458 samplenames <- colnames(data$counts)
457 sampleanno <- data.frame("sampleID"=samplenames, factors) 459 sampleanno <- data.frame("sampleID" = samplenames, factors)
458 row.names(factors) <- samplenames # for "Summary of experimental data" table 460 row.names(factors) <- samplenames # for "Summary of experimental data" table
459 461
460 # Creating colours for the groups 462 # Creating colours for the groups
461 cols <- as.numeric(factors[, 1]) 463 cols <- as.numeric(factors[, 1])
462 col.group <- palette()[cols] 464 col.group <- palette()[cols]
463 465
464 # If filter crieteria set, filter out genes that do not have a required cpm/counts in a required number of 466 # If filter crieteria set, filter out genes that do not have a required cpm/counts in a required number of
465 # samples. Default is no filtering 467 # samples. Default is no filtering
466 preFilterCount <- nrow(data$counts) 468 prefilter_count <- nrow(data$counts)
467 nsamples <- ncol(data$counts) 469 nsamples <- ncol(data$counts)
468 470
469 if (filtCPM || filtSmpCount || filtTotCount) { 471 if (filt_cpm || filt_smpcount || filt_totcount) {
470 472
471 if (filtTotCount) { 473 if (filt_totcount) {
472 keep <- rowSums(data$counts) >= opt$cntReq 474 keep <- rowSums(data$counts) >= opt$cntReq
473 } else if (filtSmpCount) { 475 } else if (filt_smpcount) {
474 keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq 476 keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq
475 } else if (filtCPM) { 477 } else if (filt_cpm) {
476 myCPM <- cpm(data$counts) 478
477 thresh <- myCPM >= opt$cpmReq 479 keep <- rowSums(cpm(data$counts) >= opt$cpmReq) >= opt$sampleReq
478 keep <- rowSums(thresh) >= opt$sampleReq
479 480
480 if ("c" %in% plots) { 481 if ("c" %in% plots) {
481 # Plot CPM vs raw counts (to check threshold) 482 # Plot CPM vs raw counts (to check threshold)
482 pdf(cpmOutPdf, width=6.5, height=10) 483 pdf(cpm_pdf, width = 6.5, height = 10)
483 par(mfrow=c(3, 2)) 484 par(mfrow = c(3, 2))
484 for (i in 1:nsamples) { 485 for (i in seq_len(nsamples)) {
485 plot(data$counts[, i], myCPM[, i], xlim=c(0,50), ylim=c(0,3), main=samplenames[i], xlab="Raw counts", ylab="CPM") 486 plot(data$counts[, i], myCPM[, i], xlim = c(0, 50), ylim = c(0, 3), main = samplenames[i], xlab = "Raw counts", ylab = "CPM")
486 abline(v=10, col="red", lty=2, lwd=2) 487 abline(v = 10, col = "red", lty = 2, lwd = 2)
487 abline(h=opt$cpmReq, col=4) 488 abline(h = opt$cpmReq, col = 4)
488 } 489 }
489 linkName <- "CpmPlots.pdf" 490 link_name <- "CpmPlots.pdf"
490 linkAddr <- "cpmplots.pdf" 491 link_addr <- "cpmplots.pdf"
491 linkData <- rbind(linkData, data.frame(Label=linkName, Link=linkAddr, stringsAsFactors=FALSE)) 492 link_data <- rbind(link_data, data.frame(Label = link_name, Link = link_addr, stringsAsFactors = FALSE))
492 invisible(dev.off()) 493 invisible(dev.off())
493 } 494 }
494 } 495 }
495 496
496 data$counts <- data$counts[keep, ] 497 data$counts <- data$counts[keep, ]
497 data$genes <- data$genes[keep, , drop=FALSE] 498 data$genes <- data$genes[keep, , drop = FALSE]
498 499
499 if (wantFilt) { 500 if (want_filt) {
500 print("Outputting filtered counts") 501 print("Outputting filtered counts")
501 filt_counts <- data.frame(data$genes, data$counts, check.names=FALSE) 502 filt_counts <- data.frame(data$genes, data$counts, check.names = FALSE)
502 write.table(filt_counts, file=filtOut, row.names=FALSE, sep="\t", quote=FALSE) 503 write.table(filt_counts, file = filt_out, row.names = FALSE, sep = "\t", quote = FALSE)
503 linkData <- rbind(linkData, data.frame(Label=paste0(deMethod, "_", "filtcounts.tsv"), Link=paste0(deMethod, "_", "filtcounts"), stringsAsFactors=FALSE)) 504 link_data <- rbind(link_data, data.frame(Label = paste0(de_method, "_", "filtcounts.tsv"), Link = paste0(de_method, "_", "filtcounts"), stringsAsFactors = FALSE))
504 } 505 }
505 506
506 # Plot Density 507 # Plot Density
507 if ("d" %in% plots) { 508 if ("d" %in% plots) {
508 # PNG 509 # PNG
509 png(denOutPng, width=1000, height=500) 510 png(den_png, width = 1000, height = 500)
510 par(mfrow=c(1,2), cex.axis=0.8) 511 par(mfrow = c(1, 2), cex.axis = 0.8)
511 512
512 # before filtering 513 # before filtering
513 lcpm1 <- cpm(counts, log=TRUE) 514 lcpm1 <- cpm(counts, log = TRUE)
514 plot(density(lcpm1[, 1]), col=col.group[1], lwd=2, las=2, main="", xlab="") 515 plot(density(lcpm1[, 1]), col = col.group[1], lwd = 2, las = 2, main = "", xlab = "")
515 title(main="Density Plot: Raw counts", xlab="Log-cpm") 516 title(main = "Density Plot: Raw counts", xlab = "Log-cpm")
516 for (i in 2:nsamples){ 517 for (i in 2:nsamples) {
517 den <- density(lcpm1[, i]) 518 den <- density(lcpm1[, i])
518 lines(den$x, den$y, col=col.group[i], lwd=2) 519 lines(den$x, den$y, col = col.group[i], lwd = 2)
519 } 520 }
520 521
521 # after filtering 522 # after filtering
522 lcpm2 <- cpm(data$counts, log=TRUE) 523 lcpm2 <- cpm(data$counts, log = TRUE)
523 plot(density(lcpm2[,1]), col=col.group[1], lwd=2, las=2, main="", xlab="") 524 plot(density(lcpm2[, 1]), col = col.group[1], lwd = 2, las = 2, main = "", xlab = "")
524 title(main="Density Plot: Filtered counts", xlab="Log-cpm") 525 title(main = "Density Plot: Filtered counts", xlab = "Log-cpm")
525 for (i in 2:nsamples){ 526 for (i in 2:nsamples) {
526 den <- density(lcpm2[, i]) 527 den <- density(lcpm2[, i])
527 lines(den$x, den$y, col=col.group[i], lwd=2) 528 lines(den$x, den$y, col = col.group[i], lwd = 2)
528 } 529 }
529 legend("topright", samplenames, text.col=col.group, bty="n") 530 legend("topright", samplenames, text.col = col.group, bty = "n")
530 imgName <- "Densityplots.png" 531 img_name <- "Densityplots.png"
531 imgAddr <- "densityplots.png" 532 img_addr <- "densityplots.png"
532 imageData <- rbind(imageData, data.frame(Label=imgName, Link=imgAddr, stringsAsFactors=FALSE)) 533 image_data <- rbind(image_data, data.frame(Label = img_name, Link = img_addr, stringsAsFactors = FALSE))
533 invisible(dev.off()) 534 invisible(dev.off())
534 535
535 # PDF 536 # PDF
536 pdf(denOutPdf, width=14) 537 pdf(den_pdf, width = 14)
537 par(mfrow=c(1,2), cex.axis=0.8) 538 par(mfrow = c(1, 2), cex.axis = 0.8)
538 plot(density(lcpm1[, 1]), col=col.group[1], lwd=2, las=2, main="", xlab="") 539 plot(density(lcpm1[, 1]), col = col.group[1], lwd = 2, las = 2, main = "", xlab = "")
539 title(main="Density Plot: Raw counts", xlab="Log-cpm") 540 title(main = "Density Plot: Raw counts", xlab = "Log-cpm")
540 for (i in 2:nsamples){ 541 for (i in 2:nsamples) {
541 den <- density(lcpm1[, i]) 542 den <- density(lcpm1[, i])
542 lines(den$x, den$y, col=col.group[i], lwd=2) 543 lines(den$x, den$y, col = col.group[i], lwd = 2)
543 } 544 }
544 plot(density(lcpm2[, 1]), col=col.group[1], lwd=2, las=2, main="", xlab="") 545 plot(density(lcpm2[, 1]), col = col.group[1], lwd = 2, las = 2, main = "", xlab = "")
545 title(main="Density Plot: Filtered counts", xlab="Log-cpm") 546 title(main = "Density Plot: Filtered counts", xlab = "Log-cpm")
546 for (i in 2:nsamples){ 547 for (i in 2:nsamples) {
547 den <- density(lcpm2[, i]) 548 den <- density(lcpm2[, i])
548 lines(den$x, den$y, col=col.group[i], lwd=2) 549 lines(den$x, den$y, col = col.group[i], lwd = 2)
549 } 550 }
550 legend("topright", samplenames, text.col=col.group, bty="n") 551 legend("topright", samplenames, text.col = col.group, bty = "n")
551 linkName <- "DensityPlots.pdf" 552 link_name <- "DensityPlots.pdf"
552 linkAddr <- "densityplots.pdf" 553 link_addr <- "densityplots.pdf"
553 linkData <- rbind(linkData, data.frame(Label=linkName, Link=linkAddr, stringsAsFactors=FALSE)) 554 link_data <- rbind(link_data, data.frame(Label = link_name, Link = link_addr, stringsAsFactors = FALSE))
554 invisible(dev.off()) 555 invisible(dev.off())
555 } 556 }
556 } 557 }
557 558
558 postFilterCount <- nrow(data$counts) 559 postfilter_count <- nrow(data$counts)
559 filteredCount <- preFilterCount-postFilterCount 560 filtered_count <- prefilter_count - postfilter_count
560 561
561 # Generating the DGEList object "y" 562 # Generating the DGEList object "y"
562 print("Generating DGEList object") 563 print("Generating DGEList object")
563 data$samples <- sampleanno 564 data$samples <- sampleanno
564 data$samples$lib.size <- colSums(data$counts) 565 data$samples$lib.size <- colSums(data$counts) # nolint
565 data$samples$norm.factors <- 1 566 data$samples$norm.factors <- 1
566 row.names(data$samples) <- colnames(data$counts) 567 row.names(data$samples) <- colnames(data$counts)
567 y <- new("DGEList", data) 568 y <- new("DGEList", data)
568 569
569 print("Generating Design") 570 print("Generating Design")
570 factorList <- sapply(names(factors), pasteListName) 571 factor_list <- sapply(names(factors), paste_listname)
571 formula <- "~0" 572 formula <- "~0"
572 for (i in 1:length(factorList)) { 573 for (i in seq_along(factor_list)) {
573 formula <- paste(formula,factorList[i], sep="+") 574 formula <- paste(formula, factor_list[i], sep = "+")
574 } 575 }
575 formula <- formula(formula) 576 formula <- formula(formula)
576 design <- model.matrix(formula) 577 design <- model.matrix(formula)
577 for (i in 1:length(factorList)) { 578 for (i in seq_along(factor_list)) {
578 colnames(design) <- gsub(factorList[i], "", colnames(design), fixed=TRUE) 579 colnames(design) <- gsub(factor_list[i], "", colnames(design), fixed = TRUE)
579 } 580 }
580 581
581 # Calculating normalising factors 582 # Calculating normalising factors
582 print("Calculating Normalisation Factors") 583 print("Calculating Normalisation Factors")
583 logcounts <- y #store for plots 584 logcounts <- y #store for plots
584 y <- calcNormFactors(y, method=opt$normOpt) 585 y <- calcNormFactors(y, method = opt$normOpt)
585 586
586 # Generate contrasts information 587 # Generate contrasts information
587 print("Generating Contrasts") 588 print("Generating Contrasts")
588 contrasts <- makeContrasts(contrasts=cons, levels=design) 589 contrasts <- makeContrasts(contrasts = cons, levels = design)
589 590
590 ################################################################################ 591 ################################################################################
591 ### Data Output 592 ### Data Output
592 ################################################################################ 593 ################################################################################
593 594
594 # Plot Box plots (before and after normalisation) 595 # Plot Box plots (before and after normalisation)
595 if (opt$normOpt != "none" & "b" %in% plots) { 596 if (opt$normOpt != "none" & "b" %in% plots) {
596 png(boxOutPng, width=1000, height=500) 597 png(box_png, width = 1000, height = 500)
597 par(mfrow=c(1,2), mar=c(6,4,2,2)+0.1) 598 par(mfrow = c(1, 2), mar = c(6, 4, 2, 2) + 0.1)
598 labels <- colnames(counts) 599 labels <- colnames(counts)
599 600
600 lcpm1 <- cpm(y$counts, log=TRUE) 601 lcpm1 <- cpm(y$counts, log = TRUE)
601 boxplot(lcpm1, las=2, col=col.group, xaxt="n", xlab="") 602 boxplot(lcpm1, las = 2, col = col.group, xaxt = "n", xlab = "")
602 axis(1, at=seq_along(labels), labels = FALSE) 603 axis(1, at = seq_along(labels), labels = FALSE)
603 abline(h=median(lcpm1), col=4) 604 abline(h = median(lcpm1), col = 4)
604 text(x=seq_along(labels), y=par("usr")[3]-1, srt=45, adj=1, labels=labels, xpd=TRUE) 605 text(x = seq_along(labels), y = par("usr")[3] - 1, srt = 45, adj = 1, labels = labels, xpd = TRUE)
605 title(main="Box Plot: Unnormalised counts", ylab="Log-cpm") 606 title(main = "Box Plot: Unnormalised counts", ylab = "Log-cpm")
606 607
607 lcpm2 <- cpm(y, log=TRUE) 608 lcpm2 <- cpm(y, log = TRUE)
608 boxplot(lcpm2, las=2, col=col.group, xaxt="n", xlab="") 609 boxplot(lcpm2, las = 2, col = col.group, xaxt = "n", xlab = "")
609 axis(1, at=seq_along(labels), labels = FALSE) 610 axis(1, at = seq_along(labels), labels = FALSE)
610 text(x=seq_along(labels), y=par("usr")[3]-1, srt=45, adj=1, labels=labels, xpd=TRUE) 611 text(x = seq_along(labels), y = par("usr")[3] - 1, srt = 45, adj = 1, labels = labels, xpd = TRUE)
611 abline(h=median(lcpm2), col=4) 612 abline(h = median(lcpm2), col = 4)
612 title(main="Box Plot: Normalised counts", ylab="Log-cpm") 613 title(main = "Box Plot: Normalised counts", ylab = "Log-cpm")
613 614
614 imgName <- "Boxplots.png" 615 img_name <- "Boxplots.png"
615 imgAddr <- "boxplots.png" 616 img_addr <- "boxplots.png"
616 imageData <- rbind(imageData, data.frame(Label=imgName, Link=imgAddr, stringsAsFactors=FALSE)) 617 image_data <- rbind(image_data, data.frame(Label = img_name, Link = img_addr, stringsAsFactors = FALSE))
617 invisible(dev.off()) 618 invisible(dev.off())
618 619
619 pdf(boxOutPdf, width=14) 620 pdf(box_pdf, width = 14)
620 par(mfrow=c(1,2), mar=c(6,4,2,2)+0.1) 621 par(mfrow = c(1, 2), mar = c(6, 4, 2, 2) + 0.1)
621 boxplot(lcpm1, las=2, col=col.group, xaxt="n", xlab="") 622 boxplot(lcpm1, las = 2, col = col.group, xaxt = "n", xlab = "")
622 axis(1, at=seq_along(labels), labels = FALSE) 623 axis(1, at = seq_along(labels), labels = FALSE)
623 abline(h=median(lcpm1), col=4) 624 abline(h = median(lcpm1), col = 4)
624 text(x=seq_along(labels), y=par("usr")[3]-1, srt=45, adj=1, labels=labels, xpd=TRUE) 625 text(x = seq_along(labels), y = par("usr")[3] - 1, srt = 45, adj = 1, labels = labels, xpd = TRUE)
625 title(main="Box Plot: Unnormalised counts", ylab="Log-cpm") 626 title(main = "Box Plot: Unnormalised counts", ylab = "Log-cpm")
626 boxplot(lcpm2, las=2, col=col.group, xaxt="n", xlab="") 627 boxplot(lcpm2, las = 2, col = col.group, xaxt = "n", xlab = "")
627 axis(1, at=seq_along(labels), labels = FALSE) 628 axis(1, at = seq_along(labels), labels = FALSE)
628 text(x=seq_along(labels), y=par("usr")[3]-1, srt=45, adj=1, labels=labels, xpd=TRUE) 629 text(x = seq_along(labels), y = par("usr")[3] - 1, srt = 45, adj = 1, labels = labels, xpd = TRUE)
629 abline(h=median(lcpm2), col=4) 630 abline(h = median(lcpm2), col = 4)
630 title(main="Box Plot: Normalised counts", ylab="Log-cpm") 631 title(main = "Box Plot: Normalised counts", ylab = "Log-cpm")
631 linkName <- "BoxPlots.pdf" 632 link_name <- "BoxPlots.pdf"
632 linkAddr <- "boxplots.pdf" 633 link_addr <- "boxplots.pdf"
633 linkData <- rbind(linkData, data.frame(Label=linkName, Link=linkAddr, stringsAsFactors=FALSE)) 634 link_data <- rbind(link_data, data.frame(Label = link_name, Link = link_addr, stringsAsFactors = FALSE))
634 invisible(dev.off()) 635 invisible(dev.off())
635 } 636 }
636 637
637 # Plot MDS 638 # Plot MDS
638 print("Generating MDS plot") 639 print("Generating MDS plot")
639 labels <- names(counts) 640 labels <- names(counts)
640 641
641 # Scree plot (Variance Explained) code copied from Glimma 642 # Scree plot (Variance Explained) code copied from Glimma
642 643
643 # get column of matrix 644 # get column of matrix
644 getCols <- function(x, inds) { 645 get_cols <- function(x, inds) {
645 x[, inds, drop=FALSE] 646 x[, inds, drop = FALSE]
646 } 647 }
647 648
648 x <- cpm(y, log=TRUE) 649 x <- cpm(y, log = TRUE)
649 ndim <- nsamples - 1 650 ndim <- nsamples - 1
650 nprobes <- nrow(x) 651 nprobes <- nrow(x)
651 top <- 500 652 top <- 500
652 top <- min(top, nprobes) 653 top <- min(top, nprobes)
653 cn <- colnames(x) 654 cn <- colnames(x)
654 bad <- rowSums(is.finite(x)) < nsamples 655 bad <- rowSums(is.finite(x)) < nsamples
655 656
656 if (any(bad)) { 657 if (any(bad)) {
657 warning("Rows containing infinite values have been removed") 658 warning("Rows containing infinite values have been removed")
658 x <- x[!bad, , drop=FALSE] 659 x <- x[!bad, , drop = FALSE]
659 } 660 }
660 661
661 dd <- matrix(0, nrow=nsamples, ncol=nsamples, dimnames=list(cn, cn)) 662 dd <- matrix(0, nrow = nsamples, ncol = nsamples, dimnames = list(cn, cn))
662 topindex <- nprobes - top + 1L 663 topindex <- nprobes - top + 1L
663 for (i in 2L:(nsamples)) { 664 for (i in 2L:(nsamples)) {
664 for (j in 1L:(i - 1L)) { 665 for (j in 1L:(i - 1L)) {
665 dists <- (getCols(x, i) - getCols(x, j))^2 666 dists <- (get_cols(x, i) - get_cols(x, j))^2
666 dists <- sort.int(dists, partial = topindex ) 667 dists <- sort.int(dists, partial = topindex)
667 topdist <- dists[topindex:nprobes] 668 topdist <- dists[topindex:nprobes]
668 dd[i, j] <- sqrt(mean(topdist)) 669 dd[i, j] <- sqrt(mean(topdist))
669 } 670 }
670 } 671 }
671 672
672 a1 <- suppressWarnings(cmdscale(as.dist(dd), k=min(ndim, 8), eig=TRUE)) 673 a1 <- suppressWarnings(cmdscale(as.dist(dd), k = min(ndim, 8), eig = TRUE))
673 eigen <- data.frame(name = 1:min(ndim, 8), eigen = round(a1$eig[1:min(ndim, 8)]/sum(a1$eig), 2)) 674 eigen <- data.frame(name = 1:min(ndim, 8), eigen = round(a1$eig[1:min(ndim, 8)] / sum(a1$eig), 2))
674 675
675 png(mdsscreeOutPng, width=1000, height=500) 676 png(mdsscree_png, width = 1000, height = 500)
676 par(mfrow=c(1, 2)) 677 par(mfrow = c(1, 2))
677 plotMDS(y, labels=samplenames, col=as.numeric(factors[, 1]), main="MDS Plot: Dims 1 and 2") 678 plotMDS(y, labels = samplenames, col = as.numeric(factors[, 1]), main = "MDS Plot: Dims 1 and 2")
678 barplot(eigen$eigen, names.arg=eigen$name, main = "Scree Plot: Variance Explained", xlab = "Dimension", ylab = "Proportion", las=1) 679 barplot(eigen$eigen, names.arg = eigen$name, main = "Scree Plot: Variance Explained", xlab = "Dimension", ylab = "Proportion", las = 1)
679 imgName <- paste0("MDSPlot_", names(factors)[1], ".png") 680 img_name <- paste0("MDSPlot_", names(factors)[1], ".png")
680 imgAddr <- "mdsscree.png" 681 img_addr <- "mdsscree.png"
681 imageData <- rbind(imageData, data.frame(Label=imgName, Link=imgAddr, stringsAsFactors=FALSE)) 682 image_data <- rbind(image_data, data.frame(Label = img_name, Link = img_addr, stringsAsFactors = FALSE))
682 invisible(dev.off()) 683 invisible(dev.off())
683 684
684 pdf(mdsscreeOutPdf, width=14) 685 pdf(mdsscree_pdf, width = 14)
685 par(mfrow=c(1, 2)) 686 par(mfrow = c(1, 2))
686 plotMDS(y, labels=samplenames, col=as.numeric(factors[, 1]), main="MDS Plot: Dims 1 and 2") 687 plotMDS(y, labels = samplenames, col = as.numeric(factors[, 1]), main = "MDS Plot: Dims 1 and 2")
687 barplot(eigen$eigen, names.arg=eigen$name, main = "Scree Plot: Variance Explained", xlab = "Dimension", ylab = "Proportion", las=1) 688 barplot(eigen$eigen, names.arg = eigen$name, main = "Scree Plot: Variance Explained", xlab = "Dimension", ylab = "Proportion", las = 1)
688 linkName <- paste0("MDSPlot_", names(factors)[1], ".pdf") 689 link_name <- paste0("MDSPlot_", names(factors)[1], ".pdf")
689 linkAddr <- "mdsscree.pdf" 690 link_addr <- "mdsscree.pdf"
690 linkData <- rbind(linkData, data.frame(Label=linkName, Link=linkAddr, stringsAsFactors=FALSE)) 691 link_data <- rbind(link_data, data.frame(Label = link_name, Link = link_addr, stringsAsFactors = FALSE))
691 invisible(dev.off()) 692 invisible(dev.off())
692 693
693 # generate Glimma interactive MDS Plot 694 # generate Glimma interactive MDS Plot
694 if ("i" %in% plots) { 695 if ("i" %in% plots) {
695 Glimma::glMDSPlot(y, labels=samplenames, groups=factors[, 1], 696 Glimma::glMDSPlot(y, labels = samplenames, groups = factors[, 1],
696 folder="glimma_MDS", launch=FALSE) 697 folder = "glimma_MDS", launch = FALSE)
697 linkName <- "Glimma_MDSPlot.html" 698 link_name <- "Glimma_MDSPlot.html"
698 linkAddr <- "glimma_MDS/MDS-Plot.html" 699 link_addr <- "glimma_MDS/MDS-Plot.html"
699 linkData <- rbind(linkData, c(linkName, linkAddr)) 700 link_data <- rbind(link_data, c(link_name, link_addr))
700 } 701 }
701 702
702 if ("x" %in% plots) { 703 if ("x" %in% plots) {
703 png(mdsxOutPng, width=1000, height=500) 704 png(mdsx_png, width = 1000, height = 500)
704 par(mfrow=c(1, 2)) 705 par(mfrow = c(1, 2))
705 for (i in 2:3) { 706 for (i in 2:3) {
706 dim1 <- i 707 dim1 <- i
707 dim2 <- i + 1 708 dim2 <- i + 1
708 plotMDS(y, dim=c(dim1, dim2), labels=samplenames, col=as.numeric(factors[, 1]), main=paste("MDS Plot: Dims", dim1, "and", dim2)) 709 plotMDS(y, dim = c(dim1, dim2), labels = samplenames, col = as.numeric(factors[, 1]), main = paste("MDS Plot: Dims", dim1, "and", dim2))
709 } 710 }
710 imgName <- paste0("MDSPlot_extra.png") 711 img_name <- paste0("MDSPlot_extra.png")
711 imgAddr <- paste0("mdsplot_extra.png") 712 img_addr <- paste0("mdsplot_extra.png")
712 imageData <- rbind(imageData, data.frame(Label=imgName, Link=imgAddr, stringsAsFactors=FALSE)) 713 image_data <- rbind(image_data, data.frame(Label = img_name, Link = img_addr, stringsAsFactors = FALSE))
713 invisible(dev.off()) 714 invisible(dev.off())
714 715
715 pdf(mdsxOutPdf, width=14) 716 pdf(mdsx_pdf, width = 14)
716 par(mfrow=c(1, 2)) 717 par(mfrow = c(1, 2))
717 for (i in 2:3) { 718 for (i in 2:3) {
718 dim1 <- i 719 dim1 <- i
719 dim2 <- i + 1 720 dim2 <- i + 1
720 plotMDS(y, dim=c(dim1, dim2), labels=samplenames, col=as.numeric(factors[, 1]), main=paste("MDS Plot: Dims", dim1, "and", dim2)) 721 plotMDS(y, dim = c(dim1, dim2), labels = samplenames, col = as.numeric(factors[, 1]), main = paste("MDS Plot: Dims", dim1, "and", dim2))
721 } 722 }
722 linkName <- "MDSPlot_extra.pdf" 723 link_name <- "MDSPlot_extra.pdf"
723 linkAddr <- "mdsplot_extra.pdf" 724 link_addr <- "mdsplot_extra.pdf"
724 linkData <- rbind(linkData, data.frame(Label=linkName, Link=linkAddr, stringsAsFactors=FALSE)) 725 link_data <- rbind(link_data, data.frame(Label = link_name, Link = link_addr, stringsAsFactors = FALSE))
725 invisible(dev.off()) 726 invisible(dev.off())
726 } 727 }
727 728
728 if ("m" %in% plots) { 729 if ("m" %in% plots) {
729 # Plot MD plots for individual samples 730 # Plot MD plots for individual samples
730 print("Generating MD plots for samples") 731 print("Generating MD plots for samples")
731 pdf(mdsamOutPdf, width=6.5, height=10) 732 pdf(mdsam_pdf, width = 6.5, height = 10)
732 par(mfrow=c(3, 2)) 733 par(mfrow = c(3, 2))
733 for (i in 1:nsamples) { 734 for (i in 1:nsamples) {
734 if (opt$normOpt != "none") { 735 if (opt$normOpt != "none") {
735 plotMD(logcounts, column=i, main=paste(colnames(logcounts)[i], "(before)")) 736 plotMD(logcounts, column = i, main = paste(colnames(logcounts)[i], "(before)"))
736 abline(h=0, col="red", lty=2, lwd=2) 737 abline(h = 0, col = "red", lty = 2, lwd = 2)
737 } 738 }
738 plotMD(y, column=i) 739 plotMD(y, column = i)
739 abline(h=0, col="red", lty=2, lwd=2) 740 abline(h = 0, col = "red", lty = 2, lwd = 2)
740 } 741 }
741 linkName <- "MDPlots_Samples.pdf" 742 link_name <- "MDPlots_Samples.pdf"
742 linkAddr <- "mdplots_samples.pdf" 743 link_addr <- "mdplots_samples.pdf"
743 linkData <- rbind(linkData, c(linkName, linkAddr)) 744 link_data <- rbind(link_data, c(link_name, link_addr))
744 invisible(dev.off()) 745 invisible(dev.off())
745 } 746 }
746 747
747 748
748 if (wantTrend) { 749 if (want_trend) {
749 # limma-trend approach 750 # limma-trend approach
750 logCPM <- cpm(y, log=TRUE, prior.count=opt$trend) 751 logcpm <- cpm(y, log = TRUE, prior.count = opt$trend)
751 fit <- lmFit(logCPM, design) 752 fit <- lmFit(logcpm, design)
752 fit$genes <- y$genes 753 fit$genes <- y$genes
753 fit <- contrasts.fit(fit, contrasts) 754 fit <- contrasts.fit(fit, contrasts)
754 if (wantRobust) { 755 if (want_robust) {
755 fit <- eBayes(fit, trend=TRUE, robust=TRUE) 756 fit <- eBayes(fit, trend = TRUE, robust = TRUE)
756 } else { 757 } else {
757 fit <- eBayes(fit, trend=TRUE, robust=FALSE) 758 fit <- eBayes(fit, trend = TRUE, robust = FALSE)
758 } 759 }
759 760
760 plotData <- logCPM 761 plot_data <- logcpm
761 762
762 # Save normalised counts (log2cpm) 763 # Save normalised counts (log2cpm)
763 if (wantNorm) { 764 if (want_norm) {
764 write.table(logCPM, file=normOut, row.names=TRUE, sep="\t", quote=FALSE) 765 write.table(logcpm, file = norm_out, row.names = TRUE, sep = "\t", quote = FALSE)
765 linkData <- rbind(linkData, c((paste0(deMethod, "_", "normcounts.tsv")), (paste0(deMethod, "_", "normcounts")))) 766 link_data <- rbind(link_data, c((paste0(de_method, "_", "normcounts.tsv")), (paste0(de_method, "_", "normcounts"))))
766 } 767 }
767 } else { 768 } else {
768 # limma-voom approach 769 # limma-voom approach
769 770
770 if (wantWeight) { 771 if (want_weight) {
771 voomWtsOutPdf <- makeOut("voomwtsplot.pdf") 772 voomwts_pdf <- make_out("voomwtsplot.pdf")
772 voomWtsOutPng <- makeOut("voomwtsplot.png") 773 voomwts_png <- make_out("voomwtsplot.png")
773 # Creating voom data object and plot 774 # Creating voom data object and plot
774 png(voomWtsOutPng, width=1000, height=500) 775 png(voomwts_png, width = 1000, height = 500)
775 vData <- voomWithQualityWeights(y, design=design, plot=TRUE) 776 vdata <- voomWithQualityWeights(y, design = design, plot = TRUE)
776 imgName <- "VoomWithQualityWeightsPlot.png" 777 img_name <- "VoomWithQualityWeightsPlot.png"
777 imgAddr <- "voomwtsplot.png" 778 img_addr <- "voomwtsplot.png"
778 imageData <- rbind(imageData, c(imgName, imgAddr)) 779 image_data <- rbind(image_data, c(img_name, img_addr))
779 invisible(dev.off()) 780 invisible(dev.off())
780 781
781 pdf(voomWtsOutPdf, width=14) 782 pdf(voomwts_pdf, width = 14)
782 vData <- voomWithQualityWeights(y, design=design, plot=TRUE) 783 vdata <- voomWithQualityWeights(y, design = design, plot = TRUE)
783 linkName <- "VoomWithQualityWeightsPlot.pdf" 784 link_name <- "VoomWithQualityWeightsPlot.pdf"
784 linkAddr <- "voomwtsplot.pdf" 785 link_addr <- "voomwtsplot.pdf"
785 linkData <- rbind(linkData, c(linkName, linkAddr)) 786 link_data <- rbind(link_data, c(link_name, link_addr))
786 invisible(dev.off()) 787 invisible(dev.off())
787 788
788 # Generating fit data and top table with weights 789 # Generating fit data and top table with weights
789 wts <- vData$weights 790 wts <- vdata$weights
790 voomFit <- lmFit(vData, design, weights=wts) 791 voomfit <- lmFit(vdata, design, weights = wts)
791 792
792 } else { 793 } else {
793 voomOutPdf <- makeOut("voomplot.pdf") 794 voom_pdf <- make_out("voomplot.pdf")
794 voomOutPng <- makeOut("voomplot.png") 795 voom_png <- make_out("voomplot.png")
795 # Creating voom data object and plot 796 # Creating voom data object and plot
796 png(voomOutPng, width=500, height=500) 797 png(voom_png, width = 500, height = 500)
797 vData <- voom(y, design=design, plot=TRUE) 798 vdata <- voom(y, design = design, plot = TRUE)
798 imgName <- "VoomPlot" 799 img_name <- "VoomPlot"
799 imgAddr <- "voomplot.png" 800 img_addr <- "voomplot.png"
800 imageData <- rbind(imageData, c(imgName, imgAddr)) 801 image_data <- rbind(image_data, c(img_name, img_addr))
801 invisible(dev.off()) 802 invisible(dev.off())
802 803
803 pdf(voomOutPdf) 804 pdf(voom_pdf)
804 vData <- voom(y, design=design, plot=TRUE) 805 vdata <- voom(y, design = design, plot = TRUE)
805 linkName <- "VoomPlot.pdf" 806 link_name <- "VoomPlot.pdf"
806 linkAddr <- "voomplot.pdf" 807 link_addr <- "voomplot.pdf"
807 linkData <- rbind(linkData, c(linkName, linkAddr)) 808 link_data <- rbind(link_data, c(link_name, link_addr))
808 invisible(dev.off()) 809 invisible(dev.off())
809 810
810 # Generate voom fit 811 # Generate voom fit
811 voomFit <- lmFit(vData, design) 812 voomfit <- lmFit(vdata, design)
812 } 813 }
813 814
814 # Save normalised counts (log2cpm) 815 # Save normalised counts (log2cpm)
815 if (wantNorm) { 816 if (want_norm) {
816 norm_counts <- data.frame(vData$genes, vData$E, check.names=FALSE) 817 norm_counts <- data.frame(vdata$genes, vdata$E, check.names = FALSE)
817 write.table(norm_counts, file=normOut, row.names=FALSE, sep="\t", quote=FALSE) 818 write.table(norm_counts, file = norm_out, row.names = FALSE, sep = "\t", quote = FALSE)
818 linkData <- rbind(linkData, c((paste0(deMethod, "_", "normcounts.tsv")), (paste0(deMethod, "_", "normcounts")))) 819 link_data <- rbind(link_data, c((paste0(de_method, "_", "normcounts.tsv")), (paste0(de_method, "_", "normcounts"))))
819 } 820 }
820 821
821 # Fit linear model and estimate dispersion with eBayes 822 # Fit linear model and estimate dispersion with eBayes
822 voomFit <- contrasts.fit(voomFit, contrasts) 823 voomfit <- contrasts.fit(voomfit, contrasts)
823 if (wantRobust) { 824 if (want_robust) {
824 fit <- eBayes(voomFit, robust=TRUE) 825 fit <- eBayes(voomfit, robust = TRUE)
825 } else { 826 } else {
826 fit <- eBayes(voomFit, robust=FALSE) 827 fit <- eBayes(voomfit, robust = FALSE)
827 } 828 }
828 plotData <- vData 829 plot_data <- vdata
829 } 830 }
830 831
831 # plot final model mean-variance trend with plotSA 832 # plot final model mean-variance trend with plotSA
832 saOutPng <- makeOut("saplot.png") 833 sa_png <- make_out("saplot.png")
833 saOutPdf <- makeOut("saplot.pdf") 834 sa_pdf <- make_out("saplot.pdf")
834 835
835 png(saOutPng, width=500, height=500) 836 png(sa_png, width = 500, height = 500)
836 plotSA(fit, main="Final model: Mean-variance trend (SA Plot)") 837 plotSA(fit, main = "Final model: Mean-variance trend (SA Plot)")
837 imgName <- "SAPlot.png" 838 img_name <- "SAPlot.png"
838 imgAddr <- "saplot.png" 839 img_addr <- "saplot.png"
839 imageData <- rbind(imageData, c(imgName, imgAddr)) 840 image_data <- rbind(image_data, c(img_name, img_addr))
840 invisible(dev.off()) 841 invisible(dev.off())
841 842
842 pdf(saOutPdf) 843 pdf(sa_pdf)
843 plotSA(fit, main="Final model: Mean-variance trend (SA Plot)") 844 plotSA(fit, main = "Final model: Mean-variance trend (SA Plot)")
844 linkName <- "SAPlot.pdf" 845 link_name <- "SAPlot.pdf"
845 linkAddr <- "saplot.pdf" 846 link_addr <- "saplot.pdf"
846 linkData <- rbind(linkData, c(linkName, linkAddr)) 847 link_data <- rbind(link_data, c(link_name, link_addr))
847 invisible(dev.off()) 848 invisible(dev.off())
848 849
849 # Save library size info 850 # Save library size info
850 if (wantLibinfo) { 851 if (want_libinfo) {
851 efflibsize <- round(y$samples$lib.size * y$samples$norm.factors) 852 efflibsize <- round(y$samples$lib.size * y$samples$norm.factors)
852 libsizeinfo <- cbind(y$samples, EffectiveLibrarySize=efflibsize) 853 libsizeinfo <- cbind(y$samples, EffectiveLibrarySize = efflibsize)
853 libsizeinfo$lib.size <- round(libsizeinfo$lib.size) 854 libsizeinfo$lib.size <- round(libsizeinfo$lib.size) # nolint
854 names(libsizeinfo)[names(libsizeinfo)=="sampleID"] <- "SampleID" 855 names(libsizeinfo)[names(libsizeinfo) == "sampleID"] <- "SampleID"
855 names(libsizeinfo)[names(libsizeinfo)=="lib.size"] <- "LibrarySize" 856 names(libsizeinfo)[names(libsizeinfo) == "lib.size"] <- "LibrarySize"
856 names(libsizeinfo)[names(libsizeinfo)=="norm.factors"] <- "NormalisationFactor" 857 names(libsizeinfo)[names(libsizeinfo) == "norm.factors"] <- "NormalisationFactor"
857 write.table(libsizeinfo, file="libsizeinfo", row.names=FALSE, sep="\t", quote=FALSE) 858 write.table(libsizeinfo, file = "libsizeinfo", row.names = FALSE, sep = "\t", quote = FALSE)
858 } 859 }
859 860
860 print("Generating DE results") 861 print("Generating DE results")
861 862
862 if (wantTreat) { 863 if (want_treat) {
863 print("Applying TREAT method") 864 print("Applying TREAT method")
864 if (wantRobust) { 865 if (want_robust) {
865 fit <- treat(fit, lfc=opt$lfcReq, robust=TRUE) 866 fit <- treat(fit, lfc = opt$lfcReq, robust = TRUE)
866 } else { 867 } else {
867 fit <- treat(fit, lfc=opt$lfcReq, robust=FALSE) 868 fit <- treat(fit, lfc = opt$lfcReq, robust = FALSE)
868 } 869 }
869 } 870 }
870 871
871 status = decideTests(fit, adjust.method=opt$pAdjOpt, p.value=opt$pValReq, 872 status <- decideTests(fit, adjust.method = opt$pAdjOpt, p.value = opt$pValReq,
872 lfc=opt$lfcReq) 873 lfc = opt$lfcReq)
873 sumStatus <- summary(status) 874 sum_status <- summary(status)
874 875
875 for (i in 1:length(cons)) { 876 for (i in seq_along(cons)) {
876 con_name <- cons[i] 877 con_name <- cons[i]
877 con <- cons[i] 878 con <- cons[i]
878 con <- gsub("\\(|\\)", "", con) 879 con <- gsub("\\(|\\)", "", con)
879 # Collect counts for differential expression 880 # Collect counts for differential expression
880 upCount[i] <- sumStatus["Up", i] 881 up_count[i] <- sum_status["Up", i]
881 downCount[i] <- sumStatus["Down", i] 882 down_count[i] <- sum_status["Down", i]
882 flatCount[i] <- sumStatus["NotSig", i] 883 flat_count[i] <- sum_status["NotSig", i]
883 884
884 # Write top expressions table 885 # Write top expressions table
885 if (wantTreat) { 886 if (want_treat) {
886 top <- topTreat(fit, coef=i, adjust.method=opt$pAdjOpt, number=Inf, sort.by="P") 887 top <- topTreat(fit, coef = i, adjust.method = opt$pAdjOpt, number = Inf, sort.by = "P")
887 } else{ 888 } else{
888 top <- topTable(fit, coef=i, adjust.method=opt$pAdjOpt, number=Inf, sort.by="P") 889 top <- topTable(fit, coef = i, adjust.method = opt$pAdjOpt, number = Inf, sort.by = "P")
889 } 890 }
890 write.table(top, file=topOut[i], row.names=FALSE, sep="\t", quote=FALSE) 891 write.table(top, file = top_out[i], row.names = FALSE, sep = "\t", quote = FALSE)
891 linkName <- paste0(deMethod, "_", con, ".tsv") 892 link_name <- paste0(de_method, "_", con, ".tsv")
892 linkAddr <- paste0(deMethod, "_", con, ".tsv") 893 link_addr <- paste0(de_method, "_", con, ".tsv")
893 linkData <- rbind(linkData, c(linkName, linkAddr)) 894 link_data <- rbind(link_data, c(link_name, link_addr))
894 895
895 # Plot MD (log ratios vs mean average) using limma package on weighted 896 # Plot MD (log ratios vs mean average) using limma package on weighted
896 pdf(mdOutPdf[i]) 897 pdf(md_pdf[i])
897 limma::plotMD(fit, status=status[, i], coef=i, 898 limma::plotMD(fit, status = status[, i], coef = i,
898 main=paste("MD Plot:", unmake.names(con)), 899 main = paste("MD Plot:", unmake_names(con)),
899 hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1), 900 hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1),
900 xlab="Average Expression", ylab="logFC") 901 xlab = "Average Expression", ylab = "logFC")
901 abline(h=0, col="grey", lty=2) 902 abline(h = 0, col = "grey", lty = 2)
902 linkName <- paste0("MDPlot_", con, ".pdf") 903 link_name <- paste0("MDPlot_", con, ".pdf")
903 linkAddr <- paste0("mdplot_", con, ".pdf") 904 link_addr <- paste0("mdplot_", con, ".pdf")
904 linkData <- rbind(linkData, c(linkName, linkAddr)) 905 link_data <- rbind(link_data, c(link_name, link_addr))
905 invisible(dev.off()) 906 invisible(dev.off())
906 907
907 # Generate Glimma interactive Volcano, MD plot and tables, requires annotation file (assumes gene labels/symbols in 2nd column) 908 # Generate Glimma interactive Volcano, MD plot and tables, requires annotation file (assumes gene labels/symbols in 2nd column)
908 if ("i" %in% plots & haveAnno) { 909 if ("i" %in% plots & have_anno) {
909 # make gene labels unique to handle NAs 910 # make gene labels unique to handle NAs
910 geneanno <- y$genes 911 geneanno <- y$genes
911 geneanno[, 2] <- make.unique(geneanno[, 2]) 912 geneanno[, 2] <- make.unique(geneanno[, 2])
912 913
913 # use the logCPMS for the counts 914 # use the logcpms for the counts
914 if (wantTrend) { 915 if (want_trend) {
915 cnts <- logCPM 916 cnts <- logcpm
916 } else{ 917 } else{
917 cnts <- vData$E 918 cnts <- vdata$E
918 } 919 }
919 920
920 # MD plot 921 # MD plot
921 Glimma::glMDPlot(fit, coef=i, counts=cnts, anno=geneanno, groups=factors[, 1], 922 Glimma::glMDPlot(fit, coef = i, counts = cnts, anno = geneanno, groups = factors[, 1],
922 status=status[, i], sample.cols=col.group, 923 status = status[, i], sample.cols = col.group,
923 main=paste("MD Plot:", unmake.names(con)), side.main=colnames(y$genes)[2], 924 main = paste("MD Plot:", unmake_names(con)), side.main = colnames(y$genes)[2],
924 folder=paste0("glimma_", unmake.names(con)), launch=FALSE) 925 folder = paste0("glimma_", unmake_names(con)), launch = FALSE)
925 linkName <- paste0("Glimma_MDPlot_", con, ".html") 926 link_name <- paste0("Glimma_MDPlot_", con, ".html")
926 linkAddr <- paste0("glimma_", con, "/MD-Plot.html") 927 link_addr <- paste0("glimma_", con, "/MD-Plot.html")
927 linkData <- rbind(linkData, c(linkName, linkAddr)) 928 link_data <- rbind(link_data, c(link_name, link_addr))
928 929
929 # Volcano plot 930 # Volcano plot
930 Glimma::glXYPlot(x=fit$coefficients[, i], y=-log10(fit$p.value[, i]), counts=cnts, anno=geneanno, groups=factors[, 1], 931 Glimma::glXYPlot(x = fit$coefficients[, i], y = -log10(fit$p.value[, i]), counts = cnts, anno = geneanno, groups = factors[, 1],
931 status=status[, i], sample.cols=col.group, 932 status = status[, i], sample.cols = col.group,
932 main=paste("Volcano Plot:", unmake.names(con)), side.main=colnames(y$genes)[2], 933 main = paste("Volcano Plot:", unmake_names(con)), side.main = colnames(y$genes)[2],
933 xlab="logFC", ylab="-log10(P-value)", 934 xlab = "logFC", ylab = "-log10(P-value)",
934 folder=paste0("glimma_volcano_", unmake.names(con)), launch=FALSE) 935 folder = paste0("glimma_volcano_", unmake_names(con)), launch = FALSE)
935 linkName <- paste0("Glimma_VolcanoPlot_", con, ".html") 936 link_name <- paste0("Glimma_VolcanoPlot_", con, ".html")
936 linkAddr <- paste0("glimma_volcano_", con, "/XY-Plot.html") 937 link_addr <- paste0("glimma_volcano_", con, "/XY-Plot.html")
937 linkData <- rbind(linkData, c(linkName, linkAddr)) 938 link_data <- rbind(link_data, c(link_name, link_addr))
938 } 939 }
939 940
940 # Plot Volcano 941 # Plot Volcano
941 pdf(volOutPdf[i]) 942 pdf(vol_pdf[i])
942 if (haveAnno) { 943 if (have_anno) {
943 # labels must be in second column currently 944 # labels must be in second column currently
944 labels <- fit$genes[, 2] 945 labels <- fit$genes[, 2]
945 } else { 946 } else {
946 labels <- fit$genes$GeneID 947 labels <- fit$genes$GeneID
947 } 948 }
948 limma::volcanoplot(fit, coef=i, 949 limma::volcanoplot(fit, coef = i,
949 main=paste("Volcano Plot:", unmake.names(con)), 950 main = paste("Volcano Plot:", unmake_names(con)),
950 highlight=opt$topgenes, 951 highlight = opt$topgenes,
951 names=labels) 952 names = labels)
952 linkName <- paste0("VolcanoPlot_", con, ".pdf") 953 link_name <- paste0("VolcanoPlot_", con, ".pdf")
953 linkAddr <- paste0("volplot_", con, ".pdf") 954 link_addr <- paste0("volplot_", con, ".pdf")
954 linkData <- rbind(linkData, c(linkName, linkAddr)) 955 link_data <- rbind(link_data, c(link_name, link_addr))
955 invisible(dev.off()) 956 invisible(dev.off())
956 957
957 # PNG of MD and Volcano 958 # PNG of MD and Volcano
958 png(mdvolOutPng[i], width=1000, height=500) 959 png(mdvol_png[i], width = 1000, height = 500)
959 par(mfrow=c(1, 2), mar=c(5,4,2,2)+0.1, oma=c(0,0,3,0)) 960 par(mfrow = c(1, 2), mar = c(5, 4, 2, 2) + 0.1, oma = c(0, 0, 3, 0))
960 961
961 # MD plot 962 # MD plot
962 limma::plotMD(fit, status=status[, i], coef=i, main="MD Plot", 963 limma::plotMD(fit, status = status[, i], coef = i, main = "MD Plot",
963 hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1), 964 hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1),
964 xlab="Average Expression", ylab="logFC") 965 xlab = "Average Expression", ylab = "logFC")
965 abline(h=0, col="grey", lty=2) 966 abline(h = 0, col = "grey", lty = 2)
966 967
967 # Volcano 968 # Volcano
968 if (haveAnno) { 969 if (have_anno) {
969 # labels must be in second column currently 970 # labels must be in second column currently
970 limma::volcanoplot(fit, coef=i, main="Volcano Plot", 971 limma::volcanoplot(fit, coef = i, main = "Volcano Plot",
971 highlight=opt$topgenes, 972 highlight = opt$topgenes,
972 names=fit$genes[, 2]) 973 names = fit$genes[, 2])
973 } else { 974 } else {
974 limma::volcanoplot(fit, coef=i, main="Volcano Plot", 975 limma::volcanoplot(fit, coef = i, main = "Volcano Plot",
975 highlight=opt$topgenes, 976 highlight = opt$topgenes,
976 names=fit$genes$GeneID) 977 names = fit$genes$GeneID)
977 } 978 }
978 979
979 imgName <- paste0("MDVolPlot_", con) 980 img_name <- paste0("MDVolPlot_", con)
980 imgAddr <- paste0("mdvolplot_", con, ".png") 981 img_addr <- paste0("mdvolplot_", con, ".png")
981 imageData <- rbind(imageData, c(imgName, imgAddr)) 982 image_data <- rbind(image_data, c(img_name, img_addr))
982 title(paste0("Contrast: ", con_name), outer=TRUE, cex.main=1.5) 983 title(paste0("Contrast: ", con_name), outer = TRUE, cex.main = 1.5)
983 invisible(dev.off()) 984 invisible(dev.off())
984 985
985 if ("h" %in% plots) { 986 if ("h" %in% plots) {
986 # Plot Heatmap 987 # Plot Heatmap
987 topgenes <- rownames(top[1:opt$topgenes, ]) 988 topgenes <- rownames(top[1:opt$topgenes, ])
988 if (wantTrend) { 989 if (want_trend) {
989 topexp <- plotData[topgenes, ] 990 topexp <- plot_data[topgenes, ]
990 } else { 991 } else {
991 topexp <- plotData$E[topgenes, ] 992 topexp <- plot_data$E[topgenes, ]
992 } 993 }
993 pdf(heatOutPdf[i]) 994 pdf(heat_pdf[i])
994 mycol <- colorpanel(1000,"blue","white","red") 995 mycol <- colorpanel(1000, "blue", "white", "red")
995 if (haveAnno) { 996 if (have_anno) {
996 # labels must be in second column currently 997 # labels must be in second column currently
997 labels <- top[topgenes, 2] 998 labels <- top[topgenes, 2]
998 } else { 999 } else {
999 labels <- rownames(topexp) 1000 labels <- rownames(topexp)
1000 } 1001 }
1001 heatmap.2(topexp, scale="row", Colv=FALSE, Rowv=FALSE, dendrogram="none", 1002 heatmap.2(topexp, scale = "row", Colv = FALSE, Rowv = FALSE, dendrogram = "none",
1002 main=paste("Contrast:", unmake.names(con), "\nTop", opt$topgenes, "genes by adj.P.Val"), 1003 main = paste("Contrast:", unmake_names(con), "\nTop", opt$topgenes, "genes by adj.P.Val"),
1003 trace="none", density.info="none", lhei=c(2,10), margin=c(8, 6), labRow=labels, cexRow=0.7, srtCol=45, 1004 trace = "none", density.info = "none", lhei = c(2, 10), margin = c(8, 6), labRow = labels, cexRow = 0.7, srtCol = 45,
1004 col=mycol, ColSideColors=col.group) 1005 col = mycol, ColSideColors = col.group)
1005 linkName <- paste0("Heatmap_", con, ".pdf") 1006 link_name <- paste0("Heatmap_", con, ".pdf")
1006 linkAddr <- paste0("heatmap_", con, ".pdf") 1007 link_addr <- paste0("heatmap_", con, ".pdf")
1007 linkData <- rbind(linkData, c(linkName, linkAddr)) 1008 link_data <- rbind(link_data, c(link_name, link_addr))
1008 invisible(dev.off()) 1009 invisible(dev.off())
1009 } 1010 }
1010 1011
1011 if ("s" %in% plots) { 1012 if ("s" %in% plots) {
1012 # Plot Stripcharts of top genes 1013 # Plot Stripcharts of top genes
1013 pdf(stripOutPdf[i], title=paste("Contrast:", unmake.names(con))) 1014 pdf(strip_pdf[i], title = paste("Contrast:", unmake_names(con)))
1014 par(mfrow = c(3,2), cex.main=0.8, cex.axis=0.8) 1015 par(mfrow = c(3, 2), cex.main = 0.8, cex.axis = 0.8)
1015 cols <- unique(col.group) 1016 cols <- unique(col.group)
1016 1017
1017 for (j in 1:length(topgenes)) { 1018 for (j in seq_along(topgenes)) {
1018 lfc <- round(top[topgenes[j], "logFC"], 2) 1019 lfc <- round(top[topgenes[j], "logFC"], 2)
1019 pval <- round(top[topgenes[j], "adj.P.Val"], 5) 1020 pval <- round(top[topgenes[j], "adj.P.Val"], 5)
1020 if (wantTrend) { 1021 if (want_trend) {
1021 stripchart(plotData[topgenes[j], ] ~ factors[, 1], vertical=TRUE, las=2, pch=16, cex=0.8, cex.lab=0.8, col=cols, 1022 stripchart(plot_data[topgenes[j], ] ~ factors[, 1], vertical = TRUE, las = 2, pch = 16, cex = 0.8, cex.lab = 0.8, col = cols,
1022 method="jitter", ylab="Normalised log2 expression", main=paste0(labels[j], "\nlogFC=", lfc, ", adj.P.Val=", pval)) 1023 method = "jitter", ylab = "Normalised log2 expression", main = paste0(labels[j], "\nlogFC=", lfc, ", adj.P.Val=", pval))
1023 } else { 1024 } else {
1024 stripchart(plotData$E[topgenes[j], ] ~ factors[, 1], vertical=TRUE, las=2, pch=16, cex=0.8, cex.lab=0.8, col=cols, 1025 stripchart(plot_data$E[topgenes[j], ] ~ factors[, 1], vertical = TRUE, las = 2, pch = 16, cex = 0.8, cex.lab = 0.8, col = cols,
1025 method="jitter", ylab="Normalised log2 expression", main=paste0(labels[j], "\nlogFC=", lfc, ", adj.P.Val=", pval)) 1026 method = "jitter", ylab = "Normalised log2 expression", main = paste0(labels[j], "\nlogFC=", lfc, ", adj.P.Val=", pval))
1026 } 1027 }
1027 } 1028 }
1028 linkName <- paste0("Stripcharts_", con, ".pdf") 1029 link_name <- paste0("Stripcharts_", con, ".pdf")
1029 linkAddr <- paste0("stripcharts_", con, ".pdf") 1030 link_addr <- paste0("stripcharts_", con, ".pdf")
1030 linkData <- rbind(linkData, c(linkName, linkAddr)) 1031 link_data <- rbind(link_data, c(link_name, link_addr))
1031 invisible(dev.off()) 1032 invisible(dev.off())
1032 } 1033 }
1033 } 1034 }
1034 sigDiff <- data.frame(Up=upCount, Flat=flatCount, Down=downCount) 1035 sig_diff <- data.frame(Up = up_count, Flat = flat_count, Down = down_count)
1035 row.names(sigDiff) <- cons 1036 row.names(sig_diff) <- cons
1036 1037
1037 # Save relevant items as rda object 1038 # Save relevant items as rda object
1038 if (wantRda) { 1039 if (want_rda) {
1039 print("Saving RData") 1040 print("Saving RData")
1040 if (wantWeight) { 1041 if (want_weight) {
1041 save(counts, data, y, status, plotData, labels, factors, wts, fit, top, contrastData, contrasts, design, 1042 save(counts, data, y, status, plot_data, labels, factors, wts, fit, top, contrast_data, contrasts, design,
1042 file=rdaOut, ascii=TRUE) 1043 file = rda_out, ascii = TRUE)
1043 } else { 1044 } else {
1044 save(counts, data, y, status, plotData, labels, factors, fit, top, contrastData, contrasts, design, 1045 save(counts, data, y, status, plot_data, labels, factors, fit, top, contrast_data, contrasts, design,
1045 file=rdaOut, ascii=TRUE) 1046 file = rda_out, ascii = TRUE)
1046 } 1047 }
1047 linkData <- rbind(linkData, c((paste0(deMethod, "_analysis.RData")), (paste0(deMethod, "_analysis.RData")))) 1048 link_data <- rbind(link_data, c((paste0(de_method, "_analysis.RData")), (paste0(de_method, "_analysis.RData"))))
1048 } 1049 }
1049 1050
1050 # Record session info 1051 # Record session info
1051 writeLines(capture.output(sessionInfo()), sessionOut) 1052 writeLines(capture.output(sessionInfo()), session_out)
1052 linkData <- rbind(linkData, c("Session Info", "session_info.txt")) 1053 link_data <- rbind(link_data, c("Session Info", "session_info.txt"))
1053 1054
1054 # Record ending time and calculate total run time 1055 # Record ending time and calculate total run time
1055 timeEnd <- as.character(Sys.time()) 1056 time_end <- as.character(Sys.time())
1056 timeTaken <- capture.output(round(difftime(timeEnd,timeStart), digits=3)) 1057 time_taken <- capture.output(round(difftime(time_end, time_start), digits = 3))
1057 timeTaken <- gsub("Time difference of ", "", timeTaken, fixed=TRUE) 1058 time_taken <- gsub("Time difference of ", "", time_taken, fixed = TRUE)
1058 ################################################################################ 1059 ################################################################################
1059 ### HTML Generation 1060 ### HTML Generation
1060 ################################################################################ 1061 ################################################################################
1061 1062
1062 # Clear file 1063 # Clear file
1063 cat("", file=opt$htmlPath) 1064 cat("", file = opt$htmlPath)
1064 1065
1065 cata("<html>\n") 1066 cata("<html>\n")
1066 1067
1067 cata("<body>\n") 1068 cata("<body>\n")
1068 cata("<h3>Limma Analysis Output:</h3>\n") 1069 cata("<h3>Limma Analysis Output:</h3>\n")
1069 cata("Links to PDF copies of plots are in 'Plots' section below <br />\n") 1070 cata("Links to PDF copies of plots are in 'Plots' section below <br />\n")
1070 1071
1071 for (i in 1:nrow(imageData)) { 1072 for (i in seq_len(nrow(image_data))) {
1072 if (grepl("density|box|mds|mdvol|wts", imageData$Link[i])) { 1073 if (grepl("density|box|mds|mdvol|wts", image_data$Link[i])) {
1073 HtmlImage(imageData$Link[i], imageData$Label[i], width=1000) 1074 html_image(image_data$Link[i], image_data$Label[i], width = 1000)
1074 } else { 1075 } else {
1075 HtmlImage(imageData$Link[i], imageData$Label[i]) 1076 html_image(image_data$Link[i], image_data$Label[i])
1076 } 1077 }
1077 } 1078 }
1078 1079
1079 cata("<h4>Differential Expression Counts:</h4>\n") 1080 cata("<h4>Differential Expression Counts:</h4>\n")
1080 1081
1081 cata("<table border=\"1\" cellpadding=\"4\">\n") 1082 cata("<table border=\"1\" cellpadding=\"4\">\n")
1082 cata("<tr>\n") 1083 cata("<tr>\n")
1083 TableItem() 1084 table_item()
1084 for (i in colnames(sigDiff)) { 1085 for (i in colnames(sig_diff)) {
1085 TableHeadItem(i) 1086 table_head_item(i)
1086 } 1087 }
1087 cata("</tr>\n") 1088 cata("</tr>\n")
1088 for (i in 1:nrow(sigDiff)) { 1089 for (i in seq_len(nrow(sig_diff))) {
1089 cata("<tr>\n") 1090 cata("<tr>\n")
1090 TableHeadItem(unmake.names(row.names(sigDiff)[i])) 1091 table_head_item(unmake_names(row.names(sig_diff)[i]))
1091 for (j in 1:ncol(sigDiff)) { 1092 for (j in seq_len(ncol(sig_diff))) {
1092 TableItem(as.character(sigDiff[i, j])) 1093 table_item(as.character(sig_diff[i, j]))
1093 } 1094 }
1094 cata("</tr>\n") 1095 cata("</tr>\n")
1095 } 1096 }
1096 cata("</table>") 1097 cata("</table>")
1097 1098
1098 cata("<h4>Plots:</h4>\n") 1099 cata("<h4>Plots:</h4>\n")
1099 #PDFs 1100 #PDFs
1100 for (i in 1:nrow(linkData)) { 1101 for (i in seq_len(nrow(link_data))) {
1101 if (grepl(".pdf", linkData$Link[i]) & grepl("density|cpm|boxplot|mds|mdplots|voom|saplot", linkData$Link[i])) { 1102 if (grepl(".pdf", link_data$Link[i]) & grepl("density|cpm|boxplot|mds|mdplots|voom|saplot", link_data$Link[i])) {
1102 HtmlLink(linkData$Link[i], linkData$Label[i]) 1103 html_link(link_data$Link[i], link_data$Label[i])
1103 } 1104 }
1104 } 1105 }
1105 1106
1106 for (i in 1:nrow(linkData)) { 1107 for (i in seq_len(nrow(link_data))) {
1107 if (grepl("mdplot_", linkData$Link[i])) { 1108 if (grepl("mdplot_", link_data$Link[i])) {
1108 HtmlLink(linkData$Link[i], linkData$Label[i]) 1109 html_link(link_data$Link[i], link_data$Label[i])
1109 } 1110 }
1110 } 1111 }
1111 1112
1112 for (i in 1:nrow(linkData)) { 1113 for (i in seq_len(nrow(link_data))) {
1113 if (grepl("volplot", linkData$Link[i])) { 1114 if (grepl("volplot", link_data$Link[i])) {
1114 HtmlLink(linkData$Link[i], linkData$Label[i]) 1115 html_link(link_data$Link[i], link_data$Label[i])
1115 } 1116 }
1116 } 1117 }
1117 1118
1118 for (i in 1:nrow(linkData)) { 1119 for (i in seq_len(nrow(link_data))) {
1119 if (grepl("heatmap", linkData$Link[i])) { 1120 if (grepl("heatmap", link_data$Link[i])) {
1120 HtmlLink(linkData$Link[i], linkData$Label[i]) 1121 html_link(link_data$Link[i], link_data$Label[i])
1121 } 1122 }
1122 } 1123 }
1123 1124
1124 for (i in 1:nrow(linkData)) { 1125 for (i in seq_len(nrow(link_data))) {
1125 if (grepl("stripcharts", linkData$Link[i])) { 1126 if (grepl("stripcharts", link_data$Link[i])) {
1126 HtmlLink(linkData$Link[i], linkData$Label[i]) 1127 html_link(link_data$Link[i], link_data$Label[i])
1127 } 1128 }
1128 } 1129 }
1129 1130
1130 cata("<h4>Tables:</h4>\n") 1131 cata("<h4>Tables:</h4>\n")
1131 for (i in 1:nrow(linkData)) { 1132 for (i in seq_len(nrow(link_data))) {
1132 if (grepl("counts$", linkData$Link[i])) { 1133 if (grepl("counts$", link_data$Link[i])) {
1133 HtmlLink(linkData$Link[i], linkData$Label[i]) 1134 html_link(link_data$Link[i], link_data$Label[i])
1134 } else if (grepl(".tsv", linkData$Link[i])) { 1135 } else if (grepl(".tsv", link_data$Link[i])) {
1135 HtmlLink(linkData$Link[i], linkData$Label[i]) 1136 html_link(link_data$Link[i], link_data$Label[i])
1136 } 1137 }
1137 } 1138 }
1138 1139
1139 if (wantRda) { 1140 if (want_rda) {
1140 cata("<h4>R Data Object:</h4>\n") 1141 cata("<h4>R Data Object:</h4>\n")
1141 for (i in 1:nrow(linkData)) { 1142 for (i in seq_len(nrow(link_data))) {
1142 if (grepl(".RData", linkData$Link[i])) { 1143 if (grepl(".RData", link_data$Link[i])) {
1143 HtmlLink(linkData$Link[i], linkData$Label[i]) 1144 html_link(link_data$Link[i], link_data$Label[i])
1144 } 1145 }
1145 } 1146 }
1146 } 1147 }
1147 1148
1148 if ("i" %in% plots) { 1149 if ("i" %in% plots) {
1149 cata("<h4>Glimma Interactive Results:</h4>\n") 1150 cata("<h4>Glimma Interactive Results:</h4>\n")
1150 for (i in 1:nrow(linkData)) { 1151 for (i in seq_len(nrow(link_data))) {
1151 if (grepl("glimma", linkData$Link[i])) { 1152 if (grepl("glimma", link_data$Link[i])) {
1152 HtmlLink(linkData$Link[i], linkData$Label[i]) 1153 html_link(link_data$Link[i], link_data$Label[i])
1153 } 1154 }
1154 } 1155 }
1155 } 1156 }
1156 1157
1157 cata("<p>Alt-click links to download file.</p>\n") 1158 cata("<p>Alt-click links to download file.</p>\n")
1160 cata("<p>.tsv files can be viewed in Excel or any spreadsheet program.</p>\n") 1161 cata("<p>.tsv files can be viewed in Excel or any spreadsheet program.</p>\n")
1161 1162
1162 cata("<h4>Additional Information</h4>\n") 1163 cata("<h4>Additional Information</h4>\n")
1163 cata("<ul>\n") 1164 cata("<ul>\n")
1164 1165
1165 if (filtCPM || filtSmpCount || filtTotCount) { 1166 if (filt_cpm || filt_smpcount || filt_totcount) {
1166 if (filtCPM) { 1167 if (filt_cpm) {
1167 tempStr <- paste("Genes without more than", opt$cpmReq, 1168 temp_str <- paste("Genes without more than", opt$cpmReq,
1168 "CPM in at least", opt$sampleReq, "samples are insignificant", 1169 "CPM in at least", opt$sampleReq, "samples are insignificant",
1169 "and filtered out.") 1170 "and filtered out.")
1170 } else if (filtSmpCount) { 1171 } else if (filt_smpcount) {
1171 tempStr <- paste("Genes without more than", opt$cntReq, 1172 temp_str <- paste("Genes without more than", opt$cntReq,
1172 "counts in at least", opt$sampleReq, "samples are insignificant", 1173 "counts in at least", opt$sampleReq, "samples are insignificant",
1173 "and filtered out.") 1174 "and filtered out.")
1174 } else if (filtTotCount) { 1175 } else if (filt_totcount) {
1175 tempStr <- paste("Genes without more than", opt$cntReq, 1176 temp_str <- paste("Genes without more than", opt$cntReq,
1176 "counts, after summing counts for all samples, are insignificant", 1177 "counts, after summing counts for all samples, are insignificant",
1177 "and filtered out.") 1178 "and filtered out.")
1178 } 1179 }
1179 1180
1180 ListItem(tempStr) 1181 list_item(temp_str)
1181 filterProp <- round(filteredCount/preFilterCount*100, digits=2) 1182 filter_prop <- round(filtered_count / prefilter_count * 100, digits = 2)
1182 tempStr <- paste0(filteredCount, " of ", preFilterCount," (", filterProp, 1183 temp_str <- paste0(filtered_count, " of ", prefilter_count, " (", filter_prop,
1183 "%) genes were filtered out for low expression.") 1184 "%) genes were filtered out for low expression.")
1184 ListItem(tempStr) 1185 list_item(temp_str)
1185 } 1186 }
1186 ListItem(opt$normOpt, " was the method used to normalise library sizes.") 1187 list_item(opt$normOpt, " was the method used to normalise library sizes.")
1187 if (wantTrend) { 1188 if (want_trend) {
1188 ListItem("The limma-trend method was used.") 1189 list_item("The limma-trend method was used.")
1189 } else { 1190 } else {
1190 ListItem("The limma-voom method was used.") 1191 list_item("The limma-voom method was used.")
1191 } 1192 }
1192 if (wantWeight) { 1193 if (want_weight) {
1193 ListItem("Weights were applied to samples.") 1194 list_item("Weights were applied to samples.")
1194 } else { 1195 } else {
1195 ListItem("Weights were not applied to samples.") 1196 list_item("Weights were not applied to samples.")
1196 } 1197 }
1197 if (wantTreat) { 1198 if (want_treat) {
1198 ListItem(paste("Testing significance relative to a fold-change threshold (TREAT) was performed using a threshold of log2 =", opt$lfcReq, "at FDR of", opt$pValReq, ".")) 1199 list_item(paste("Testing significance relative to a fold-change threshold (TREAT) was performed using a threshold of log2 =", opt$lfcReq, "at FDR of", opt$pValReq, "."))
1199 } 1200 }
1200 if (wantRobust) { 1201 if (want_robust) {
1201 if (wantTreat) { 1202 if (want_treat) {
1202 ListItem("TREAT was used with robust settings (robust=TRUE).") 1203 list_item("TREAT was used with robust settings (robust = TRUE).")
1203 } else { 1204 } else {
1204 ListItem("eBayes was used with robust settings (robust=TRUE).") 1205 list_item("eBayes was used with robust settings (robust = TRUE).")
1205 } 1206 }
1206 } 1207 }
1207 if (opt$pAdjOpt!="none") { 1208 if (opt$pAdjOpt != "none") {
1208 if (opt$pAdjOpt=="BH" || opt$pAdjOpt=="BY") { 1209 if (opt$pAdjOpt == "BH" || opt$pAdjOpt == "BY") {
1209 tempStr <- paste0("MD Plot highlighted genes are significant at FDR ", 1210 temp_str <- paste0("MD Plot highlighted genes are significant at FDR ",
1210 "of ", opt$pValReq," and exhibit log2-fold-change of at ", 1211 "of ", opt$pValReq, " and exhibit log2-fold-change of at ",
1211 "least ", opt$lfcReq, ".") 1212 "least ", opt$lfcReq, ".")
1212 ListItem(tempStr) 1213 list_item(temp_str)
1213 } else if (opt$pAdjOpt=="holm") { 1214 } else if (opt$pAdjOpt == "holm") {
1214 tempStr <- paste0("MD Plot highlighted genes are significant at adjusted ", 1215 temp_str <- paste0("MD Plot highlighted genes are significant at adjusted ",
1215 "p-value of ", opt$pValReq," by the Holm(1979) ", 1216 "p-value of ", opt$pValReq, " by the Holm(1979) ",
1216 "method, and exhibit log2-fold-change of at least ", 1217 "method, and exhibit log2-fold-change of at least ",
1217 opt$lfcReq, ".") 1218 opt$lfcReq, ".")
1218 ListItem(tempStr) 1219 list_item(temp_str)
1219 } 1220 }
1220 } else { 1221 } else {
1221 tempStr <- paste0("MD Plot highlighted genes are significant at p-value ", 1222 temp_str <- paste0("MD Plot highlighted genes are significant at p-value ",
1222 "of ", opt$pValReq," and exhibit log2-fold-change of at ", 1223 "of ", opt$pValReq, " and exhibit log2-fold-change of at ",
1223 "least ", opt$lfcReq, ".") 1224 "least ", opt$lfcReq, ".")
1224 ListItem(tempStr) 1225 list_item(temp_str)
1225 } 1226 }
1226 cata("</ul>\n") 1227 cata("</ul>\n")
1227 1228
1228 cata("<h4>Summary of experimental data:</h4>\n") 1229 cata("<h4>Summary of experimental data:</h4>\n")
1229 1230
1230 cata("<p>*CHECK THAT SAMPLES ARE ASSOCIATED WITH CORRECT GROUP(S)*</p>\n") 1231 cata("<p>*CHECK THAT SAMPLES ARE ASSOCIATED WITH CORRECT GROUP(S)*</p>\n")
1231 1232
1232 cata("<table border=\"1\" cellpadding=\"3\">\n") 1233 cata("<table border=\"1\" cellpadding=\"3\">\n")
1233 cata("<tr>\n") 1234 cata("<tr>\n")
1234 TableHeadItem("SampleID") 1235 table_head_item("SampleID")
1235 TableHeadItem(names(factors)[1]," (Primary Factor)") 1236 table_head_item(names(factors)[1], " (Primary Factor)")
1236 1237
1237 if (ncol(factors) > 1) { 1238 if (ncol(factors) > 1) {
1238 for (i in names(factors)[2:length(names(factors))]) { 1239 for (i in names(factors)[2:length(names(factors))]) {
1239 TableHeadItem(i) 1240 table_head_item(i)
1240 } 1241 }
1241 cata("</tr>\n") 1242 cata("</tr>\n")
1242 } 1243 }
1243 1244
1244 for (i in 1:nrow(factors)) { 1245 for (i in seq_len(nrow(factors))) {
1245 cata("<tr>\n") 1246 cata("<tr>\n")
1246 TableHeadItem(row.names(factors)[i]) 1247 table_head_item(row.names(factors)[i])
1247 for (j in 1:ncol(factors)) { 1248 for (j in seq_len(ncol(factors))) {
1248 TableItem(as.character(unmake.names(factors[i, j]))) 1249 table_item(as.character(unmake_names(factors[i, j])))
1249 } 1250 }
1250 cata("</tr>\n") 1251 cata("</tr>\n")
1251 } 1252 }
1252 cata("</table>") 1253 cata("</table>")
1253 1254
1311 cata(cit[2], "\n") 1312 cata(cit[2], "\n")
1312 1313
1313 cata("<h4>limma</h4>\n") 1314 cata("<h4>limma</h4>\n")
1314 cata(cit[3], "\n") 1315 cata(cit[3], "\n")
1315 cata("<ol>\n") 1316 cata("<ol>\n")
1316 ListItem(cit[4]) 1317 list_item(cit[4])
1317 ListItem(cit[10]) 1318 list_item(cit[10])
1318 ListItem(cit[11]) 1319 list_item(cit[11])
1319 cata("</ol>\n") 1320 cata("</ol>\n")
1320 1321
1321 cata("<h4>edgeR</h4>\n") 1322 cata("<h4>edgeR</h4>\n")
1322 cata(cit[5], "\n") 1323 cata(cit[5], "\n")
1323 cata("<ol>\n") 1324 cata("<ol>\n")
1324 ListItem(cit[6]) 1325 list_item(cit[6])
1325 ListItem(cit[7]) 1326 list_item(cit[7])
1326 ListItem(cit[8]) 1327 list_item(cit[8])
1327 ListItem(cit[9]) 1328 list_item(cit[9])
1328 cata("</ol>\n") 1329 cata("</ol>\n")
1329 1330
1330 cata("<p>Please report problems or suggestions to: su.s@wehi.edu.au</p>\n") 1331 cata("<p>Please report problems or suggestions to: su.s@wehi.edu.au</p>\n")
1331 1332
1332 for (i in 1:nrow(linkData)) { 1333 for (i in seq_len(nrow(link_data))) {
1333 if (grepl("session_info", linkData$Link[i])) { 1334 if (grepl("session_info", link_data$Link[i])) {
1334 HtmlLink(linkData$Link[i], linkData$Label[i]) 1335 html_link(link_data$Link[i], link_data$Label[i])
1335 } 1336 }
1336 } 1337 }
1337 1338
1338 cata("<table border=\"0\">\n") 1339 cata("<table border=\"0\">\n")
1339 cata("<tr>\n") 1340 cata("<tr>\n")
1340 TableItem("Task started at:"); TableItem(timeStart) 1341 table_item("Task started at:"); table_item(time_start)
1341 cata("</tr>\n") 1342 cata("</tr>\n")
1342 cata("<tr>\n") 1343 cata("<tr>\n")
1343 TableItem("Task ended at:"); TableItem(timeEnd) 1344 table_item("Task ended at:"); table_item(time_end)
1344 cata("</tr>\n") 1345 cata("</tr>\n")
1345 cata("<tr>\n") 1346 cata("<tr>\n")
1346 TableItem("Task run time:"); TableItem(timeTaken) 1347 table_item("Task run time:"); table_item(time_taken)
1347 cata("<tr>\n") 1348 cata("<tr>\n")
1348 cata("</table>\n") 1349 cata("</table>\n")
1349 1350
1350 cata("</body>\n") 1351 cata("</body>\n")
1351 cata("</html>") 1352 cata("</html>")