comparison limma_voom.R @ 22:708348a17fa1 draft

"planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/limma_voom commit 8c4b8f8711df0a4d6fa98fa8a6f91b977395d62c"

21:58c35179ebf0

    counts <- counts[, -1]
    countsrows <- nrow(counts)

    # Process factors
    if (is.null(opt$factInput)) {
        factordata <- read.table(opt$factFile, header = TRUE, sep = "\t", strip.white = TRUE)
        if (!setequal(factordata[, 1], colnames(counts)))
            stop("Sample IDs in counts and factors files don't match")
        # order samples as in counts matrix
        factordata <- factordata[match(colnames(counts), factordata[, 1]), ]
        factors <- factordata[, -1, drop = FALSE]

if (nrow(factors) != ncol(counts)) {
    stop("There are a different number of samples in the counts files and factors")
}
# make groups valid R names, required for makeContrasts
factors <- sapply(factors, make.names)
factors <- data.frame(factors)

 # if annotation file provided
if (have_anno) {
    geneanno <- read.table(opt$annoPath, header = TRUE, sep = "\t", quote = "", strip.white = TRUE, stringsAsFactors = FALSE)
}

    if (filt_totcount) {
        keep <- rowSums(data$counts) >= opt$cntReq
    } else if (filt_smpcount) {
        keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq
    } else if (filt_cpm) {

        keep <- rowSums(cpm(data$counts) >= opt$cpmReq) >= opt$sampleReq

        if ("c" %in% plots) {
            # Plot CPM vs raw counts (to check threshold)
            pdf(cpm_pdf, width = 6.5, height = 10)
            par(mfrow = c(3, 2))
            for (i in seq_len(nsamples)) {
                plot(data$counts[, i], myCPM[, i], xlim = c(0, 50), ylim = c(0, 3), main = samplenames[i], xlab = "Raw counts", ylab = "CPM")
                abline(v = 10, col = "red", lty = 2, lwd = 2)
                abline(h = opt$cpmReq, col = 4)
            }
            link_name <- "CpmPlots.pdf"
            link_addr <- "cpmplots.pdf"

  }
}

a1 <- suppressWarnings(cmdscale(as.dist(dd), k = min(ndim, 8), eig = TRUE))
eigen <- data.frame(name = 1:min(ndim, 8), eigen = round(a1$eig[1:min(ndim, 8)] / sum(a1$eig), 2))

png(mdsscree_png, width = 1000, height = 500)
par(mfrow = c(1, 2))
plotMDS(y, labels = samplenames, col = as.numeric(factors[, 1]), main = "MDS Plot: Dims 1 and 2")
barplot(eigen$eigen, names.arg = eigen$name,  main = "Scree Plot: Variance Explained", xlab = "Dimension", ylab = "Proportion", las = 1)
img_name <- paste0("MDSPlot_", names(factors)[1], ".png")

22:708348a17fa1

    counts <- counts[, -1]
    countsrows <- nrow(counts)

    # Process factors
    if (is.null(opt$factInput)) {
        factordata <- read.table(opt$factFile, header = TRUE, sep = "\t", strip.white = TRUE, stringsAsFactors = TRUE)
        if (!setequal(factordata[, 1], colnames(counts)))
            stop("Sample IDs in counts and factors files don't match")
        # order samples as in counts matrix
        factordata <- factordata[match(colnames(counts), factordata[, 1]), ]
        factors <- factordata[, -1, drop = FALSE]

if (nrow(factors) != ncol(counts)) {
    stop("There are a different number of samples in the counts files and factors")
}
# make groups valid R names, required for makeContrasts
factors <- sapply(factors, make.names)
factors <- data.frame(factors, stringsAsFactors = TRUE)

 # if annotation file provided
if (have_anno) {
    geneanno <- read.table(opt$annoPath, header = TRUE, sep = "\t", quote = "", strip.white = TRUE, stringsAsFactors = FALSE)
}

    if (filt_totcount) {
        keep <- rowSums(data$counts) >= opt$cntReq
    } else if (filt_smpcount) {
        keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq
    } else if (filt_cpm) {
        cpms <- cpm(data$counts)
        thresh <- cpms >= opt$cpmReq
        keep <- rowSums(thresh) >= opt$sampleReq

        if ("c" %in% plots) {
            # Plot CPM vs raw counts (to check threshold)
            pdf(cpm_pdf, width = 6.5, height = 10)
            par(mfrow = c(3, 2))
            for (i in seq_len(nsamples)) {
                plot(data$counts[, i], cpms[, i], xlim = c(0, 50), ylim = c(0, 3), main = samplenames[i], xlab = "Raw counts", ylab = "CPM")
                abline(v = 10, col = "red", lty = 2, lwd = 2)
                abline(h = opt$cpmReq, col = 4)
            }
            link_name <- "CpmPlots.pdf"
            link_addr <- "cpmplots.pdf"

  }
}

a1 <- suppressWarnings(cmdscale(as.dist(dd), k = min(ndim, 8), eig = TRUE))
eigen <- data.frame(name = 1:min(ndim, 8), eigen = round(a1$eig[1:min(ndim, 8)] / sum(a1$eig), 2))
png(mdsscree_png, width = 1000, height = 500)
par(mfrow = c(1, 2))
plotMDS(y, labels = samplenames, col = as.numeric(factors[, 1]), main = "MDS Plot: Dims 1 and 2")
barplot(eigen$eigen, names.arg = eigen$name,  main = "Scree Plot: Variance Explained", xlab = "Dimension", ylab = "Proportion", las = 1)
img_name <- paste0("MDSPlot_", names(factors)[1], ".png")