| Previous changeset 1:c24765926774 (2019-03-26) Next changeset 3:d1f7fa5bb3cb (2021-07-23) |
|
Commit message:
"planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415" |
|
modified:
get_deseq_dataset.R ruvseq.R ruvseq.xml |
| b |
| diff -r c24765926774 -r fed9d0350d72 get_deseq_dataset.R --- a/get_deseq_dataset.R Tue Mar 26 06:25:38 2019 -0400 +++ b/get_deseq_dataset.R Fri Jan 15 17:53:15 2021 +0000 |
| [ |
| @@ -1,76 +1,80 @@ -get_deseq_dataset <- function(sampleTable, header, designFormula, tximport, txtype, tx2gene) { +get_deseq_dataset <- function(sample_table, header, design_formula, tximport, txtype, tx2gene) { dir <- "" - if (!is.null(header)) { - hasHeader <- TRUE - } else { - hasHeader <- FALSE - } - - if (!is.null(tximport)) { + has_header <- !is.null(header) + use_txi <- !is.null(tximport) + if (use_txi) { if (is.null(tx2gene)) stop("A transcript-to-gene map or a GTF/GFF3 file is required for tximport") if (tolower(file_ext(tx2gene)) == "gff") { - gffFile <-tx2gene + gff_file <- tx2gene } else { - gffFile <- NULL - tx2gene <- read.table(tx2gene, header=hasHeader) + gff_file <- NULL + tx2gene <- read.table(tx2gene, header = has_header) } - useTXI <- TRUE - } else { - useTXI <- FALSE } - if (!useTXI & hasHeader) { - countfiles <- lapply(as.character(sampleTable$filename), function(x){read.delim(x, row.names=1)}) + if (!use_txi & has_header) { + countfiles <- lapply(as.character(sample_table$filename), read.delim, row.names = 1) tbl <- do.call("cbind", countfiles) - colnames(tbl) <- rownames(sampleTable) # take sample ids from header + colnames(tbl) <- rownames(sample_table) # take sample ids from header # check for htseq report lines (from DESeqDataSetFromHTSeqCount function) - oldSpecialNames <- c("no_feature", "ambiguous", "too_low_aQual", - "not_aligned", "alignment_not_unique") - specialRows <- (substr(rownames(tbl), 1, 1) == "_") | rownames(tbl) %in% oldSpecialNames - tbl <- tbl[!specialRows, , drop = FALSE] + old_special_names <- c( + "no_feature", + "ambiguous", + "too_low_aQual", + "not_aligned", + "alignment_not_unique" + ) + special_rows <- (substr(rownames(tbl), 1, 1) == "_") | rownames(tbl) %in% old_special_names + tbl <- tbl[!special_rows, , drop = FALSE] - dds <- DESeqDataSetFromMatrix(countData = tbl, - colData = subset(sampleTable, select=-(filename)), - design = designFormula) - } else if (!useTXI & !hasHeader) { + dds <- DESeqDataSetFromMatrix( + countData = tbl, + colData = subset(sample_table, select = -filename), + design = design_formula + ) + } else if (!use_txi & !has_header) { # construct the object from HTSeq files - dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, - directory = dir, - design = designFormula) - colnames(dds) <- row.names(sampleTable) + dds <- DESeqDataSetFromHTSeqCount( + sampleTable = sample_table, + directory = dir, + design = design_formula + ) + colnames(dds) <- row.names(sample_table) } else { # construct the object using tximport library("tximport") - txiFiles <- as.character(sampleTable$filename) - labs <- row.names(sampleTable) - names(txiFiles) <- labs - if (!is.null(gffFile)) { + txi_files <- as.character(sample_table$filename) + labs <- row.names(sample_table) + names(txi_files) <- labs + if (!is.null(gff_file)) { # first need to make the tx2gene table # this takes ~2-3 minutes using Bioconductor functions suppressPackageStartupMessages({ library("GenomicFeatures") }) - txdb <- makeTxDbFromGFF(gffFile) + txdb <- makeTxDbFromGFF(gff_file) k <- keys(txdb, keytype = "TXNAME") - tx2gene <- select(txdb, keys=k, columns="GENEID", keytype="TXNAME") - # Remove 'transcript:' from transcript IDs (when gffFile is a GFF3 from Ensembl and the transcript does not have a Name) - tx2gene$TXNAME <- sub('^transcript:', '', tx2gene$TXNAME) + tx2gene <- select(txdb, keys = k, columns = "GENEID", keytype = "TXNAME") + # Remove 'transcript:' from transcript IDs (when gff_file is a GFF3 from Ensembl and the transcript does not have a Name) + tx2gene$TXNAME <- sub("^transcript:", "", tx2gene$TXNAME) # nolint } - try(txi <- tximport(txiFiles, type=txtype, tx2gene=tx2gene)) + try(txi <- tximport(txi_files, type = txtype, tx2gene = tx2gene)) if (!exists("txi")) { # Remove version from transcript IDs in tx2gene... - tx2gene$TXNAME <- sub('\\.[0-9]+$', '', tx2gene$TXNAME) - # ...and in txiFiles - txi <- tximport(txiFiles, type=txtype, tx2gene=tx2gene, ignoreTxVersion=TRUE) + tx2gene$TXNAME <- sub("\\.[0-9]+$", "", tx2gene$TXNAME) # nolint + # ...and in txi_files + txi <- tximport(txi_files, type = txtype, tx2gene = tx2gene, ignoreTxVersion = TRUE) } - dds <- DESeqDataSetFromTximport(txi, - subset(sampleTable, select=-c(filename)), - designFormula) + dds <- DESeqDataSetFromTximport( + txi, + subset(sample_table, select = -c(filename)), + design_formula + ) } return(dds) } |
| b |
| diff -r c24765926774 -r fed9d0350d72 ruvseq.R --- a/ruvseq.R Tue Mar 26 06:25:38 2019 -0400 +++ b/ruvseq.R Fri Jan 15 17:53:15 2021 +0000 |
| [ |
| b'@@ -1,6 +1,8 @@\n # setup R error handling to go to stderr\n library("getopt")\n-options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )\n+options(show.error.messages = F, error = function() {\n+ cat(geterrmessage(), file = stderr()); q("no", 1, F)\n+})\n options(stringAsFactors = FALSE, useFancyQuotes = FALSE)\n \n setup_cmdline_options <- function() {\n@@ -12,18 +14,18 @@\n "min_k", "min_k", 1, "double",\n "max_k", "max_k", 1, "double",\n "sample_json", "s", 1, "character",\n- "plots" , "p", 1, "character",\n+ "plots", "p", 1, "character",\n "header", "H", 0, "logical",\n "txtype", "y", 1, "character",\n "tx2gene", "x", 1, "character"), # a space-sep tx-to-gene map or GTF file (auto detect .gtf/.GTF)\n- byrow=TRUE, ncol=4)\n+ byrow = TRUE, ncol = 4)\n \n opt <- getopt(spec)\n # if help was asked for print a friendly message\n # and exit with a non-zero error code\n if (!is.null(opt$help)) {\n- cat(getopt(spec, usage=TRUE))\n- q(status=1)\n+ cat(getopt(spec, usage = TRUE))\n+ q(status = 1)\n } else {\n load_libraries()\n }\n@@ -42,117 +44,112 @@\n library("ggrepel")\n }\n \n-source_local <- function(fname){\n+source_local <- function(fname) {\n argv <- commandArgs(trailingOnly = FALSE)\n base_dir <- dirname(substring(argv[grep("--file=", argv)], 8))\n- source(paste(base_dir, fname, sep="/"))\n+ source(paste(base_dir, fname, sep = "/"))\n }\n \n # Source get_deseq_dataset.R for getting deseq dataset from htseq/featurecounts/tximport\n-source_local(\'get_deseq_dataset.R\')\n+source_local("get_deseq_dataset.R")\n \n # RUVseq function definitions\n \n-plot_pca_rle <- function (set, title) {\n- x <- pData(set)[,1]\n+plot_pca_rle <- function(set, title) {\n+ x <- pData(set)[, 1]\n colors <- brewer.pal(3, "Set2")\n- label <- paste0(\' for \', title)\n- plotRLE(set, outline=FALSE, ylim=c(-4, 4), col=colors[x])\n- title(main=paste0("RLE", label))\n- plotPCA(set, col=colors[x], cex=1.2)\n- title(main=paste0("PCA", label))\n+ label <- paste0(" for ", title)\n+ plotRLE(set, outline = FALSE, ylim = c(-4, 4), col = colors[x])\n+ title(main = paste0("RLE", label))\n+ plotPCA(set, col = colors[x], cex = 1.2)\n+ title(main = paste0("PCA", label))\n }\n \n-plot_factors_of_unwanted_variation <- function(set, method, k){\n+plot_factors_of_unwanted_var <- function(set, method, k) {\n pd <- pData(set)\n- pd[\'sample\'] <- row.names(pd)\n- colnames(pd)[1] <- \'condition\'\n- d = melt(pd, id.vars = c(\'sample\', \'condition\'))\n- d[\'x\'] <- 1 # There is no information on the X, so we just fake it to be able to do a scatterplot\n- print(ggplot(d, aes(x=x, y=value, color=condition, label=sample)) +\n+ pd["sample"] <- row.names(pd)\n+ colnames(pd)[1] <- "condition"\n+ d <- melt(pd, id.vars = c("sample", "condition"))\n+ d["x"] <- 1 # There is no information on the X, so we just fake it to be able to do a scatterplot\n+ print(ggplot(d, aes(x = x, y = value, color = condition, label = sample)) +\n geom_point() +\n- ggtitle(paste0(\'Factors of unwanted variation for method: \', method, ", k=", k)) +\n- facet_wrap( ~ variable, scales = "free_x") +\n+ ggtitle(paste0("Factors of unwanted variation for method: ", method, ", k=", k)) +\n+ facet_wrap(~ variable, scales = "free_x") +\n geom_text_repel() +\n- theme(axis.title.x=element_blank(),\n- axis.text.x=element_blank(),\n- axis.ticks.x=element_blank(),\n+ theme(axis.title.x = element_blank(),\n+ axis.text.x = element_blank(),\n+ axis.ticks.x = element_blank(),\n plot.title = element_text(hjust = 0.5))\n )\n }\n \n-create_seq_expression_set <- function (dds, min_mean_count) {\n+create_seq_expression_set <- function(dds, min_mean_count) {\n count_values <- counts(dds)\n- print(paste0("feature count before filtering :",nrow(count_values),"\\n"))\n+ print(paste0("feature count before filtering :", nrow(count_values), "\\n"))\n print(paste0("Filtering features which mean expression is less or eq. than ", min_me'..b'\n+ plot_factors_of_unwanted_var(set, method = "RUVg with empirical control genes", k = k)\n return(set)\n }\n \n ruv_residual_method <- function(set, k) {\n genes <- rownames(counts(set))\n- x <- pData(set)[,1]\n+ x <- pData(set)[, 1]\n # Initial edger residuals\n- design <- model.matrix(~x, data=pData(set))\n- y <- DGEList(counts=counts(set), group=x)\n- y <- calcNormFactors(y, method="upperquartile")\n+ design <- model.matrix(~x, data = pData(set))\n+ y <- DGEList(counts = counts(set), group = x)\n+ y <- calcNormFactors(y, method = "upperquartile")\n y <- estimateGLMCommonDisp(y, design)\n y <- estimateGLMTagwiseDisp(y, design)\n fit <- glmFit(y, design)\n- res <- residuals(fit, type="deviance")\n- set <- RUVr(set, genes, k=k, res)\n- plot_pca_rle(set = set, title = paste0(\'RUVr using residuals, k=\', k))\n- plot_factors_of_unwanted_variation(set, method="RUVr using residuals", k=k)\n+ res <- residuals(fit, type = "deviance")\n+ set <- RUVr(set, genes, k = k, res)\n+ plot_pca_rle(set = set, title = paste0("RUVr using residuals, k=", k))\n+ plot_factors_of_unwanted_var(set, method = "RUVr using residuals", k = k)\n return(set)\n }\n \n-ruv_replicate_method <- function (set, k) {\n+ruv_replicate_method <- function(set, k) {\n genes <- rownames(counts(set))\n- x <- pData(set)[,1]\n+ x <- pData(set)[, 1]\n differences <- makeGroups(x)\n- set <- RUVs(set, genes, k=k, differences)\n- plot_pca_rle(set, paste0(\'RUVs with replicate samples, k=\', k))\n- plot_factors_of_unwanted_variation(set, method="RUVs using replicates", k=k)\n+ set <- RUVs(set, genes, k = k, differences)\n+ plot_pca_rle(set, paste0("RUVs with replicate samples, k=", k))\n+ plot_factors_of_unwanted_var(set, method = "RUVs using replicates", k = k)\n return(set)\n }\n \n-get_differentially_expressed_genes <- function(dds, contrast, alpha=0.01) {\n- r <- results(dds, contrast=contrast, alpha=alpha)\n- return(rownames(r[which(r$padj < alpha),]))\n-}\n-\n opt <- setup_cmdline_options()\n alpha <- opt$alpha\n min_k <- opt$min_k\n@@ -162,12 +159,10 @@\n sample_paths <- sample_json$path\n sample_names <- sample_json$label\n condition <- as.factor(sample_json$condition)\n-sampleTable <- data.frame(samplename=sample_names,\n- filename = sample_paths,\n- condition=condition)\n-rownames(sampleTable) <- sample_names\n+sample_table <- data.frame(samplename = sample_names, filename = sample_paths, condition = condition)\n+rownames(sample_table) <- sample_names\n \n-dds <- get_deseq_dataset(sampleTable, header=opt$header, designFormula= ~ condition, tximport=opt$txtype, txtype=opt$txtype, tx2gene=opt$tx2gene)\n+dds <- get_deseq_dataset(sample_table, header = opt$header, design_formula = ~ condition, tximport = opt$txtype, txtype = opt$txtype, tx2gene = opt$tx2gene)\n if (!is.null(opt$plots)) {\n pdf(opt$plots)\n }\n@@ -176,9 +171,9 @@\n set <- create_seq_expression_set(dds, min_mean_count = min_c)\n result <- list(no_correction = set)\n for (k in seq(min_k, max_k)) {\n- result[[paste0(\'residual_method_k\', k)]] <- ruv_residual_method(set, k=k)\n- result[[paste0(\'replicate_method_k\', k)]] <- ruv_replicate_method(set, k=k)\n- result[[paste0(\'control_method_k\', k)]] <- ruv_control_gene_method(set, k=k, cutoff_p=0.5)\n+ result[[paste0("residual_method_k", k)]] <- ruv_residual_method(set, k = k)\n+ result[[paste0("replicate_method_k", k)]] <- ruv_replicate_method(set, k = k)\n+ result[[paste0("control_method_k", k)]] <- ruv_control_gene_method(set, k = k, cutoff_p = 0.5)\n }\n \n for (name in names(result)) {\n@@ -187,8 +182,8 @@\n unwanted_variation <- pData(set)\n df <- data.frame(identifier = rownames(unwanted_variation))\n df <- cbind(df, unwanted_variation)\n- colnames(df)[2] <- \'condition\'\n- write.table(df, file=paste0("batch_effects_", name, ".tabular"), sep="\\t", quote=F, row.names=F)\n+ colnames(df)[2] <- "condition"\n+ write.table(df, file = paste0("batch_effects_", name, ".tabular"), sep = "\\t", quote = F, row.names = F)\n }\n }\n \n' |
| b |
| diff -r c24765926774 -r fed9d0350d72 ruvseq.xml --- a/ruvseq.xml Tue Mar 26 06:25:38 2019 -0400 +++ b/ruvseq.xml Fri Jan 15 17:53:15 2021 +0000 |
| b |
| @@ -1,7 +1,11 @@ -<tool id="ruvseq" name="Remove Unwanted Variation" version="1.16.0"> +<tool id="ruvseq" name="Remove Unwanted Variation" version="@TOOL_VERSION@+galaxy@WRAPPER_VERSION@"> <description>from RNA-seq data</description> + <macros> + <token name="@TOOL_VERSION@">1.16.0</token> + <token name="@WRAPPER_VERSION@">1</token> + </macros> <requirements> - <requirement type="package" version="1.16.0">bioconductor-ruvseq</requirement> + <requirement type="package" version="@TOOL_VERSION@">bioconductor-ruvseq</requirement> <requirement type="package" version="1.22.1">bioconductor-deseq2</requirement> <requirement type="package" version="1.10.0">bioconductor-tximport</requirement> <requirement type="package" version="1.34.1">bioconductor-genomicfeatures</requirement> |