Mercurial > repos > iuc > edger

comparison edger.R @ 15:5bf899c13979 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/edger commit 025a0d70be99e50ae89ad78a70cc8c51571ea1a2
author iuc
date Wed, 22 Nov 2023 03:57:37 +0000
parents 070900306913
children ae2aad0a6d50
comparison
equal deleted inserted replaced
14:070900306913 15:5bf899c13979
83 unmake_names <- function(string) { 83 unmake_names <- function(string) {
84 string <- gsub(".", " ", string, fixed = TRUE) 84 string <- gsub(".", " ", string, fixed = TRUE)
85 return(string) 85 return(string)
86 } 86 }
87 87
88 # Sanitise file base names coming from factors or contrasts
89 sanitise_basename <- function(string) {
90 string <- gsub("[/^]", "_", string)
91 return(string)
92 }
93
88 # Generate output folder and paths 94 # Generate output folder and paths
89 make_out <- function(filename) { 95 make_out <- function(filename) {
90 return(paste0(out_path, "/", filename)) 96 return(paste0(out_path, "/", filename))
91 } 97 }
92 98
329 ql_pdf <- make_out("qlplot.pdf") 335 ql_pdf <- make_out("qlplot.pdf")
330 ql_png <- make_out("qlplot.png") 336 ql_png <- make_out("qlplot.png")
331 mds_pdf <- character() # Initialise character vector 337 mds_pdf <- character() # Initialise character vector
332 mds_png <- character() 338 mds_png <- character()
333 for (i in seq_len(ncol(factors))) { 339 for (i in seq_len(ncol(factors))) {
334 mds_pdf[i] <- make_out(paste0("mdsplot_", names(factors)[i], ".pdf")) 340 mds_pdf[i] <- make_out(paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".pdf"))
335 mds_png[i] <- make_out(paste0("mdsplot_", names(factors)[i], ".png")) 341 mds_png[i] <- make_out(paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".png"))
336 } 342 }
337 md_pdf <- character() 343 md_pdf <- character()
338 md_png <- character() 344 md_png <- character()
339 top_out <- character() 345 top_out <- character()
340 for (i in seq_along(contrast_data)) { 346 for (i in seq_along(contrast_data)) {
341 md_pdf[i] <- make_out(paste0("mdplot_", contrast_data[i], ".pdf")) 347 md_pdf[i] <- make_out(paste0("mdplot_", sanitise_basename(contrast_data[i]), ".pdf"))
342 md_png[i] <- make_out(paste0("mdplot_", contrast_data[i], ".png")) 348 md_png[i] <- make_out(paste0("mdplot_", sanitise_basename(contrast_data[i]), ".png"))
343 top_out[i] <- make_out(paste0("edgeR_", contrast_data[i], ".tsv")) 349 top_out[i] <- make_out(paste0("edgeR_", sanitise_basename(contrast_data[i]), ".tsv"))
344 } # Save output paths for each contrast as vectors 350 } # Save output paths for each contrast as vectors
345 norm_out <- make_out("edgeR_normcounts.tsv") 351 norm_out <- make_out("edgeR_normcounts.tsv")
346 rda_out <- make_out("edgeR_analysis.RData") 352 rda_out <- make_out("edgeR_analysis.RData")
347 session_out <- make_out("session_info.txt") 353 session_out <- make_out("session_info.txt")
348 354
444 labels <- names(counts) 450 labels <- names(counts)
445 451
446 # MDS plot 452 # MDS plot
447 png(mds_png, width = 600, height = 600) 453 png(mds_png, width = 600, height = 600)
448 plotMDS(data, labels = labels, col = as.numeric(factors[, 1]), cex = 0.8, main = paste("MDS Plot:", names(factors)[1])) 454 plotMDS(data, labels = labels, col = as.numeric(factors[, 1]), cex = 0.8, main = paste("MDS Plot:", names(factors)[1]))
449 img_name <- paste0("MDS Plot_", names(factors)[1], ".png") 455 img_name <- paste0("MDS Plot_", sanitise_basename(names(factors)[1]), ".png")
450 img_addr <- paste0("mdsplot_", names(factors)[1], ".png") 456 img_addr <- paste0("mdsplot_", sanitise_basename(names(factors)[1]), ".png")
451 image_data[1, ] <- c(img_name, img_addr) 457 image_data[1, ] <- c(img_name, img_addr)
452 invisible(dev.off()) 458 invisible(dev.off())
453 459
454 pdf(mds_pdf) 460 pdf(mds_pdf)
455 plotMDS(data, labels = labels, col = as.numeric(factors[, 1]), cex = 0.8, main = paste("MDS Plot:", names(factors)[1])) 461 plotMDS(data, labels = labels, col = as.numeric(factors[, 1]), cex = 0.8, main = paste("MDS Plot:", names(factors)[1]))
456 link_name <- paste0("MDS Plot_", names(factors)[1], ".pdf") 462 link_name <- paste0("MDS Plot_", sanitise_basename(names(factors)[1]), ".pdf")
457 link_addr <- paste0("mdsplot_", names(factors)[1], ".pdf") 463 link_addr <- paste0("mdsplot_", sanitise_basename(names(factors)[1]), ".pdf")
458 link_data[1, ] <- c(link_name, link_addr) 464 link_data[1, ] <- c(link_name, link_addr)
459 invisible(dev.off()) 465 invisible(dev.off())
460 466
461 # If additional factors create additional MDS plots coloured by factor 467 # If additional factors create additional MDS plots coloured by factor
462 if (ncol(factors) > 1) { 468 if (ncol(factors) > 1) {
463 for (i in 2:ncol(factors)) { 469 for (i in 2:ncol(factors)) {
464 png(mds_png[i], width = 600, height = 600) 470 png(mds_png[i], width = 600, height = 600)
465 plotMDS(data, labels = labels, col = as.numeric(factors[, i]), cex = 0.8, main = paste("MDS Plot:", names(factors)[i])) 471 plotMDS(data, labels = labels, col = as.numeric(factors[, i]), cex = 0.8, main = paste("MDS Plot:", names(factors)[i]))
466 img_name <- paste0("MDS Plot_", names(factors)[i], ".png") 472 img_name <- paste0("MDS Plot_", sanitise_basename(names(factors)[i]), ".png")
467 img_addr <- paste0("mdsplot_", names(factors)[i], ".png") 473 img_addr <- paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".png")
468 image_data <- rbind(image_data, c(img_name, img_addr)) 474 image_data <- rbind(image_data, c(img_name, img_addr))
469 invisible(dev.off()) 475 invisible(dev.off())
470 476
471 pdf(mds_pdf[i]) 477 pdf(mds_pdf[i])
472 plotMDS(data, labels = labels, col = as.numeric(factors[, i]), cex = 0.8, main = paste("MDS Plot:", names(factors)[i])) 478 plotMDS(data, labels = labels, col = as.numeric(factors[, i]), cex = 0.8, main = paste("MDS Plot:", names(factors)[i]))
473 link_name <- paste0("MDS Plot_", names(factors)[i], ".pdf") 479 link_name <- paste0("MDS Plot_", sanitise_basename(names(factors)[i]), ".pdf")
474 link_addr <- paste0("mdsplot_", names(factors)[i], ".pdf") 480 link_addr <- paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".pdf")
475 link_data <- rbind(link_data, c(link_name, link_addr)) 481 link_data <- rbind(link_data, c(link_name, link_addr))
476 invisible(dev.off()) 482 invisible(dev.off())
477 } 483 }
478 } 484 }
479 485
547 553
548 # Write top expressions table 554 # Write top expressions table
549 top <- topTags(res, adjust.method = opt$pAdjOpt, n = Inf, sort.by = "PValue") 555 top <- topTags(res, adjust.method = opt$pAdjOpt, n = Inf, sort.by = "PValue")
550 write.table(top, file = top_out[i], row.names = FALSE, sep = "\t", quote = FALSE) 556 write.table(top, file = top_out[i], row.names = FALSE, sep = "\t", quote = FALSE)
551 557
552 link_name <- paste0("edgeR_", contrast_data[i], ".tsv") 558 link_name <- paste0("edgeR_", sanitise_basename(contrast_data[i]), ".tsv")
553 link_addr <- paste0("edgeR_", contrast_data[i], ".tsv") 559 link_addr <- paste0("edgeR_", sanitise_basename(contrast_data[i]), ".tsv")
554 link_data <- rbind(link_data, c(link_name, link_addr)) 560 link_data <- rbind(link_data, c(link_name, link_addr))
555 561
556 # Plot MD (log ratios vs mean difference) using limma package 562 # Plot MD (log ratios vs mean difference) using limma package
557 pdf(md_pdf[i]) 563 pdf(md_pdf[i])
558 limma::plotMD(res, 564 limma::plotMD(res,
562 xlab = "Average Expression", ylab = "logFC" 568 xlab = "Average Expression", ylab = "logFC"
563 ) 569 )
564 570
565 abline(h = 0, col = "grey", lty = 2) 571 abline(h = 0, col = "grey", lty = 2)
566 572
567 link_name <- paste0("MD Plot_", contrast_data[i], ".pdf") 573 link_name <- paste0("MD Plot_", sanitise_basename(contrast_data[i]), ".pdf")
568 link_addr <- paste0("mdplot_", contrast_data[i], ".pdf") 574 link_addr <- paste0("mdplot_", sanitise_basename(contrast_data[i]), ".pdf")
569 link_data <- rbind(link_data, c(link_name, link_addr)) 575 link_data <- rbind(link_data, c(link_name, link_addr))
570 invisible(dev.off()) 576 invisible(dev.off())
571 577
572 png(md_png[i], height = 600, width = 600) 578 png(md_png[i], height = 600, width = 600)
573 limma::plotMD(res, 579 limma::plotMD(res,
577 xlab = "Average Expression", ylab = "logFC" 583 xlab = "Average Expression", ylab = "logFC"
578 ) 584 )
579 585
580 abline(h = 0, col = "grey", lty = 2) 586 abline(h = 0, col = "grey", lty = 2)
581 587
582 img_name <- paste0("MD Plot_", contrast_data[i], ".png") 588 img_name <- paste0("MD Plot_", sanitise_basename(contrast_data[i]), ".png")
583 img_addr <- paste0("mdplot_", contrast_data[i], ".png") 589 img_addr <- paste0("mdplot_", sanitise_basename(contrast_data[i]), ".png")
584 image_data <- rbind(image_data, c(img_name, img_addr)) 590 image_data <- rbind(image_data, c(img_name, img_addr))
585 invisible(dev.off()) 591 invisible(dev.off())
586 } 592 }
587 sig_diff <- data.frame(Up = up_count, Flat = flat_count, Down = down_count) 593 sig_diff <- data.frame(Up = up_count, Flat = flat_count, Down = down_count)
588 row.names(sig_diff) <- contrast_data 594 row.names(sig_diff) <- contrast_data