Mercurial > repos > iuc > limma_voom

comparison limma_voom.R @ 22:708348a17fa1 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 8c4b8f8711df0a4d6fa98fa8a6f91b977395d62c"
author iuc
date Sat, 26 Jun 2021 19:10:16 +0000
parents 58c35179ebf0
children 119b069fc845
comparison
equal deleted inserted replaced
21:58c35179ebf0 22:708348a17fa1
314 counts <- counts[, -1] 314 counts <- counts[, -1]
315 countsrows <- nrow(counts) 315 countsrows <- nrow(counts)
316 316
317 # Process factors 317 # Process factors
318 if (is.null(opt$factInput)) { 318 if (is.null(opt$factInput)) {
319 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, stringsAsFactors = TRUE)
320 if (!setequal(factordata[, 1], colnames(counts))) 320 if (!setequal(factordata[, 1], colnames(counts)))
321 stop("Sample IDs in counts and factors files don't match") 321 stop("Sample IDs in counts and factors files don't match")
322 # order samples as in counts matrix 322 # order samples as in counts matrix
323 factordata <- factordata[match(colnames(counts), factordata[, 1]), ] 323 factordata <- factordata[match(colnames(counts), factordata[, 1]), ]
324 factors <- factordata[, -1, drop = FALSE] 324 factors <- factordata[, -1, drop = FALSE]
343 if (nrow(factors) != ncol(counts)) { 343 if (nrow(factors) != ncol(counts)) {
344 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")
345 } 345 }
346 # make groups valid R names, required for makeContrasts 346 # make groups valid R names, required for makeContrasts
347 factors <- sapply(factors, make.names) 347 factors <- sapply(factors, make.names)
348 factors <- data.frame(factors) 348 factors <- data.frame(factors, stringsAsFactors = TRUE)
349 349
350 # if annotation file provided 350 # if annotation file provided
351 if (have_anno) { 351 if (have_anno) {
352 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)
353 } 353 }
473 if (filt_totcount) { 473 if (filt_totcount) {
474 keep <- rowSums(data$counts) >= opt$cntReq 474 keep <- rowSums(data$counts) >= opt$cntReq
475 } else if (filt_smpcount) { 475 } else if (filt_smpcount) {
476 keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq 476 keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq
477 } else if (filt_cpm) { 477 } else if (filt_cpm) {
478 478 cpms <- cpm(data$counts)
479 keep <- rowSums(cpm(data$counts) >= opt$cpmReq) >= opt$sampleReq 479 thresh <- cpms >= opt$cpmReq
480 keep <- rowSums(thresh) >= opt$sampleReq
480 481
481 if ("c" %in% plots) { 482 if ("c" %in% plots) {
482 # Plot CPM vs raw counts (to check threshold) 483 # Plot CPM vs raw counts (to check threshold)
483 pdf(cpm_pdf, width = 6.5, height = 10) 484 pdf(cpm_pdf, width = 6.5, height = 10)
484 par(mfrow = c(3, 2)) 485 par(mfrow = c(3, 2))
485 for (i in seq_len(nsamples)) { 486 for (i in seq_len(nsamples)) {
486 plot(data$counts[, i], myCPM[, i], xlim = c(0, 50), ylim = c(0, 3), main = samplenames[i], xlab = "Raw counts", ylab = "CPM") 487 plot(data$counts[, i], cpms[, i], xlim = c(0, 50), ylim = c(0, 3), main = samplenames[i], xlab = "Raw counts", ylab = "CPM")
487 abline(v = 10, col = "red", lty = 2, lwd = 2) 488 abline(v = 10, col = "red", lty = 2, lwd = 2)
488 abline(h = opt$cpmReq, col = 4) 489 abline(h = opt$cpmReq, col = 4)
489 } 490 }
490 link_name <- "CpmPlots.pdf" 491 link_name <- "CpmPlots.pdf"
491 link_addr <- "cpmplots.pdf" 492 link_addr <- "cpmplots.pdf"
670 } 671 }
671 } 672 }
672 673
673 a1 <- suppressWarnings(cmdscale(as.dist(dd), k = min(ndim, 8), eig = TRUE)) 674 a1 <- suppressWarnings(cmdscale(as.dist(dd), k = min(ndim, 8), eig = TRUE))
674 eigen <- data.frame(name = 1:min(ndim, 8), eigen = round(a1$eig[1:min(ndim, 8)] / sum(a1$eig), 2)) 675 eigen <- data.frame(name = 1:min(ndim, 8), eigen = round(a1$eig[1:min(ndim, 8)] / sum(a1$eig), 2))
675
676 png(mdsscree_png, width = 1000, height = 500) 676 png(mdsscree_png, width = 1000, height = 500)
677 par(mfrow = c(1, 2)) 677 par(mfrow = c(1, 2))
678 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")
679 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)
680 img_name <- paste0("MDSPlot_", names(factors)[1], ".png") 680 img_name <- paste0("MDSPlot_", names(factors)[1], ".png")