Mercurial > repos > iuc > limma_voom

comparison limma_voom.R @ 11:7e8af58c8052 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 4bcbd83574ecf7194f3370aa883a9573563afdbc
author iuc
date Mon, 11 Jun 2018 08:18:25 -0400
parents e26047c4562d
children 81796eb60bd0
comparison
equal deleted inserted replaced
10:e26047c4562d 11:7e8af58c8052
175 "robOpt", "b", 0, "logical", 175 "robOpt", "b", 0, "logical",
176 "trend", "t", 1, "double", 176 "trend", "t", 1, "double",
177 "weightOpt", "w", 0, "logical", 177 "weightOpt", "w", 0, "logical",
178 "topgenes", "G", 1, "integer", 178 "topgenes", "G", 1, "integer",
179 "treatOpt", "T", 0, "logical", 179 "treatOpt", "T", 0, "logical",
180 "plots", "P", 1, "character"), 180 "plots", "P", 1, "character",
181 "libinfoOpt", "L", 0, "logical"),
181 byrow=TRUE, ncol=4) 182 byrow=TRUE, ncol=4)
182 opt <- getopt(spec) 183 opt <- getopt(spec)
183 184
184 185
185 if (is.null(opt$matrixPath) & is.null(opt$filesPath)) { 186 if (is.null(opt$matrixPath) & is.null(opt$filesPath)) {
246 247
247 if (is.null(opt$treatOpt)) { 248 if (is.null(opt$treatOpt)) {
248 wantTreat <- FALSE 249 wantTreat <- FALSE
249 } else { 250 } else {
250 wantTreat <- TRUE 251 wantTreat <- TRUE
252 }
253
254 if (is.null(opt$libinfoOpt)) {
255 wantLibinfo <- FALSE
256 } else {
257 wantLibinfo <- TRUE
251 } 258 }
252 259
253 260
254 if (!is.null(opt$filesPath)) { 261 if (!is.null(opt$filesPath)) {
255 # Process the separate count files (adapted from DESeq2 wrapper) 262 # Process the separate count files (adapted from DESeq2 wrapper)
281 countfiles <- lapply(sampleTable$filename, function(x){read.delim(x, row.names=1)}) 288 countfiles <- lapply(sampleTable$filename, function(x){read.delim(x, row.names=1)})
282 counts <- do.call("cbind", countfiles) 289 counts <- do.call("cbind", countfiles)
283 290
284 } else { 291 } else {
285 # Process the single count matrix 292 # Process the single count matrix
286 counts <- read.table(opt$matrixPath, header=TRUE, sep="\t", stringsAsFactors=FALSE) 293 counts <- read.table(opt$matrixPath, header=TRUE, sep="\t", strip.white=TRUE, stringsAsFactors=FALSE)
287 row.names(counts) <- counts[, 1] 294 row.names(counts) <- counts[, 1]
288 counts <- counts[ , -1] 295 counts <- counts[ , -1]
289 countsRows <- nrow(counts) 296 countsRows <- nrow(counts)
290 297
291 # Process factors 298 # Process factors
292 if (is.null(opt$factInput)) { 299 if (is.null(opt$factInput)) {
293 factorData <- read.table(opt$factFile, header=TRUE, sep="\t") 300 factorData <- read.table(opt$factFile, header=TRUE, sep="\t", strip.white=TRUE)
294 factors <- factorData[, -1, drop=FALSE] 301 factors <- factorData[, -1, drop=FALSE]
295 } else { 302 } else {
296 factors <- unlist(strsplit(opt$factInput, "|", fixed=TRUE)) 303 factors <- unlist(strsplit(opt$factInput, "|", fixed=TRUE))
297 factorData <- list() 304 factorData <- list()
298 for (fact in factors) { 305 for (fact in factors) {
311 } 318 }
312 } 319 }
313 320
314 # if annotation file provided 321 # if annotation file provided
315 if (haveAnno) { 322 if (haveAnno) {
316 geneanno <- read.table(opt$annoPath, header=TRUE, sep="\t", stringsAsFactors=FALSE) 323 geneanno <- read.table(opt$annoPath, header=TRUE, sep="\t", quote= "", strip.white=TRUE, stringsAsFactors=FALSE)
317 } 324 }
318 325
319 #Create output directory 326 #Create output directory
320 dir.create(opt$outPath, showWarnings=FALSE) 327 dir.create(opt$outPath, showWarnings=FALSE)
321 328
510 colnames(design) <- gsub(factorList[i], "", colnames(design), fixed=TRUE) 517 colnames(design) <- gsub(factorList[i], "", colnames(design), fixed=TRUE)
511 } 518 }
512 519
513 # Calculating normalising factors 520 # Calculating normalising factors
514 print("Calculating Normalisation Factors") 521 print("Calculating Normalisation Factors")
522 logcounts <- y #store for plots
515 y <- calcNormFactors(y, method=opt$normOpt) 523 y <- calcNormFactors(y, method=opt$normOpt)
516 524
517 # Generate contrasts information 525 # Generate contrasts information
518 print("Generating Contrasts") 526 print("Generating Contrasts")
519 contrasts <- makeContrasts(contrasts=contrastData, levels=design) 527 contrasts <- makeContrasts(contrasts=contrastData, levels=design)
651 # Plot MD plots for individual samples 659 # Plot MD plots for individual samples
652 print("Generating MD plots for samples") 660 print("Generating MD plots for samples")
653 pdf(mdsamOutPdf, width=6.5, height=10) 661 pdf(mdsamOutPdf, width=6.5, height=10)
654 par(mfrow=c(3, 2)) 662 par(mfrow=c(3, 2))
655 for (i in 1:nsamples) { 663 for (i in 1:nsamples) {
656 plotMD(y, column = i) 664 if (opt$normOpt != "none") {
665 plotMD(logcounts, column=i, main=paste(colnames(logcounts)[i], "(before)"))
666 abline(h=0, col="red", lty=2, lwd=2)
667 }
668 plotMD(y, column=i)
657 abline(h=0, col="red", lty=2, lwd=2) 669 abline(h=0, col="red", lty=2, lwd=2)
658 } 670 }
659 linkName <- "MDPlots_Samples.pdf" 671 linkName <- "MDPlots_Samples.pdf"
660 linkAddr <- "mdplots_samples.pdf" 672 linkAddr <- "mdplots_samples.pdf"
661 linkData <- rbind(linkData, c(linkName, linkAddr)) 673 linkData <- rbind(linkData, c(linkName, linkAddr))
759 fit <- eBayes(voomFit, robust=FALSE) 771 fit <- eBayes(voomFit, robust=FALSE)
760 } 772 }
761 plotData <- vData 773 plotData <- vData
762 } 774 }
763 775
776 # Save library size info
777 if (wantLibinfo) {
778 efflibsize <- round(y$samples$lib.size * y$samples$norm.factors)
779 libsizeinfo <- cbind(y$samples, EffectiveLibrarySize=efflibsize)
780 libsizeinfo$lib.size <- round(libsizeinfo$lib.size)
781 names(libsizeinfo)[names(libsizeinfo)=="sampleID"] <- "SampleID"
782 names(libsizeinfo)[names(libsizeinfo)=="lib.size"] <- "LibrarySize"
783 names(libsizeinfo)[names(libsizeinfo)=="norm.factors"] <- "NormalisationFactor"
784 write.table(libsizeinfo, file="libsizeinfo", row.names=FALSE, sep="\t", quote=FALSE)
785 }
764 786
765 print("Generating DE results") 787 print("Generating DE results")
766 788
767 if (wantTreat) { 789 if (wantTreat) {
768 print("Applying TREAT method") 790 print("Applying TREAT method")