Mercurial > repos > iuc > edger

comparison edger.R @ 3:2aa9dd87aad3 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 e646be741e315df9332b5206cec1e47c11370ff1
author iuc
date Sun, 06 May 2018 13:38:22 -0400
parents 3de72f1f78aa
children 9a62dbdb122b
comparison
equal deleted inserted replaced
2:3de72f1f78aa 3:2aa9dd87aad3
304 304
305 bcvOutPdf <- makeOut("bcvplot.pdf") 305 bcvOutPdf <- makeOut("bcvplot.pdf")
306 bcvOutPng <- makeOut("bcvplot.png") 306 bcvOutPng <- makeOut("bcvplot.png")
307 qlOutPdf <- makeOut("qlplot.pdf") 307 qlOutPdf <- makeOut("qlplot.pdf")
308 qlOutPng <- makeOut("qlplot.png") 308 qlOutPng <- makeOut("qlplot.png")
309 mdsOutPdf <- makeOut("mdsplot.pdf") 309 mdsOutPdf <- character() # Initialise character vector
310 mdsOutPng <- makeOut("mdsplot.png") 310 mdsOutPng <- character()
311 mdOutPdf <- character() # Initialise character vector 311 for (i in 1:ncol(factors)) {
312 mdsOutPdf[i] <- makeOut(paste0("mdsplot_", names(factors)[i], ".pdf"))
313 mdsOutPng[i] <- makeOut(paste0("mdsplot_", names(factors)[i], ".png"))
314 }
315 mdOutPdf <- character()
312 mdOutPng <- character() 316 mdOutPng <- character()
313 topOut <- character() 317 topOut <- character()
314 for (i in 1:length(contrastData)) { 318 for (i in 1:length(contrastData)) {
315 mdOutPdf[i] <- makeOut(paste0("mdplot_", contrastData[i], ".pdf")) 319 mdOutPdf[i] <- makeOut(paste0("mdplot_", contrastData[i], ".pdf"))
316 mdOutPng[i] <- makeOut(paste0("mdplot_", contrastData[i], ".png")) 320 mdOutPng[i] <- makeOut(paste0("mdplot_", contrastData[i], ".png"))
336 340
337 # Extract counts and annotation data 341 # Extract counts and annotation data
338 data <- list() 342 data <- list()
339 data$counts <- counts 343 data$counts <- counts
340 if (haveAnno) { 344 if (haveAnno) {
341 data$genes <- geneanno 345 # order annotation by genes in counts (assumes gene ids are in 1st column of geneanno)
346 annoord <- geneanno[match(row.names(counts), geneanno[,1]), ]
347 data$genes <- annoord
342 } else { 348 } else {
343 data$genes <- data.frame(GeneID=row.names(counts)) 349 data$genes <- data.frame(GeneID=row.names(counts))
344 } 350 }
345 351
346 # If filter crieteria set, filter out genes that do not have a required cpm/counts in a required number of 352 # If filter crieteria set, filter out genes that do not have a required cpm/counts in a required number of
408 ### Data Output 414 ### Data Output
409 ################################################################################ 415 ################################################################################
410 416
411 # Plot MDS 417 # Plot MDS
412 labels <- names(counts) 418 labels <- names(counts)
419
420 # MDS plot
413 png(mdsOutPng, width=600, height=600) 421 png(mdsOutPng, width=600, height=600)
414 # Currently only using a single factor 422 plotMDS(data, labels=labels, col=as.numeric(factors[, 1]), cex=0.8, main=paste("MDS Plot:", names(factors)[1]))
415 plotMDS(data, labels=labels, col=as.numeric(factors[, 1]), cex=0.8, main="MDS Plot") 423 imgName <- paste0("MDS Plot_", names(factors)[1], ".png")
416 imageData[1, ] <- c("MDS Plot", "mdsplot.png") 424 imgAddr <- paste0("mdsplot_", names(factors)[1], ".png")
425 imageData[1, ] <- c(imgName, imgAddr)
417 invisible(dev.off()) 426 invisible(dev.off())
418 427
419 pdf(mdsOutPdf) 428 pdf(mdsOutPdf)
420 plotMDS(data, labels=labels, cex=0.5) 429 plotMDS(data, labels=labels, col=as.numeric(factors[, 1]), cex=0.8, main=paste("MDS Plot:", names(factors)[1]))
421 linkData[1, ] <- c("MDS Plot.pdf", "mdsplot.pdf") 430 linkName <- paste0("MDS Plot_", names(factors)[1], ".pdf")
431 linkAddr <- paste0("mdsplot_", names(factors)[1], ".pdf")
432 linkData[1, ] <- c(linkName, linkAddr)
422 invisible(dev.off()) 433 invisible(dev.off())
434
435 # If additional factors create additional MDS plots coloured by factor
436 if (ncol(factors) > 1) {
437 for (i in 2:ncol(factors)) {
438 png(mdsOutPng[i], width=600, height=600)
439 plotMDS(data, labels=labels, col=as.numeric(factors[, i]), cex=0.8, main=paste("MDS Plot:", names(factors)[i]))
440 imgName <- paste0("MDS Plot_", names(factors)[i], ".png")
441 imgAddr <- paste0("mdsplot_", names(factors)[i], ".png")
442 imageData <- rbind(imageData, c(imgName, imgAddr))
443 invisible(dev.off())
444
445 pdf(mdsOutPdf[i])
446 plotMDS(data, labels=labels, col=as.numeric(factors[, i]), cex=0.8, main=paste("MDS Plot:", names(factors)[i]))
447 linkName <- paste0("MDS Plot_", names(factors)[i], ".pdf")
448 linkAddr <- paste0("mdsplot_", names(factors)[i], ".pdf")
449 linkData <- rbind(linkData, c(linkName, linkAddr))
450 invisible(dev.off())
451 }
452 }
423 453
424 # BCV Plot 454 # BCV Plot
425 png(bcvOutPng, width=600, height=600) 455 png(bcvOutPng, width=600, height=600)
426 plotBCV(data, main="BCV Plot") 456 plotBCV(data, main="BCV Plot")
427 imgName <- "BCV Plot" 457 imgName <- "BCV Plot"
467 497
468 # Save normalised counts (log2cpm) 498 # Save normalised counts (log2cpm)
469 if (wantNorm) { 499 if (wantNorm) {
470 normalisedCounts <- cpm(data, normalized.lib.sizes=TRUE, log=TRUE) 500 normalisedCounts <- cpm(data, normalized.lib.sizes=TRUE, log=TRUE)
471 normalisedCounts <- data.frame(data$genes, normalisedCounts) 501 normalisedCounts <- data.frame(data$genes, normalisedCounts)
472 write.table (normalisedCounts, file=normOut, row.names=FALSE, sep="\t") 502 write.table (normalisedCounts, file=normOut, row.names=FALSE, sep="\t", quote=FALSE)
473 linkData <- rbind(linkData, c("edgeR_normcounts.tsv", "edgeR_normcounts.tsv")) 503 linkData <- rbind(linkData, c("edgeR_normcounts.tsv", "edgeR_normcounts.tsv"))
474 } 504 }
475 505
476 506
477 for (i in 1:length(contrastData)) { 507 for (i in 1:length(contrastData)) {
490 downCount[i] <- sumStatus["Down", ] 520 downCount[i] <- sumStatus["Down", ]
491 flatCount[i] <- sumStatus["NotSig", ] 521 flatCount[i] <- sumStatus["NotSig", ]
492 522
493 # Write top expressions table 523 # Write top expressions table
494 top <- topTags(res, n=Inf, sort.by="PValue") 524 top <- topTags(res, n=Inf, sort.by="PValue")
495 write.table(top, file=topOut[i], row.names=FALSE, sep="\t") 525 write.table(top, file=topOut[i], row.names=FALSE, sep="\t", quote=FALSE)
496 526
497 linkName <- paste0("edgeR_", contrastData[i], ".tsv") 527 linkName <- paste0("edgeR_", contrastData[i], ".tsv")
498 linkAddr <- paste0("edgeR_", contrastData[i], ".tsv") 528 linkAddr <- paste0("edgeR_", contrastData[i], ".tsv")
499 linkData <- rbind(linkData, c(linkName, linkAddr)) 529 linkData <- rbind(linkData, c(linkName, linkAddr))
500 530
501 # Plot MD (log ratios vs mean difference) using limma package 531 # Plot MD (log ratios vs mean difference) using limma package
502 pdf(mdOutPdf[i]) 532 pdf(mdOutPdf[i])
503 limma::plotMD(res, status=status, 533 limma::plotMD(res, status=status,
504 main=paste("MD Plot:", unmake.names(contrastData[i])), 534 main=paste("MD Plot:", unmake.names(contrastData[i])),
505 col=alpha(c("firebrick", "blue"), 0.4), values=c("1", "-1"), 535 hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1),
506 xlab="Average Expression", ylab="logFC") 536 xlab="Average Expression", ylab="logFC")
507 537
508 abline(h=0, col="grey", lty=2) 538 abline(h=0, col="grey", lty=2)
509 539
510 linkName <- paste0("MD Plot_", contrastData[i], ".pdf") 540 linkName <- paste0("MD Plot_", contrastData[i], ".pdf")
513 invisible(dev.off()) 543 invisible(dev.off())
514 544
515 png(mdOutPng[i], height=600, width=600) 545 png(mdOutPng[i], height=600, width=600)
516 limma::plotMD(res, status=status, 546 limma::plotMD(res, status=status,
517 main=paste("MD Plot:", unmake.names(contrastData[i])), 547 main=paste("MD Plot:", unmake.names(contrastData[i])),
518 col=alpha(c("firebrick", "blue"), 0.4), values=c("1", "-1"), 548 hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1),
519 xlab="Average Expression", ylab="logFC") 549 xlab="Average Expression", ylab="logFC")
520 550
521 abline(h=0, col="grey", lty=2) 551 abline(h=0, col="grey", lty=2)
522 552
523 imgName <- paste0("MD Plot_", contrastData[i], ".png") 553 imgName <- paste0("MD Plot_", contrastData[i], ".png")
618 cata("<h4>Additional Information</h4>\n") 648 cata("<h4>Additional Information</h4>\n")
619 cata("<ul>\n") 649 cata("<ul>\n")
620 650
621 if (filtCPM || filtSmpCount || filtTotCount) { 651 if (filtCPM || filtSmpCount || filtTotCount) {
622 if (filtCPM) { 652 if (filtCPM) {
623 tempStr <- paste("Genes without more than", opt$cmpReq, 653 tempStr <- paste("Genes without more than", opt$cpmReq,
624 "CPM in at least", opt$sampleReq, "samples are insignificant", 654 "CPM in at least", opt$sampleReq, "samples are insignificant",
625 "and filtered out.") 655 "and filtered out.")
626 } else if (filtSmpCount) { 656 } else if (filtSmpCount) {
627 tempStr <- paste("Genes without more than", opt$cntReq, 657 tempStr <- paste("Genes without more than", opt$cntReq,
628 "counts in at least", opt$sampleReq, "samples are insignificant", 658 "counts in at least", opt$sampleReq, "samples are insignificant",