Mercurial > repos > iuc > limma_voom

comparison limma_voom.R @ 26:119b069fc845 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/galaxyproject/tools-iuc/tree/master/tools/limma_voom commit b4d788c0f159d507159117a063b1f867b243c738
author iuc
date Fri, 09 Feb 2024 17:06:25 +0000
parents 708348a17fa1
children
comparison
equal deleted inserted replaced
25:d6f5fa4ee473 26:119b069fc845
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 time_start <- as.character(Sys.time()) 49 time_start <- Sys.time()
50
51 # Setup R error handling to go to stderr
52 options(
53 show.error.messages = FALSE,
54 error = function() {
55 cat(geterrmessage(), file = stderr())
56 q("no", 1, FALSE)
57 }
58 )
59
60 # Unify locale settings
61 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
62
63 warnings()
50 64
51 # Load all required libraries 65 # Load all required libraries
52 library(methods, quietly = TRUE, warn.conflicts = FALSE) 66 library(methods, quietly = TRUE, warn.conflicts = FALSE)
53 library(statmod, quietly = TRUE, warn.conflicts = FALSE) 67 library(statmod, quietly = TRUE, warn.conflicts = FALSE)
54 library(splines, quietly = TRUE, warn.conflicts = FALSE) 68 library(splines, quietly = TRUE, warn.conflicts = FALSE)
104 } 118 }
105 119
106 # Create cata function: default path set, default seperator empty and appending 120 # Create cata function: default path set, default seperator empty and appending
107 # true by default (Ripped straight from the cat function with altered argument 121 # true by default (Ripped straight from the cat function with altered argument
108 # defaults) 122 # defaults)
109 cata <- function(..., file = opt$htmlPath, sep = "", fill = FALSE, labels = NULL, 123 cata <- function(..., file = opt$htmlPath, sep = "", fill = FALSE,
110 append = TRUE) { 124 labels = NULL, append = TRUE) {
111 if (is.character(file)) 125 if (is.character(file)) {
112 if (file == "") 126 if (file == "") {
113 file <- stdout() 127 file <- stdout()
114 else if (substring(file, 1L, 1L) == "|") { 128 } else if (substring(file, 1L, 1L) == "|") {
115 file <- pipe(substring(file, 2L), "w") 129 file <- pipe(substring(file, 2L), "w")
116 on.exit(close(file)) 130 on.exit(close(file))
117 } 131 } else {
118 else { 132 file <- file(file, ifelse(append, "a", "w"))
119 file <- file(file, ifelse(append, "a", "w")) 133 on.exit(close(file))
120 on.exit(close(file)) 134 }
121 } 135 }
122 .Internal(cat(list(...), file, sep, fill, labels, append)) 136 .Internal(cat(list(...), file, sep, fill, labels, append))
123 } 137 }
124 138
125 # Function to write code for html head and title 139 # Function to write code for html head and title
126 html_head <- function(title) { 140 html_head <- function(title) {
160 # Collect arguments from command line 174 # Collect arguments from command line
161 args <- commandArgs(trailingOnly = TRUE) 175 args <- commandArgs(trailingOnly = TRUE)
162 176
163 # Get options, using the spec as defined by the enclosed list. 177 # Get options, using the spec as defined by the enclosed list.
164 # Read the options from the default: commandArgs(TRUE). 178 # Read the options from the default: commandArgs(TRUE).
165 spec <- matrix(c( 179 spec <- matrix(
166 "htmlPath", "R", 1, "character", 180 c(
167 "outPath", "o", 1, "character", 181 "htmlPath", "R", 1, "character",
168 "filesPath", "j", 2, "character", 182 "outPath", "o", 1, "character",
169 "matrixPath", "m", 2, "character", 183 "filesPath", "j", 2, "character",
170 "factFile", "f", 2, "character", 184 "matrixPath", "m", 2, "character",
171 "factInput", "i", 2, "character", 185 "factFile", "f", 2, "character",
172 "annoPath", "a", 2, "character", 186 "factInput", "i", 2, "character",
173 "contrastFile", "C", 1, "character", 187 "annoPath", "a", 2, "character",
174 "contrastInput", "D", 1, "character", 188 "contrastFile", "C", 1, "character",
175 "cpmReq", "c", 1, "double", 189 "contrastInput", "D", 1, "character",
176 "totReq", "y", 0, "logical", 190 "cpmReq", "c", 1, "double",
177 "cntReq", "z", 1, "integer", 191 "totReq", "y", 0, "logical",
178 "sampleReq", "s", 1, "integer", 192 "cntReq", "z", 1, "integer",
179 "filtCounts", "F", 0, "logical", 193 "sampleReq", "s", 1, "integer",
180 "normCounts", "x", 0, "logical", 194 "filtCounts", "F", 0, "logical",
181 "rdaOpt", "r", 0, "logical", 195 "normCounts", "x", 0, "logical",
182 "lfcReq", "l", 1, "double", 196 "rdaOpt", "r", 0, "logical",
183 "pValReq", "p", 1, "double", 197 "lfcReq", "l", 1, "double",
184 "pAdjOpt", "d", 1, "character", 198 "pValReq", "p", 1, "double",
185 "normOpt", "n", 1, "character", 199 "pAdjOpt", "d", 1, "character",
186 "robOpt", "b", 0, "logical", 200 "normOpt", "n", 1, "character",
187 "trend", "t", 1, "double", 201 "robOpt", "b", 0, "logical",
188 "weightOpt", "w", 0, "logical", 202 "trend", "t", 1, "double",
189 "topgenes", "G", 1, "integer", 203 "weightOpt", "w", 0, "logical",
190 "treatOpt", "T", 0, "logical", 204 "topgenes", "G", 1, "integer",
191 "plots", "P", 1, "character", 205 "treatOpt", "T", 0, "logical",
192 "libinfoOpt", "L", 0, "logical"), 206 "plots", "P", 1, "character",
193 byrow = TRUE, ncol = 4) 207 "libinfoOpt", "L", 0, "logical"
208 ),
209 byrow = TRUE, ncol = 4
210 )
194 opt <- getopt(spec) 211 opt <- getopt(spec)
195 212
196 213
197 if (is.null(opt$matrixPath) & is.null(opt$filesPath)) { 214 if (is.null(opt$matrixPath) && is.null(opt$filesPath)) {
198 cat("A counts matrix (or a set of counts files) is required.\n") 215 cat("A counts matrix (or a set of counts files) is required.\n")
199 q(status = 1) 216 q(status = 1)
200 } 217 }
201 218
202 if (is.null(opt$cpmReq)) { 219 if (is.null(opt$cpmReq)) {
281 parser <- newJSONParser() 298 parser <- newJSONParser()
282 parser$addData(opt$filesPath) 299 parser$addData(opt$filesPath)
283 factor_list <- parser$getObject() 300 factor_list <- parser$getObject()
284 factors <- sapply(factor_list, function(x) x[[1]]) 301 factors <- sapply(factor_list, function(x) x[[1]])
285 filenames_in <- unname(unlist(factor_list[[1]][[2]])) 302 filenames_in <- unname(unlist(factor_list[[1]][[2]]))
286 sampletable <- data.frame(sample = basename(filenames_in), 303 sampletable <- data.frame(
287 filename = filenames_in, 304 sample = basename(filenames_in),
288 row.names = filenames_in, 305 filename = filenames_in,
289 stringsAsFactors = FALSE) 306 row.names = filenames_in,
307 stringsAsFactors = FALSE
308 )
290 for (factor in factor_list) { 309 for (factor in factor_list) {
291 factorname <- factor[[1]] 310 factorname <- factor[[1]]
292 sampletable[[factorname]] <- character(nrow(sampletable)) 311 sampletable[[factorname]] <- character(nrow(sampletable))
293 lvls <- sapply(factor[[2]], function(x) names(x)) 312 lvls <- sapply(factor[[2]], function(x) names(x))
294 for (i in seq_along(factor[[2]])) { 313 for (i in seq_along(factor[[2]])) {
299 } 318 }
300 rownames(sampletable) <- sampletable$sample 319 rownames(sampletable) <- sampletable$sample
301 rem <- c("sample", "filename") 320 rem <- c("sample", "filename")
302 factors <- sampletable[, !(names(sampletable) %in% rem), drop = FALSE] 321 factors <- sampletable[, !(names(sampletable) %in% rem), drop = FALSE]
303 322
304 #read in count files and create single table 323 # read in count files and create single table
305 countfiles <- lapply(sampletable$filename, function(x) { 324 countfiles <- lapply(sampletable$filename, function(x) {
306 read.delim(x, row.names = 1) 325 read.delim(x, row.names = 1)
307 }) 326 })
308 counts <- do.call("cbind", countfiles) 327 counts <- do.call("cbind", countfiles)
309
310 } else { 328 } else {
311 # Process the single count matrix 329 # Process the single count matrix
312 counts <- read.table(opt$matrixPath, header = TRUE, sep = "\t", strip.white = TRUE, stringsAsFactors = FALSE, check.names = FALSE) 330 counts <- read.table(opt$matrixPath, header = TRUE, sep = "\t", strip.white = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
313 row.names(counts) <- counts[, 1] 331 row.names(counts) <- counts[, 1]
314 counts <- counts[, -1] 332 counts <- counts[, -1]
315 countsrows <- nrow(counts) 333 countsrows <- nrow(counts)
316 334
317 # Process factors 335 # Process factors
318 if (is.null(opt$factInput)) { 336 if (is.null(opt$factInput)) {
319 factordata <- read.table(opt$factFile, header = TRUE, sep = "\t", strip.white = TRUE, stringsAsFactors = TRUE) 337 factordata <- read.table(opt$factFile, header = TRUE, sep = "\t", strip.white = TRUE, stringsAsFactors = TRUE)
320 if (!setequal(factordata[, 1], colnames(counts))) 338 if (!setequal(factordata[, 1], colnames(counts))) {
321 stop("Sample IDs in counts and factors files don't match") 339 stop("Sample IDs in counts and factors files don't match")
340 }
322 # order samples as in counts matrix 341 # order samples as in counts matrix
323 factordata <- factordata[match(colnames(counts), factordata[, 1]), ] 342 factordata <- factordata[match(colnames(counts), factordata[, 1]), ]
324 factors <- factordata[, -1, drop = FALSE] 343 factors <- factordata[, -1, drop = FALSE]
325 } else { 344 } else {
326 factors <- unlist(strsplit(opt$factInput, "|", fixed = TRUE)) 345 factors <- unlist(strsplit(opt$factInput, "|", fixed = TRUE))
327 factordata <- list() 346 factordata <- list()
328 for (fact in factors) { 347 for (fact in factors) {
329 newfact <- unlist(strsplit(fact, split = "::")) 348 newfact <- unlist(strsplit(fact, split = "::"))
330 factordata <- rbind(factordata, newfact) 349 factordata <- rbind(factordata, newfact)
345 } 364 }
346 # make groups valid R names, required for makeContrasts 365 # make groups valid R names, required for makeContrasts
347 factors <- sapply(factors, make.names) 366 factors <- sapply(factors, make.names)
348 factors <- data.frame(factors, stringsAsFactors = TRUE) 367 factors <- data.frame(factors, stringsAsFactors = TRUE)
349 368
350 # if annotation file provided 369 # if annotation file provided
351 if (have_anno) { 370 if (have_anno) {
352 geneanno <- read.table(opt$annoPath, header = TRUE, sep = "\t", quote = "", strip.white = TRUE, stringsAsFactors = FALSE) 371 geneanno <- read.table(opt$annoPath, header = TRUE, sep = "\t", quote = "", strip.white = TRUE, stringsAsFactors = FALSE)
353 } 372 }
354 373
355 #Create output directory 374 # Create output directory
356 dir.create(opt$outPath, showWarnings = FALSE) 375 dir.create(opt$outPath, showWarnings = FALSE)
357 376
358 # Process contrasts 377 # Process contrasts
359 if (is.null(opt$contrastInput)) { 378 if (is.null(opt$contrastInput)) {
360 contrast_data <- read.table(opt$contrastFile, header = TRUE, sep = "\t", quote = "", strip.white = TRUE, stringsAsFactors = FALSE) 379 contrast_data <- read.table(opt$contrastFile, header = TRUE, sep = "\t", quote = "", strip.white = TRUE, stringsAsFactors = FALSE)
361 contrast_data <- contrast_data[, 1, drop = TRUE] 380 contrast_data <- contrast_data[, 1, drop = TRUE]
362 } else { 381 } else {
363 # Split up contrasts seperated by comma into a vector then sanitise 382 # Split up contrasts seperated by comma into a vector then sanitise
364 contrast_data <- unlist(strsplit(opt$contrastInput, split = ",")) 383 contrast_data <- unlist(strsplit(opt$contrastInput, split = ","))
365 } 384 }
366 contrast_data <- sanitise_equation(contrast_data) 385 contrast_data <- sanitise_equation(contrast_data)
367 contrast_data <- gsub(" ", ".", contrast_data, fixed = TRUE) 386 contrast_data <- gsub(" ", ".", contrast_data, fixed = TRUE)
368 387
369 # in case input groups start with numbers make the names valid R names, required for makeContrasts 388 # in case input groups start with numbers make the names valid R names, required for makeContrasts
370 cons <- NULL 389 cons <- NULL
371 cons_d <- NULL 390 cons_d <- NULL
372 for (i in contrast_data) { 391 for (i in contrast_data) {
373
374 # if the contrast is a difference of differences e.g. (A-B)-(X-Y) 392 # if the contrast is a difference of differences e.g. (A-B)-(X-Y)
375 if (grepl("\\)-\\(", i)) { 393 if (grepl("\\)-\\(", i)) {
376 i <- unlist(strsplit(i, split = "\\)-\\(")) 394 i <- unlist(strsplit(i, split = "\\)-\\("))
377 i <- gsub("\\(|\\)", "", i) 395 i <- gsub("\\(|\\)", "", i)
378 for (j in i) { 396 for (j in i) {
379 j <- sanitise_contrast(j) 397 j <- sanitise_contrast(j)
380 j <- paste0("(", j, ")") 398 j <- paste0("(", j, ")")
381 cons_d <- append(cons_d, unlist(j)) 399 cons_d <- append(cons_d, unlist(j))
382 } 400 }
383 cons_d <- paste(cons_d, collapse = "-") 401 cons_d <- paste(cons_d, collapse = "-")
384 cons <- append(cons, unlist(cons_d)) 402 cons <- append(cons, unlist(cons_d))
385 } else { 403 } else {
386 i <- sanitise_contrast(i) 404 i <- sanitise_contrast(i)
426 rda_out <- make_out(paste0(de_method, "_analysis.RData")) 444 rda_out <- make_out(paste0(de_method, "_analysis.RData"))
427 session_out <- make_out("session_info.txt") 445 session_out <- make_out("session_info.txt")
428 446
429 # Initialise data for html links and images, data frame with columns Label and 447 # Initialise data for html links and images, data frame with columns Label and
430 # Link 448 # Link
431 link_data <- data.frame(Label = character(), Link = character(), 449 link_data <- data.frame(
432 stringsAsFactors = FALSE) 450 Label = character(), Link = character(),
433 image_data <- data.frame(Label = character(), Link = character(), 451 stringsAsFactors = FALSE
434 stringsAsFactors = FALSE) 452 )
453 image_data <- data.frame(
454 Label = character(), Link = character(),
455 stringsAsFactors = FALSE
456 )
435 457
436 # Initialise vectors for storage of up/down/neutral regulated counts 458 # Initialise vectors for storage of up/down/neutral regulated counts
437 up_count <- numeric() 459 up_count <- numeric()
438 down_count <- numeric() 460 down_count <- numeric()
439 flat_count <- numeric() 461 flat_count <- numeric()
445 # Extract counts and annotation data 467 # Extract counts and annotation data
446 print("Extracting counts") 468 print("Extracting counts")
447 data <- list() 469 data <- list()
448 data$counts <- counts 470 data$counts <- counts
449 if (have_anno) { 471 if (have_anno) {
450 # order annotation by genes in counts (assumes gene ids are in 1st column of geneanno) 472 # order annotation by genes in counts (assumes gene ids are in 1st column of geneanno)
451 annoord <- geneanno[match(row.names(counts), geneanno[, 1]), ] 473 annoord <- geneanno[match(row.names(counts), geneanno[, 1]), ]
452 data$genes <- annoord 474 data$genes <- annoord
453 } else { 475 } else {
454 data$genes <- data.frame(GeneID = row.names(counts)) 476 data$genes <- data.frame(GeneID = row.names(counts))
455 } 477 }
456 478
457 # Creating naming data 479 # Creating naming data
458 samplenames <- colnames(data$counts) 480 samplenames <- colnames(data$counts)
459 sampleanno <- data.frame("sampleID" = samplenames, factors) 481 sampleanno <- data.frame("sampleID" = samplenames, factors)
460 row.names(factors) <- samplenames # for "Summary of experimental data" table 482 row.names(factors) <- samplenames # for "Summary of experimental data" table
461 483
462 # Creating colours for the groups 484 # Creating colours for the groups
463 cols <- as.numeric(factors[, 1]) 485 cols <- as.numeric(factors[, 1])
464 col.group <- palette()[cols] 486 col_group <- palette()[cols]
465 487
466 # If filter crieteria set, filter out genes that do not have a required cpm/counts in a required number of 488 # If filter crieteria set, filter out genes that do not have a required cpm/counts in a required number of
467 # samples. Default is no filtering 489 # samples. Default is no filtering
468 prefilter_count <- nrow(data$counts) 490 prefilter_count <- nrow(data$counts)
469 nsamples <- ncol(data$counts) 491 nsamples <- ncol(data$counts)
470 492
471 if (filt_cpm || filt_smpcount || filt_totcount) { 493 if (filt_cpm || filt_smpcount || filt_totcount) {
472
473 if (filt_totcount) { 494 if (filt_totcount) {
474 keep <- rowSums(data$counts) >= opt$cntReq 495 keep <- rowSums(data$counts) >= opt$cntReq
475 } else if (filt_smpcount) { 496 } else if (filt_smpcount) {
476 keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq 497 keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq
477 } else if (filt_cpm) { 498 } else if (filt_cpm) {
511 png(den_png, width = 1000, height = 500) 532 png(den_png, width = 1000, height = 500)
512 par(mfrow = c(1, 2), cex.axis = 0.8) 533 par(mfrow = c(1, 2), cex.axis = 0.8)
513 534
514 # before filtering 535 # before filtering
515 lcpm1 <- cpm(counts, log = TRUE) 536 lcpm1 <- cpm(counts, log = TRUE)
516 plot(density(lcpm1[, 1]), col = col.group[1], lwd = 2, las = 2, main = "", xlab = "") 537 plot(density(lcpm1[, 1]), col = col_group[1], lwd = 2, las = 2, main = "", xlab = "")
517 title(main = "Density Plot: Raw counts", xlab = "Log-cpm") 538 title(main = "Density Plot: Raw counts", xlab = "Log-cpm")
518 for (i in 2:nsamples) { 539 for (i in 2:nsamples) {
519 den <- density(lcpm1[, i]) 540 den <- density(lcpm1[, i])
520 lines(den$x, den$y, col = col.group[i], lwd = 2) 541 lines(den$x, den$y, col = col_group[i], lwd = 2)
521 } 542 }
522 543
523 # after filtering 544 # after filtering
524 lcpm2 <- cpm(data$counts, log = TRUE) 545 lcpm2 <- cpm(data$counts, log = TRUE)
525 plot(density(lcpm2[, 1]), col = col.group[1], lwd = 2, las = 2, main = "", xlab = "") 546 plot(density(lcpm2[, 1]), col = col_group[1], lwd = 2, las = 2, main = "", xlab = "")
526 title(main = "Density Plot: Filtered counts", xlab = "Log-cpm") 547 title(main = "Density Plot: Filtered counts", xlab = "Log-cpm")
527 for (i in 2:nsamples) { 548 for (i in 2:nsamples) {
528 den <- density(lcpm2[, i]) 549 den <- density(lcpm2[, i])
529 lines(den$x, den$y, col = col.group[i], lwd = 2) 550 lines(den$x, den$y, col = col_group[i], lwd = 2)
530 } 551 }
531 legend("topright", samplenames, text.col = col.group, bty = "n") 552 legend("topright", samplenames, text.col = col_group, bty = "n")
532 img_name <- "Densityplots.png" 553 img_name <- "Densityplots.png"
533 img_addr <- "densityplots.png" 554 img_addr <- "densityplots.png"
534 image_data <- rbind(image_data, data.frame(Label = img_name, Link = img_addr, stringsAsFactors = FALSE)) 555 image_data <- rbind(image_data, data.frame(Label = img_name, Link = img_addr, stringsAsFactors = FALSE))
535 invisible(dev.off()) 556 invisible(dev.off())
536 557
537 # PDF 558 # PDF
538 pdf(den_pdf, width = 14) 559 pdf(den_pdf, width = 14)
539 par(mfrow = c(1, 2), cex.axis = 0.8) 560 par(mfrow = c(1, 2), cex.axis = 0.8)
540 plot(density(lcpm1[, 1]), col = col.group[1], lwd = 2, las = 2, main = "", xlab = "") 561 plot(density(lcpm1[, 1]), col = col_group[1], lwd = 2, las = 2, main = "", xlab = "")
541 title(main = "Density Plot: Raw counts", xlab = "Log-cpm") 562 title(main = "Density Plot: Raw counts", xlab = "Log-cpm")
542 for (i in 2:nsamples) { 563 for (i in 2:nsamples) {
543 den <- density(lcpm1[, i]) 564 den <- density(lcpm1[, i])
544 lines(den$x, den$y, col = col.group[i], lwd = 2) 565 lines(den$x, den$y, col = col_group[i], lwd = 2)
545 } 566 }
546 plot(density(lcpm2[, 1]), col = col.group[1], lwd = 2, las = 2, main = "", xlab = "") 567 plot(density(lcpm2[, 1]), col = col_group[1], lwd = 2, las = 2, main = "", xlab = "")
547 title(main = "Density Plot: Filtered counts", xlab = "Log-cpm") 568 title(main = "Density Plot: Filtered counts", xlab = "Log-cpm")
548 for (i in 2:nsamples) { 569 for (i in 2:nsamples) {
549 den <- density(lcpm2[, i]) 570 den <- density(lcpm2[, i])
550 lines(den$x, den$y, col = col.group[i], lwd = 2) 571 lines(den$x, den$y, col = col_group[i], lwd = 2)
551 } 572 }
552 legend("topright", samplenames, text.col = col.group, bty = "n") 573 legend("topright", samplenames, text.col = col_group, bty = "n")
553 link_name <- "DensityPlots.pdf" 574 link_name <- "DensityPlots.pdf"
554 link_addr <- "densityplots.pdf" 575 link_addr <- "densityplots.pdf"
555 link_data <- rbind(link_data, data.frame(Label = link_name, Link = link_addr, stringsAsFactors = FALSE)) 576 link_data <- rbind(link_data, data.frame(Label = link_name, Link = link_addr, stringsAsFactors = FALSE))
556 invisible(dev.off()) 577 invisible(dev.off())
557 } 578 }
580 colnames(design) <- gsub(factor_list[i], "", colnames(design), fixed = TRUE) 601 colnames(design) <- gsub(factor_list[i], "", colnames(design), fixed = TRUE)
581 } 602 }
582 603
583 # Calculating normalising factors 604 # Calculating normalising factors
584 print("Calculating Normalisation Factors") 605 print("Calculating Normalisation Factors")
585 logcounts <- y #store for plots 606 logcounts <- y # store for plots
586 y <- calcNormFactors(y, method = opt$normOpt) 607 y <- calcNormFactors(y, method = opt$normOpt)
587 608
588 # Generate contrasts information 609 # Generate contrasts information
589 print("Generating Contrasts") 610 print("Generating Contrasts")
590 contrasts <- makeContrasts(contrasts = cons, levels = design) 611 contrasts <- makeContrasts(contrasts = cons, levels = design)
592 ################################################################################ 613 ################################################################################
593 ### Data Output 614 ### Data Output
594 ################################################################################ 615 ################################################################################
595 616
596 # Plot Box plots (before and after normalisation) 617 # Plot Box plots (before and after normalisation)
597 if (opt$normOpt != "none" & "b" %in% plots) { 618 if (opt$normOpt != "none" && "b" %in% plots) {
598 png(box_png, width = 1000, height = 500) 619 png(box_png, width = 1000, height = 500)
599 par(mfrow = c(1, 2), mar = c(6, 4, 2, 2) + 0.1) 620 par(mfrow = c(1, 2), mar = c(6, 4, 2, 2) + 0.1)
600 labels <- colnames(counts) 621 labels <- colnames(counts)
601 622
602 lcpm1 <- cpm(y$counts, log = TRUE) 623 lcpm1 <- cpm(y$counts, log = TRUE)
603 boxplot(lcpm1, las = 2, col = col.group, xaxt = "n", xlab = "") 624 boxplot(lcpm1, las = 2, col = col_group, xaxt = "n", xlab = "")
604 axis(1, at = seq_along(labels), labels = FALSE) 625 axis(1, at = seq_along(labels), labels = FALSE)
605 abline(h = median(lcpm1), col = 4) 626 abline(h = median(lcpm1), col = 4)
606 text(x = seq_along(labels), y = par("usr")[3] - 1, srt = 45, adj = 1, labels = labels, xpd = TRUE) 627 text(x = seq_along(labels), y = par("usr")[3] - 1, srt = 45, adj = 1, labels = labels, xpd = TRUE)
607 title(main = "Box Plot: Unnormalised counts", ylab = "Log-cpm") 628 title(main = "Box Plot: Unnormalised counts", ylab = "Log-cpm")
608 629
609 lcpm2 <- cpm(y, log = TRUE) 630 lcpm2 <- cpm(y, log = TRUE)
610 boxplot(lcpm2, las = 2, col = col.group, xaxt = "n", xlab = "") 631 boxplot(lcpm2, las = 2, col = col_group, xaxt = "n", xlab = "")
611 axis(1, at = seq_along(labels), labels = FALSE) 632 axis(1, at = seq_along(labels), labels = FALSE)
612 text(x = seq_along(labels), y = par("usr")[3] - 1, srt = 45, adj = 1, labels = labels, xpd = TRUE) 633 text(x = seq_along(labels), y = par("usr")[3] - 1, srt = 45, adj = 1, labels = labels, xpd = TRUE)
613 abline(h = median(lcpm2), col = 4) 634 abline(h = median(lcpm2), col = 4)
614 title(main = "Box Plot: Normalised counts", ylab = "Log-cpm") 635 title(main = "Box Plot: Normalised counts", ylab = "Log-cpm")
615 636
618 image_data <- rbind(image_data, data.frame(Label = img_name, Link = img_addr, stringsAsFactors = FALSE)) 639 image_data <- rbind(image_data, data.frame(Label = img_name, Link = img_addr, stringsAsFactors = FALSE))
619 invisible(dev.off()) 640 invisible(dev.off())
620 641
621 pdf(box_pdf, width = 14) 642 pdf(box_pdf, width = 14)
622 par(mfrow = c(1, 2), mar = c(6, 4, 2, 2) + 0.1) 643 par(mfrow = c(1, 2), mar = c(6, 4, 2, 2) + 0.1)
623 boxplot(lcpm1, las = 2, col = col.group, xaxt = "n", xlab = "") 644 boxplot(lcpm1, las = 2, col = col_group, xaxt = "n", xlab = "")
624 axis(1, at = seq_along(labels), labels = FALSE) 645 axis(1, at = seq_along(labels), labels = FALSE)
625 abline(h = median(lcpm1), col = 4) 646 abline(h = median(lcpm1), col = 4)
626 text(x = seq_along(labels), y = par("usr")[3] - 1, srt = 45, adj = 1, labels = labels, xpd = TRUE) 647 text(x = seq_along(labels), y = par("usr")[3] - 1, srt = 45, adj = 1, labels = labels, xpd = TRUE)
627 title(main = "Box Plot: Unnormalised counts", ylab = "Log-cpm") 648 title(main = "Box Plot: Unnormalised counts", ylab = "Log-cpm")
628 boxplot(lcpm2, las = 2, col = col.group, xaxt = "n", xlab = "") 649 boxplot(lcpm2, las = 2, col = col_group, xaxt = "n", xlab = "")
629 axis(1, at = seq_along(labels), labels = FALSE) 650 axis(1, at = seq_along(labels), labels = FALSE)
630 text(x = seq_along(labels), y = par("usr")[3] - 1, srt = 45, adj = 1, labels = labels, xpd = TRUE) 651 text(x = seq_along(labels), y = par("usr")[3] - 1, srt = 45, adj = 1, labels = labels, xpd = TRUE)
631 abline(h = median(lcpm2), col = 4) 652 abline(h = median(lcpm2), col = 4)
632 title(main = "Box Plot: Normalised counts", ylab = "Log-cpm") 653 title(main = "Box Plot: Normalised counts", ylab = "Log-cpm")
633 link_name <- "BoxPlots.pdf" 654 link_name <- "BoxPlots.pdf"
642 663
643 # Scree plot (Variance Explained) code copied from Glimma 664 # Scree plot (Variance Explained) code copied from Glimma
644 665
645 # get column of matrix 666 # get column of matrix
646 get_cols <- function(x, inds) { 667 get_cols <- function(x, inds) {
647 x[, inds, drop = FALSE] 668 x[, inds, drop = FALSE]
648 } 669 }
649 670
650 x <- cpm(y, log = TRUE) 671 x <- cpm(y, log = TRUE)
651 ndim <- nsamples - 1 672 ndim <- nsamples - 1
652 nprobes <- nrow(x) 673 nprobes <- nrow(x)
654 top <- min(top, nprobes) 675 top <- min(top, nprobes)
655 cn <- colnames(x) 676 cn <- colnames(x)
656 bad <- rowSums(is.finite(x)) < nsamples 677 bad <- rowSums(is.finite(x)) < nsamples
657 678
658 if (any(bad)) { 679 if (any(bad)) {
659 warning("Rows containing infinite values have been removed") 680 warning("Rows containing infinite values have been removed")
660 x <- x[!bad, , drop = FALSE] 681 x <- x[!bad, , drop = FALSE]
661 } 682 }
662 683
663 dd <- matrix(0, nrow = nsamples, ncol = nsamples, dimnames = list(cn, cn)) 684 dd <- matrix(0, nrow = nsamples, ncol = nsamples, dimnames = list(cn, cn))
664 topindex <- nprobes - top + 1L 685 topindex <- nprobes - top + 1L
665 for (i in 2L:(nsamples)) { 686 for (i in 2L:(nsamples)) {
666 for (j in 1L:(i - 1L)) { 687 for (j in 1L:(i - 1L)) {
667 dists <- (get_cols(x, i) - get_cols(x, j))^2 688 dists <- (get_cols(x, i) - get_cols(x, j))^2
668 dists <- sort.int(dists, partial = topindex) 689 dists <- sort.int(dists, partial = topindex)
669 topdist <- dists[topindex:nprobes] 690 topdist <- dists[topindex:nprobes]
670 dd[i, j] <- sqrt(mean(topdist)) 691 dd[i, j] <- sqrt(mean(topdist))
671 } 692 }
672 } 693 }
673 694
674 a1 <- suppressWarnings(cmdscale(as.dist(dd), k = min(ndim, 8), eig = TRUE)) 695 a1 <- suppressWarnings(cmdscale(as.dist(dd), k = min(ndim, 8), eig = TRUE))
675 eigen <- data.frame(name = 1:min(ndim, 8), eigen = round(a1$eig[1:min(ndim, 8)] / sum(a1$eig), 2)) 696 eigen <- data.frame(name = 1:min(ndim, 8), eigen = round(a1$eig[1:min(ndim, 8)] / sum(a1$eig), 2))
676 png(mdsscree_png, width = 1000, height = 500) 697 png(mdsscree_png, width = 1000, height = 500)
677 par(mfrow = c(1, 2)) 698 par(mfrow = c(1, 2))
678 plotMDS(y, labels = samplenames, col = as.numeric(factors[, 1]), main = "MDS Plot: Dims 1 and 2") 699 plotMDS(y, labels = samplenames, col = as.numeric(factors[, 1]), main = "MDS Plot: Dims 1 and 2")
679 barplot(eigen$eigen, names.arg = eigen$name, main = "Scree Plot: Variance Explained", xlab = "Dimension", ylab = "Proportion", las = 1) 700 barplot(eigen$eigen, names.arg = eigen$name, main = "Scree Plot: Variance Explained", xlab = "Dimension", ylab = "Proportion", las = 1)
680 img_name <- paste0("MDSPlot_", names(factors)[1], ".png") 701 img_name <- paste0("MDSPlot_", names(factors)[1], ".png")
681 img_addr <- "mdsscree.png" 702 img_addr <- "mdsscree.png"
682 image_data <- rbind(image_data, data.frame(Label = img_name, Link = img_addr, stringsAsFactors = FALSE)) 703 image_data <- rbind(image_data, data.frame(Label = img_name, Link = img_addr, stringsAsFactors = FALSE))
683 invisible(dev.off()) 704 invisible(dev.off())
684 705
685 pdf(mdsscree_pdf, width = 14) 706 pdf(mdsscree_pdf, width = 14)
686 par(mfrow = c(1, 2)) 707 par(mfrow = c(1, 2))
687 plotMDS(y, labels = samplenames, col = as.numeric(factors[, 1]), main = "MDS Plot: Dims 1 and 2") 708 plotMDS(y, labels = samplenames, col = as.numeric(factors[, 1]), main = "MDS Plot: Dims 1 and 2")
688 barplot(eigen$eigen, names.arg = eigen$name, main = "Scree Plot: Variance Explained", xlab = "Dimension", ylab = "Proportion", las = 1) 709 barplot(eigen$eigen, names.arg = eigen$name, main = "Scree Plot: Variance Explained", xlab = "Dimension", ylab = "Proportion", las = 1)
689 link_name <- paste0("MDSPlot_", names(factors)[1], ".pdf") 710 link_name <- paste0("MDSPlot_", names(factors)[1], ".pdf")
690 link_addr <- "mdsscree.pdf" 711 link_addr <- "mdsscree.pdf"
691 link_data <- rbind(link_data, data.frame(Label = link_name, Link = link_addr, stringsAsFactors = FALSE)) 712 link_data <- rbind(link_data, data.frame(Label = link_name, Link = link_addr, stringsAsFactors = FALSE))
692 invisible(dev.off()) 713 invisible(dev.off())
693 714
694 # generate Glimma interactive MDS Plot 715 # generate Glimma interactive MDS Plot
695 if ("i" %in% plots) { 716 if ("i" %in% plots) {
696 Glimma::glMDSPlot(y, labels = samplenames, groups = factors[, 1], 717 Glimma::glMDSPlot(y,
697 folder = "glimma_MDS", launch = FALSE) 718 labels = samplenames, groups = factors[, 1],
719 folder = "glimma_MDS", launch = FALSE
720 )
698 link_name <- "Glimma_MDSPlot.html" 721 link_name <- "Glimma_MDSPlot.html"
699 link_addr <- "glimma_MDS/MDS-Plot.html" 722 link_addr <- "glimma_MDS/MDS-Plot.html"
700 link_data <- rbind(link_data, c(link_name, link_addr)) 723 link_data <- rbind(link_data, c(link_name, link_addr))
701 } 724 }
702 725
787 invisible(dev.off()) 810 invisible(dev.off())
788 811
789 # Generating fit data and top table with weights 812 # Generating fit data and top table with weights
790 wts <- vdata$weights 813 wts <- vdata$weights
791 voomfit <- lmFit(vdata, design, weights = wts) 814 voomfit <- lmFit(vdata, design, weights = wts)
792
793 } else { 815 } else {
794 voom_pdf <- make_out("voomplot.pdf") 816 voom_pdf <- make_out("voomplot.pdf")
795 voom_png <- make_out("voomplot.png") 817 voom_png <- make_out("voomplot.png")
796 # Creating voom data object and plot 818 # Creating voom data object and plot
797 png(voom_png, width = 500, height = 500) 819 png(voom_png, width = 500, height = 500)
810 832
811 # Generate voom fit 833 # Generate voom fit
812 voomfit <- lmFit(vdata, design) 834 voomfit <- lmFit(vdata, design)
813 } 835 }
814 836
815 # Save normalised counts (log2cpm) 837 # Save normalised counts (log2cpm)
816 if (want_norm) { 838 if (want_norm) {
817 norm_counts <- data.frame(vdata$genes, vdata$E, check.names = FALSE) 839 norm_counts <- data.frame(vdata$genes, vdata$E, check.names = FALSE)
818 write.table(norm_counts, file = norm_out, row.names = FALSE, sep = "\t", quote = FALSE) 840 write.table(norm_counts, file = norm_out, row.names = FALSE, sep = "\t", quote = FALSE)
819 link_data <- rbind(link_data, c((paste0(de_method, "_", "normcounts.tsv")), (paste0(de_method, "_", "normcounts")))) 841 link_data <- rbind(link_data, c((paste0(de_method, "_", "normcounts.tsv")), (paste0(de_method, "_", "normcounts"))))
820 } 842 }
845 link_name <- "SAPlot.pdf" 867 link_name <- "SAPlot.pdf"
846 link_addr <- "saplot.pdf" 868 link_addr <- "saplot.pdf"
847 link_data <- rbind(link_data, c(link_name, link_addr)) 869 link_data <- rbind(link_data, c(link_name, link_addr))
848 invisible(dev.off()) 870 invisible(dev.off())
849 871
850 # Save library size info 872 # Save library size info
851 if (want_libinfo) { 873 if (want_libinfo) {
852 efflibsize <- round(y$samples$lib.size * y$samples$norm.factors) 874 efflibsize <- round(y$samples$lib.size * y$samples$norm.factors)
853 libsizeinfo <- cbind(y$samples, EffectiveLibrarySize = efflibsize) 875 libsizeinfo <- cbind(y$samples, EffectiveLibrarySize = efflibsize)
854 libsizeinfo$lib.size <- round(libsizeinfo$lib.size) # nolint 876 libsizeinfo$lib.size <- round(libsizeinfo$lib.size) # nolint
855 names(libsizeinfo)[names(libsizeinfo) == "sampleID"] <- "SampleID" 877 names(libsizeinfo)[names(libsizeinfo) == "sampleID"] <- "SampleID"
867 } else { 889 } else {
868 fit <- treat(fit, lfc = opt$lfcReq, robust = FALSE) 890 fit <- treat(fit, lfc = opt$lfcReq, robust = FALSE)
869 } 891 }
870 } 892 }
871 893
872 status <- decideTests(fit, adjust.method = opt$pAdjOpt, p.value = opt$pValReq, 894 status <- decideTests(fit,
873 lfc = opt$lfcReq) 895 adjust.method = opt$pAdjOpt, p.value = opt$pValReq,
896 lfc = opt$lfcReq
897 )
874 sum_status <- summary(status) 898 sum_status <- summary(status)
875 899
876 for (i in seq_along(cons)) { 900 for (i in seq_along(cons)) {
877 con_name <- cons[i] 901 con_name <- cons[i]
878 con <- cons[i] 902 con <- cons[i]
883 flat_count[i] <- sum_status["NotSig", i] 907 flat_count[i] <- sum_status["NotSig", i]
884 908
885 # Write top expressions table 909 # Write top expressions table
886 if (want_treat) { 910 if (want_treat) {
887 top <- topTreat(fit, coef = i, adjust.method = opt$pAdjOpt, number = Inf, sort.by = "P") 911 top <- topTreat(fit, coef = i, adjust.method = opt$pAdjOpt, number = Inf, sort.by = "P")
888 } else{ 912 } else {
889 top <- topTable(fit, coef = i, adjust.method = opt$pAdjOpt, number = Inf, sort.by = "P") 913 top <- topTable(fit, coef = i, adjust.method = opt$pAdjOpt, number = Inf, sort.by = "P")
890 } 914 }
891 write.table(top, file = top_out[i], row.names = FALSE, sep = "\t", quote = FALSE) 915 write.table(top, file = top_out[i], row.names = FALSE, sep = "\t", quote = FALSE)
892 link_name <- paste0(de_method, "_", con, ".tsv") 916 link_name <- paste0(de_method, "_", con, ".tsv")
893 link_addr <- paste0(de_method, "_", con, ".tsv") 917 link_addr <- paste0(de_method, "_", con, ".tsv")
894 link_data <- rbind(link_data, c(link_name, link_addr)) 918 link_data <- rbind(link_data, c(link_name, link_addr))
895 919
896 # Plot MD (log ratios vs mean average) using limma package on weighted 920 # Plot MD (log ratios vs mean average) using limma package on weighted
897 pdf(md_pdf[i]) 921 pdf(md_pdf[i])
898 limma::plotMD(fit, status = status[, i], coef = i, 922 limma::plotMD(fit,
923 status = status[, i], coef = i,
899 main = paste("MD Plot:", unmake_names(con)), 924 main = paste("MD Plot:", unmake_names(con)),
900 hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1), 925 hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1),
901 xlab = "Average Expression", ylab = "logFC") 926 xlab = "Average Expression", ylab = "logFC"
927 )
902 abline(h = 0, col = "grey", lty = 2) 928 abline(h = 0, col = "grey", lty = 2)
903 link_name <- paste0("MDPlot_", con, ".pdf") 929 link_name <- paste0("MDPlot_", con, ".pdf")
904 link_addr <- paste0("mdplot_", con, ".pdf") 930 link_addr <- paste0("mdplot_", con, ".pdf")
905 link_data <- rbind(link_data, c(link_name, link_addr)) 931 link_data <- rbind(link_data, c(link_name, link_addr))
906 invisible(dev.off()) 932 invisible(dev.off())
907 933
908 # Generate Glimma interactive Volcano, MD plot and tables, requires annotation file (assumes gene labels/symbols in 2nd column) 934 # Generate Glimma interactive Volcano, MD plot and tables, requires annotation file (assumes gene labels/symbols in 2nd column)
909 if ("i" %in% plots & have_anno) { 935 if ("i" %in% plots && have_anno) {
910 # make gene labels unique to handle NAs 936 # make gene labels unique to handle NAs
911 geneanno <- y$genes 937 geneanno <- y$genes
912 geneanno[, 2] <- make.unique(geneanno[, 2]) 938 geneanno[, 2] <- make.unique(geneanno[, 2])
913 939
914 # use the logcpms for the counts 940 # use the logcpms for the counts
915 if (want_trend) { 941 if (want_trend) {
916 cnts <- logcpm 942 cnts <- logcpm
917 } else{ 943 } else {
918 cnts <- vdata$E 944 cnts <- vdata$E
919 } 945 }
920 946
921 # MD plot 947 # MD plot
922 Glimma::glMDPlot(fit, coef = i, counts = cnts, anno = geneanno, groups = factors[, 1], 948 Glimma::glMDPlot(fit,
923 status = status[, i], sample.cols = col.group, 949 coef = i, counts = cnts, anno = geneanno, groups = factors[, 1],
924 main = paste("MD Plot:", unmake_names(con)), side.main = colnames(y$genes)[2], 950 status = status[, i], sample.cols = col_group,
925 folder = paste0("glimma_", unmake_names(con)), launch = FALSE) 951 main = paste("MD Plot:", unmake_names(con)), side.main = colnames(y$genes)[2],
952 folder = paste0("glimma_", unmake_names(con)), launch = FALSE
953 )
926 link_name <- paste0("Glimma_MDPlot_", con, ".html") 954 link_name <- paste0("Glimma_MDPlot_", con, ".html")
927 link_addr <- paste0("glimma_", con, "/MD-Plot.html") 955 link_addr <- paste0("glimma_", con, "/MD-Plot.html")
928 link_data <- rbind(link_data, c(link_name, link_addr)) 956 link_data <- rbind(link_data, c(link_name, link_addr))
929 957
930 # Volcano plot 958 # Volcano plot
931 Glimma::glXYPlot(x = fit$coefficients[, i], y = -log10(fit$p.value[, i]), counts = cnts, anno = geneanno, groups = factors[, 1], 959 Glimma::glXYPlot(
932 status = status[, i], sample.cols = col.group, 960 x = fit$coefficients[, i], y = -log10(fit$p.value[, i]), counts = cnts, anno = geneanno, groups = factors[, 1],
961 status = status[, i], sample.cols = col_group,
933 main = paste("Volcano Plot:", unmake_names(con)), side.main = colnames(y$genes)[2], 962 main = paste("Volcano Plot:", unmake_names(con)), side.main = colnames(y$genes)[2],
934 xlab = "logFC", ylab = "-log10(P-value)", 963 xlab = "logFC", ylab = "-log10(P-value)",
935 folder = paste0("glimma_volcano_", unmake_names(con)), launch = FALSE) 964 folder = paste0("glimma_volcano_", unmake_names(con)), launch = FALSE
965 )
936 link_name <- paste0("Glimma_VolcanoPlot_", con, ".html") 966 link_name <- paste0("Glimma_VolcanoPlot_", con, ".html")
937 link_addr <- paste0("glimma_volcano_", con, "/XY-Plot.html") 967 link_addr <- paste0("glimma_volcano_", con, "/XY-Plot.html")
938 link_data <- rbind(link_data, c(link_name, link_addr)) 968 link_data <- rbind(link_data, c(link_name, link_addr))
939 } 969 }
940 970
944 # labels must be in second column currently 974 # labels must be in second column currently
945 labels <- fit$genes[, 2] 975 labels <- fit$genes[, 2]
946 } else { 976 } else {
947 labels <- fit$genes$GeneID 977 labels <- fit$genes$GeneID
948 } 978 }
949 limma::volcanoplot(fit, coef = i, 979 limma::volcanoplot(fit,
980 coef = i,
950 main = paste("Volcano Plot:", unmake_names(con)), 981 main = paste("Volcano Plot:", unmake_names(con)),
951 highlight = opt$topgenes, 982 highlight = opt$topgenes,
952 names = labels) 983 names = labels
984 )
953 link_name <- paste0("VolcanoPlot_", con, ".pdf") 985 link_name <- paste0("VolcanoPlot_", con, ".pdf")
954 link_addr <- paste0("volplot_", con, ".pdf") 986 link_addr <- paste0("volplot_", con, ".pdf")
955 link_data <- rbind(link_data, c(link_name, link_addr)) 987 link_data <- rbind(link_data, c(link_name, link_addr))
956 invisible(dev.off()) 988 invisible(dev.off())
957 989
958 # PNG of MD and Volcano 990 # PNG of MD and Volcano
959 png(mdvol_png[i], width = 1000, height = 500) 991 png(mdvol_png[i], width = 1000, height = 500)
960 par(mfrow = c(1, 2), mar = c(5, 4, 2, 2) + 0.1, oma = c(0, 0, 3, 0)) 992 par(mfrow = c(1, 2), mar = c(5, 4, 2, 2) + 0.1, oma = c(0, 0, 3, 0))
961 993
962 # MD plot 994 # MD plot
963 limma::plotMD(fit, status = status[, i], coef = i, main = "MD Plot", 995 limma::plotMD(fit,
996 status = status[, i], coef = i, main = "MD Plot",
964 hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1), 997 hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1),
965 xlab = "Average Expression", ylab = "logFC") 998 xlab = "Average Expression", ylab = "logFC"
999 )
966 abline(h = 0, col = "grey", lty = 2) 1000 abline(h = 0, col = "grey", lty = 2)
967 1001
968 # Volcano 1002 # Volcano
969 if (have_anno) { 1003 if (have_anno) {
970 # labels must be in second column currently 1004 # labels must be in second column currently
971 limma::volcanoplot(fit, coef = i, main = "Volcano Plot", 1005 limma::volcanoplot(fit,
1006 coef = i, main = "Volcano Plot",
972 highlight = opt$topgenes, 1007 highlight = opt$topgenes,
973 names = fit$genes[, 2]) 1008 names = fit$genes[, 2]
1009 )
974 } else { 1010 } else {
975 limma::volcanoplot(fit, coef = i, main = "Volcano Plot", 1011 limma::volcanoplot(fit,
1012 coef = i, main = "Volcano Plot",
976 highlight = opt$topgenes, 1013 highlight = opt$topgenes,
977 names = fit$genes$GeneID) 1014 names = fit$genes$GeneID
1015 )
978 } 1016 }
979 1017
980 img_name <- paste0("MDVolPlot_", con) 1018 img_name <- paste0("MDVolPlot_", con)
981 img_addr <- paste0("mdvolplot_", con, ".png") 1019 img_addr <- paste0("mdvolplot_", con, ".png")
982 image_data <- rbind(image_data, c(img_name, img_addr)) 1020 image_data <- rbind(image_data, c(img_name, img_addr))
997 # labels must be in second column currently 1035 # labels must be in second column currently
998 labels <- top[topgenes, 2] 1036 labels <- top[topgenes, 2]
999 } else { 1037 } else {
1000 labels <- rownames(topexp) 1038 labels <- rownames(topexp)
1001 } 1039 }
1002 heatmap.2(topexp, scale = "row", Colv = FALSE, Rowv = FALSE, dendrogram = "none", 1040 heatmap.2(topexp,
1041 scale = "row", Colv = FALSE, Rowv = FALSE, dendrogram = "none",
1003 main = paste("Contrast:", unmake_names(con), "\nTop", opt$topgenes, "genes by adj.P.Val"), 1042 main = paste("Contrast:", unmake_names(con), "\nTop", opt$topgenes, "genes by adj.P.Val"),
1004 trace = "none", density.info = "none", lhei = c(2, 10), margin = c(8, 6), labRow = labels, cexRow = 0.7, srtCol = 45, 1043 trace = "none", density.info = "none", lhei = c(2, 10), margin = c(8, 6), labRow = labels, cexRow = 0.7, srtCol = 45,
1005 col = mycol, ColSideColors = col.group) 1044 col = mycol, ColSideColors = col_group
1045 )
1006 link_name <- paste0("Heatmap_", con, ".pdf") 1046 link_name <- paste0("Heatmap_", con, ".pdf")
1007 link_addr <- paste0("heatmap_", con, ".pdf") 1047 link_addr <- paste0("heatmap_", con, ".pdf")
1008 link_data <- rbind(link_data, c(link_name, link_addr)) 1048 link_data <- rbind(link_data, c(link_name, link_addr))
1009 invisible(dev.off()) 1049 invisible(dev.off())
1010 } 1050 }
1011 1051
1012 if ("s" %in% plots) { 1052 if ("s" %in% plots) {
1013 # Plot Stripcharts of top genes 1053 # Plot Stripcharts of top genes
1014 pdf(strip_pdf[i], title = paste("Contrast:", unmake_names(con))) 1054 pdf(strip_pdf[i], title = paste("Contrast:", unmake_names(con)))
1015 par(mfrow = c(3, 2), cex.main = 0.8, cex.axis = 0.8) 1055 par(mfrow = c(3, 2), cex.main = 0.8, cex.axis = 0.8)
1016 cols <- unique(col.group) 1056 cols <- unique(col_group)
1017 1057
1018 for (j in seq_along(topgenes)) { 1058 for (j in seq_along(topgenes)) {
1019 lfc <- round(top[topgenes[j], "logFC"], 2) 1059 lfc <- round(top[topgenes[j], "logFC"], 2)
1020 pval <- round(top[topgenes[j], "adj.P.Val"], 5) 1060 pval <- round(top[topgenes[j], "adj.P.Val"], 5)
1021 if (want_trend) { 1061 if (want_trend) {
1022 stripchart(plot_data[topgenes[j], ] ~ factors[, 1], vertical = TRUE, las = 2, pch = 16, cex = 0.8, cex.lab = 0.8, col = cols, 1062 stripchart(plot_data[topgenes[j], ] ~ factors[, 1],
1023 method = "jitter", ylab = "Normalised log2 expression", main = paste0(labels[j], "\nlogFC=", lfc, ", adj.P.Val=", pval)) 1063 vertical = TRUE, las = 2, pch = 16, cex = 0.8, cex.lab = 0.8, col = cols,
1064 method = "jitter", ylab = "Normalised log2 expression", main = paste0(labels[j], "\nlogFC=", lfc, ", adj.P.Val=", pval)
1065 )
1024 } else { 1066 } else {
1025 stripchart(plot_data$E[topgenes[j], ] ~ factors[, 1], vertical = TRUE, las = 2, pch = 16, cex = 0.8, cex.lab = 0.8, col = cols, 1067 stripchart(plot_data$E[topgenes[j], ] ~ factors[, 1],
1026 method = "jitter", ylab = "Normalised log2 expression", main = paste0(labels[j], "\nlogFC=", lfc, ", adj.P.Val=", pval)) 1068 vertical = TRUE, las = 2, pch = 16, cex = 0.8, cex.lab = 0.8, col = cols,
1069 method = "jitter", ylab = "Normalised log2 expression", main = paste0(labels[j], "\nlogFC=", lfc, ", adj.P.Val=", pval)
1070 )
1027 } 1071 }
1028 } 1072 }
1029 link_name <- paste0("Stripcharts_", con, ".pdf") 1073 link_name <- paste0("Stripcharts_", con, ".pdf")
1030 link_addr <- paste0("stripcharts_", con, ".pdf") 1074 link_addr <- paste0("stripcharts_", con, ".pdf")
1031 link_data <- rbind(link_data, c(link_name, link_addr)) 1075 link_data <- rbind(link_data, c(link_name, link_addr))
1037 1081
1038 # Save relevant items as rda object 1082 # Save relevant items as rda object
1039 if (want_rda) { 1083 if (want_rda) {
1040 print("Saving RData") 1084 print("Saving RData")
1041 if (want_weight) { 1085 if (want_weight) {
1042 save(counts, data, y, status, plot_data, labels, factors, wts, fit, top, contrast_data, contrasts, design, 1086 save(counts, data, y, status, plot_data, labels, factors, wts, fit, top, contrast_data, contrasts, design,
1043 file = rda_out, ascii = TRUE) 1087 file = rda_out, ascii = TRUE
1088 )
1044 } else { 1089 } else {
1045 save(counts, data, y, status, plot_data, labels, factors, fit, top, contrast_data, contrasts, design, 1090 save(counts, data, y, status, plot_data, labels, factors, fit, top, contrast_data, contrasts, design,
1046 file = rda_out, ascii = TRUE) 1091 file = rda_out, ascii = TRUE
1092 )
1047 } 1093 }
1048 link_data <- rbind(link_data, c((paste0(de_method, "_analysis.RData")), (paste0(de_method, "_analysis.RData")))) 1094 link_data <- rbind(link_data, c((paste0(de_method, "_analysis.RData")), (paste0(de_method, "_analysis.RData"))))
1049 } 1095 }
1050 1096
1051 # Record session info 1097 # Record session info
1052 writeLines(capture.output(sessionInfo()), session_out) 1098 writeLines(capture.output(sessionInfo()), session_out)
1053 link_data <- rbind(link_data, c("Session Info", "session_info.txt")) 1099 link_data <- rbind(link_data, c("Session Info", "session_info.txt"))
1054 1100
1055 # Record ending time and calculate total run time 1101 # Record ending time and calculate total run time
1056 time_end <- as.character(Sys.time()) 1102 time_end <- Sys.time()
1057 time_taken <- capture.output(round(difftime(time_end, time_start), digits = 3)) 1103 time_taken <- capture.output(round(difftime(time_end, time_start), digits = 2))
1058 time_taken <- gsub("Time difference of ", "", time_taken, fixed = TRUE) 1104 time_taken <- gsub("Time difference of ", "", time_taken, fixed = TRUE)
1105 time_start <- format(time_start, "%A, %B %d, %Y %H:%M:%S")
1106 time_end <- format(time_end, "%A, %B %d, %Y %H:%M:%S")
1059 ################################################################################ 1107 ################################################################################
1060 ### HTML Generation 1108 ### HTML Generation
1061 ################################################################################ 1109 ################################################################################
1062 1110
1063 # Clear file 1111 # Clear file
1095 cata("</tr>\n") 1143 cata("</tr>\n")
1096 } 1144 }
1097 cata("</table>") 1145 cata("</table>")
1098 1146
1099 cata("<h4>Plots:</h4>\n") 1147 cata("<h4>Plots:</h4>\n")
1100 #PDFs 1148 # PDFs
1101 for (i in seq_len(nrow(link_data))) { 1149 for (i in seq_len(nrow(link_data))) {
1102 if (grepl(".pdf", link_data$Link[i]) & grepl("density|cpm|boxplot|mds|mdplots|voom|saplot", link_data$Link[i])) { 1150 if (grepl(".pdf", link_data$Link[i]) && grepl("density|cpm|boxplot|mds|mdplots|voom|saplot", link_data$Link[i])) {
1103 html_link(link_data$Link[i], link_data$Label[i]) 1151 html_link(link_data$Link[i], link_data$Label[i])
1104 } 1152 }
1105 } 1153 }
1106 1154
1107 for (i in seq_len(nrow(link_data))) { 1155 for (i in seq_len(nrow(link_data))) {
1108 if (grepl("mdplot_", link_data$Link[i])) { 1156 if (grepl("mdplot_", link_data$Link[i])) {
1109 html_link(link_data$Link[i], link_data$Label[i]) 1157 html_link(link_data$Link[i], link_data$Label[i])
1110 } 1158 }
1111 } 1159 }
1112 1160
1113 for (i in seq_len(nrow(link_data))) { 1161 for (i in seq_len(nrow(link_data))) {
1114 if (grepl("volplot", link_data$Link[i])) { 1162 if (grepl("volplot", link_data$Link[i])) {
1115 html_link(link_data$Link[i], link_data$Label[i]) 1163 html_link(link_data$Link[i], link_data$Label[i])
1116 } 1164 }
1117 } 1165 }
1118 1166
1119 for (i in seq_len(nrow(link_data))) { 1167 for (i in seq_len(nrow(link_data))) {
1120 if (grepl("heatmap", link_data$Link[i])) { 1168 if (grepl("heatmap", link_data$Link[i])) {
1121 html_link(link_data$Link[i], link_data$Label[i]) 1169 html_link(link_data$Link[i], link_data$Label[i])
1122 } 1170 }
1123 } 1171 }
1124 1172
1125 for (i in seq_len(nrow(link_data))) { 1173 for (i in seq_len(nrow(link_data))) {
1126 if (grepl("stripcharts", link_data$Link[i])) { 1174 if (grepl("stripcharts", link_data$Link[i])) {
1127 html_link(link_data$Link[i], link_data$Label[i]) 1175 html_link(link_data$Link[i], link_data$Label[i])
1128 } 1176 }
1129 } 1177 }
1130 1178
1131 cata("<h4>Tables:</h4>\n") 1179 cata("<h4>Tables:</h4>\n")
1132 for (i in seq_len(nrow(link_data))) { 1180 for (i in seq_len(nrow(link_data))) {
1133 if (grepl("counts$", link_data$Link[i])) { 1181 if (grepl("counts$", link_data$Link[i])) {
1146 } 1194 }
1147 } 1195 }
1148 1196
1149 if ("i" %in% plots) { 1197 if ("i" %in% plots) {
1150 cata("<h4>Glimma Interactive Results:</h4>\n") 1198 cata("<h4>Glimma Interactive Results:</h4>\n")
1151 for (i in seq_len(nrow(link_data))) { 1199 for (i in seq_len(nrow(link_data))) {
1152 if (grepl("glimma", link_data$Link[i])) { 1200 if (grepl("glimma", link_data$Link[i])) {
1153 html_link(link_data$Link[i], link_data$Label[i]) 1201 html_link(link_data$Link[i], link_data$Label[i])
1154 } 1202 }
1155 } 1203 }
1156 } 1204 }
1157 1205
1158 cata("<p>Alt-click links to download file.</p>\n") 1206 cata("<p>Alt-click links to download file.</p>\n")
1159 cata("<p>Click floppy disc icon associated history item to download ") 1207 cata("<p>Click floppy disc icon associated history item to download ")
1160 cata("all files.</p>\n") 1208 cata("all files.</p>\n")
1163 cata("<h4>Additional Information</h4>\n") 1211 cata("<h4>Additional Information</h4>\n")
1164 cata("<ul>\n") 1212 cata("<ul>\n")
1165 1213
1166 if (filt_cpm || filt_smpcount || filt_totcount) { 1214 if (filt_cpm || filt_smpcount || filt_totcount) {
1167 if (filt_cpm) { 1215 if (filt_cpm) {
1168 temp_str <- paste("Genes without more than", opt$cpmReq, 1216 temp_str <- paste(
1169 "CPM in at least", opt$sampleReq, "samples are insignificant", 1217 "Genes without more than", opt$cpmReq,
1170 "and filtered out.") 1218 "CPM in at least", opt$sampleReq, "samples are insignificant",
1219 "and filtered out."
1220 )
1171 } else if (filt_smpcount) { 1221 } else if (filt_smpcount) {
1172 temp_str <- paste("Genes without more than", opt$cntReq, 1222 temp_str <- paste(
1173 "counts in at least", opt$sampleReq, "samples are insignificant", 1223 "Genes without more than", opt$cntReq,
1174 "and filtered out.") 1224 "counts in at least", opt$sampleReq, "samples are insignificant",
1225 "and filtered out."
1226 )
1175 } else if (filt_totcount) { 1227 } else if (filt_totcount) {
1176 temp_str <- paste("Genes without more than", opt$cntReq, 1228 temp_str <- paste(
1177 "counts, after summing counts for all samples, are insignificant", 1229 "Genes without more than", opt$cntReq,
1178 "and filtered out.") 1230 "counts, after summing counts for all samples, are insignificant",
1231 "and filtered out."
1232 )
1179 } 1233 }
1180 1234
1181 list_item(temp_str) 1235 list_item(temp_str)
1182 filter_prop <- round(filtered_count / prefilter_count * 100, digits = 2) 1236 filter_prop <- round(filtered_count / prefilter_count * 100, digits = 2)
1183 temp_str <- paste0(filtered_count, " of ", prefilter_count, " (", filter_prop, 1237 temp_str <- paste0(
1184 "%) genes were filtered out for low expression.") 1238 filtered_count, " of ", prefilter_count, " (", filter_prop,
1239 "%) genes were filtered out for low expression."
1240 )
1185 list_item(temp_str) 1241 list_item(temp_str)
1186 } 1242 }
1187 list_item(opt$normOpt, " was the method used to normalise library sizes.") 1243 list_item(opt$normOpt, " was the method used to normalise library sizes.")
1188 if (want_trend) { 1244 if (want_trend) {
1189 list_item("The limma-trend method was used.") 1245 list_item("The limma-trend method was used.")
1205 list_item("eBayes was used with robust settings (robust = TRUE).") 1261 list_item("eBayes was used with robust settings (robust = TRUE).")
1206 } 1262 }
1207 } 1263 }
1208 if (opt$pAdjOpt != "none") { 1264 if (opt$pAdjOpt != "none") {
1209 if (opt$pAdjOpt == "BH" || opt$pAdjOpt == "BY") { 1265 if (opt$pAdjOpt == "BH" || opt$pAdjOpt == "BY") {
1210 temp_str <- paste0("MD Plot highlighted genes are significant at FDR ", 1266 temp_str <- paste0(
1211 "of ", opt$pValReq, " and exhibit log2-fold-change of at ", 1267 "MD Plot highlighted genes are significant at FDR ",
1212 "least ", opt$lfcReq, ".") 1268 "of ", opt$pValReq, " and exhibit log2-fold-change of at ",
1269 "least ", opt$lfcReq, "."
1270 )
1213 list_item(temp_str) 1271 list_item(temp_str)
1214 } else if (opt$pAdjOpt == "holm") { 1272 } else if (opt$pAdjOpt == "holm") {
1215 temp_str <- paste0("MD Plot highlighted genes are significant at adjusted ", 1273 temp_str <- paste0(
1216 "p-value of ", opt$pValReq, " by the Holm(1979) ", 1274 "MD Plot highlighted genes are significant at adjusted ",
1217 "method, and exhibit log2-fold-change of at least ", 1275 "p-value of ", opt$pValReq, " by the Holm(1979) ",
1218 opt$lfcReq, ".") 1276 "method, and exhibit log2-fold-change of at least ",
1277 opt$lfcReq, "."
1278 )
1219 list_item(temp_str) 1279 list_item(temp_str)
1220 } 1280 }
1221 } else { 1281 } else {
1222 temp_str <- paste0("MD Plot highlighted genes are significant at p-value ", 1282 temp_str <- paste0(
1223 "of ", opt$pValReq, " and exhibit log2-fold-change of at ", 1283 "MD Plot highlighted genes are significant at p-value ",
1224 "least ", opt$lfcReq, ".") 1284 "of ", opt$pValReq, " and exhibit log2-fold-change of at ",
1225 list_item(temp_str) 1285 "least ", opt$lfcReq, "."
1286 )
1287 list_item(temp_str)
1226 } 1288 }
1227 cata("</ul>\n") 1289 cata("</ul>\n")
1228 1290
1229 cata("<h4>Summary of experimental data:</h4>\n") 1291 cata("<h4>Summary of experimental data:</h4>\n")
1230 1292
1245 for (i in seq_len(nrow(factors))) { 1307 for (i in seq_len(nrow(factors))) {
1246 cata("<tr>\n") 1308 cata("<tr>\n")
1247 table_head_item(row.names(factors)[i]) 1309 table_head_item(row.names(factors)[i])
1248 for (j in seq_len(ncol(factors))) { 1310 for (j in seq_len(ncol(factors))) {
1249 table_item(as.character(unmake_names(factors[i, j]))) 1311 table_item(as.character(unmake_names(factors[i, j])))
1250 } 1312 }
1251 cata("</tr>\n") 1313 cata("</tr>\n")
1252 } 1314 }
1253 cata("</table>") 1315 cata("</table>")
1254 1316
1255 cit <- character() 1317 cit <- character()
1256 link <- character() 1318 link <- character()
1257 link[1] <- paste0("<a href=\"", 1319 link[1] <- paste0(
1258 "http://www.bioconductor.org/packages/release/bioc/", 1320 "<a href=\"",
1259 "vignettes/limma/inst/doc/usersguide.pdf", 1321 "http://www.bioconductor.org/packages/release/bioc/",
1260 "\">", "limma User's Guide", "</a>.") 1322 "vignettes/limma/inst/doc/usersguide.pdf",
1261 1323 "\">", "limma User's Guide", "</a>."
1262 link[2] <- paste0("<a href=\"", 1324 )
1263 "http://www.bioconductor.org/packages/release/bioc/", 1325
1264 "vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf", 1326 link[2] <- paste0(
1265 "\">", "edgeR User's Guide", "</a>") 1327 "<a href=\"",
1328 "http://www.bioconductor.org/packages/release/bioc/",
1329 "vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf",
1330 "\">", "edgeR User's Guide", "</a>"
1331 )
1266 1332
1267 cit[1] <- paste("Please cite the following paper for this tool:") 1333 cit[1] <- paste("Please cite the following paper for this tool:")
1268 1334
1269 cit[2] <- paste("Liu R, Holik AZ, Su S, Jansz N, Chen K, Leong HS, Blewitt ME,", 1335 cit[2] <- paste(
1270 "Asselin-Labat ML, Smyth GK, Ritchie ME (2015). Why weight? ", 1336 "Liu R, Holik AZ, Su S, Jansz N, Chen K, Leong HS, Blewitt ME,",
1271 "Modelling sample and observational level variability improves power ", 1337 "Asselin-Labat ML, Smyth GK, Ritchie ME (2015). Why weight? ",
1272 "in RNA-seq analyses. Nucleic Acids Research, 43(15), e97.") 1338 "Modelling sample and observational level variability improves power ",
1273 1339 "in RNA-seq analyses. Nucleic Acids Research, 43(15), e97."
1274 cit[3] <- paste("Please cite the paper below for the limma software itself.", 1340 )
1275 "Please also try to cite the appropriate methodology articles", 1341 cit[3] <- paste(
1276 "that describe the statistical methods implemented in limma,", 1342 "Please cite the paper below for the limma software itself.",
1277 "depending on which limma functions you are using. The", 1343 "Please also try to cite the appropriate methodology articles",
1278 "methodology articles are listed in Section 2.1 of the", 1344 "that describe the statistical methods implemented in limma,",
1279 link[1], 1345 "depending on which limma functions you are using. The",
1280 "Cite no. 3 only if sample weights were used.") 1346 "methodology articles are listed in Section 2.1 of the",
1281 cit[4] <- paste("Smyth GK (2005). Limma: linear models for microarray data.", 1347 link[1],
1282 "In: 'Bioinformatics and Computational Biology Solutions using", 1348 "Cite no. 3 only if sample weights were used."
1283 "R and Bioconductor'. R. Gentleman, V. Carey, S. doit,.", 1349 )
1284 "Irizarry, W. Huber (eds), Springer, New York, pages 397-420.") 1350 cit[4] <- paste(
1285 cit[5] <- paste("Please cite the first paper for the software itself and the", 1351 "Smyth GK (2005). Limma: linear models for microarray data.",
1286 "other papers for the various original statistical methods", 1352 "In: 'Bioinformatics and Computational Biology Solutions using",
1287 "implemented in edgeR. See Section 1.2 in the", link[2], 1353 "R and Bioconductor'. R. Gentleman, V. Carey, S. doit,.",
1288 "for more detail.") 1354 "Irizarry, W. Huber (eds), Springer, New York, pages 397-420."
1289 cit[6] <- paste("Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a", 1355 )
1290 "Bioconductor package for differential expression analysis", 1356 cit[5] <- paste(
1291 "of digital gene expression data. Bioinformatics 26, 139-140") 1357 "Please cite the first paper for the software itself and the",
1292 cit[7] <- paste("Robinson MD and Smyth GK (2007). Moderated statistical tests", 1358 "other papers for the various original statistical methods",
1293 "for assessing differences in tag abundance. Bioinformatics", 1359 "implemented in edgeR. See Section 1.2 in the", link[2],
1294 "23, 2881-2887") 1360 "for more detail."
1295 cit[8] <- paste("Robinson MD and Smyth GK (2008). Small-sample estimation of", 1361 )
1296 "negative binomial dispersion, with applications to SAGE data.", 1362 cit[6] <- paste(
1297 "Biostatistics, 9, 321-332") 1363 "Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a",
1298 cit[9] <- paste("McCarthy DJ, Chen Y and Smyth GK (2012). Differential", 1364 "Bioconductor package for differential expression analysis",
1299 "expression analysis of multifactor RNA-Seq experiments with", 1365 "of digital gene expression data. Bioinformatics 26, 139-140"
1300 "respect to biological variation. Nucleic Acids Research 40,", 1366 )
1301 "4288-4297") 1367 cit[7] <- paste(
1302 cit[10] <- paste("Law CW, Chen Y, Shi W, and Smyth GK (2014). Voom:", 1368 "Robinson MD and Smyth GK (2007). Moderated statistical tests",
1303 "precision weights unlock linear model analysis tools for", 1369 "for assessing differences in tag abundance. Bioinformatics",
1304 "RNA-seq read counts. Genome Biology 15, R29.") 1370 "23, 2881-2887"
1305 cit[11] <- paste("Ritchie ME, Diyagama D, Neilson J, van Laar R,", 1371 )
1306 "Dobrovic A, Holloway A and Smyth GK (2006).", 1372 cit[8] <- paste(
1307 "Empirical array quality weights for microarray data.", 1373 "Robinson MD and Smyth GK (2008). Small-sample estimation of",
1308 "BMC Bioinformatics 7, Article 261.") 1374 "negative binomial dispersion, with applications to SAGE data.",
1375 "Biostatistics, 9, 321-332"
1376 )
1377 cit[9] <- paste(
1378 "McCarthy DJ, Chen Y and Smyth GK (2012). Differential",
1379 "expression analysis of multifactor RNA-Seq experiments with",
1380 "respect to biological variation. Nucleic Acids Research 40,",
1381 "4288-4297"
1382 )
1383 cit[10] <- paste(
1384 "Law CW, Chen Y, Shi W, and Smyth GK (2014). Voom:",
1385 "precision weights unlock linear model analysis tools for",
1386 "RNA-seq read counts. Genome Biology 15, R29."
1387 )
1388 cit[11] <- paste(
1389 "Ritchie ME, Diyagama D, Neilson J, van Laar R,",
1390 "Dobrovic A, Holloway A and Smyth GK (2006).",
1391 "Empirical array quality weights for microarray data.",
1392 "BMC Bioinformatics 7, Article 261."
1393 )
1309 cata("<h3>Citations</h3>\n") 1394 cata("<h3>Citations</h3>\n")
1310 cata(cit[1], "\n") 1395 cata(cit[1], "\n")
1311 cata("<br>\n") 1396 cata("<br>\n")
1312 cata(cit[2], "\n") 1397 cata(cit[2], "\n")
1313 1398
1336 } 1421 }
1337 } 1422 }
1338 1423
1339 cata("<table border=\"0\">\n") 1424 cata("<table border=\"0\">\n")
1340 cata("<tr>\n") 1425 cata("<tr>\n")
1341 table_item("Task started at:"); table_item(time_start) 1426 table_item("Task started at:")
1427 table_item(time_start)
1342 cata("</tr>\n") 1428 cata("</tr>\n")
1343 cata("<tr>\n") 1429 cata("<tr>\n")
1344 table_item("Task ended at:"); table_item(time_end) 1430 table_item("Task ended at:")
1431 table_item(time_end)
1345 cata("</tr>\n") 1432 cata("</tr>\n")
1346 cata("<tr>\n") 1433 cata("<tr>\n")
1347 table_item("Task run time:"); table_item(time_taken) 1434 table_item("Task run time:")
1435 table_item(time_taken)
1348 cata("<tr>\n") 1436 cata("<tr>\n")
1349 cata("</table>\n") 1437 cata("</table>\n")
1350 1438
1351 cata("</body>\n") 1439 cata("</body>\n")
1352 cata("</html>") 1440 cata("</html>")