Mercurial > repos > youyuh48 > sleuth
diff sleuth.R @ 3:172091d47f5a draft default tip
planemo upload for repository https://github.com/youyuh48/galaxy-sleuth/tree/master/ commit 8d9ff6181e933951e015e67f1719588bab2b3ab5-dirty
author | youyuh48 |
---|---|
date | Tue, 12 Jun 2018 17:58:47 -0400 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sleuth.R Tue Jun 12 17:58:47 2018 -0400 @@ -0,0 +1,101 @@ +#!/usr/bin/env Rscript + +# A command-line interface to sleuth for use with Galaxy This script modified +# from https://github.com/pachterlab/bears_analyses/blob/master/sleuth.R +# https://github.com/nturaga/bioc-galaxy-integration/blob/master/README.md + +## Command to run tool: +## Rscript sleuth.R --indir test-rscript --metadata test-rscript/metadata.txt +## --full_model '~condition' --reduced_model '~1' +## --gene_anno_name 'hsapiens_gene_ensembl' + +# setup R error handling to go to stderr +options(show.error.messages = F, error = function() { + cat(geterrmessage(), file = stderr()) + q("no", 1, F) +}) + +# we need that to not crash galaxy with an UTF8 error on German LC settings. +loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") + +library("getopt") +library("tools") +options(stringAsFactors = FALSE, useFancyQuotes = FALSE) +args <- commandArgs(trailingOnly = TRUE) + +# get options, using the spec as defined by the enclosed list. +spec <- matrix(c( + "quiet", "q", 0, "logical", + "help", "h", 0, "logical", + "indir", "i", 1, "character", + "metadata", "m", 1, "character", + "full_model", "f", 1, "character", + "reduced_model", "r", 1, "character", + "gene_anno_name", "a", 2, "character"), + byrow = TRUE, ncol = 4) +opt <- getopt(spec) + +# if help was asked for print a friendly message and exit with a non-zero error +# code +if (!is.null(opt$help)) { + cat(getopt(spec, usage = TRUE)) + q(status = 1) +} + +# enforce the following required arguments +if (is.null(opt$indir)) { + cat("'indir' is required\n") + q(status = 1) +} +if (is.null(opt$metadata)) { + cat("'metadata' is required\n") + q(status = 1) +} +if (is.null(opt$full_model)) { + cat("'full_model' is required\n") + q(status = 1) +} +if (is.null(opt$reduced_model)) { + cat("'reduced_model' is required\n") + q(status = 1) +} + +verbose <- if (is.null(opt$quiet)) { + TRUE +} else { + FALSE +} + +suppressPackageStartupMessages({ + library("sleuth") + library("biomaRt") +}) + +s2c <- read.table(file.path(opt$metadata), header = TRUE, stringsAsFactors = FALSE) +run_dirs <- s2c$sample +kal_dirs <- c() + +for (dir in run_dirs) { + kal_dirs <- c(kal_dirs, file.path(opt$indir, dir, "kallisto")) +} + +s2c <- dplyr::mutate(s2c, path = kal_dirs) + +if (!is.null(opt$gene_anno_name)) { + mart <- biomaRt::useMart(biomart = "ensembl", dataset = opt$gene_anno_name) + t2g <- biomaRt::getBM(attributes = c("ensembl_transcript_id", "ensembl_gene_id", + "external_gene_name"), mart = mart) + t2g <- dplyr::rename(t2g, target_id = ensembl_transcript_id, ens_gene = ensembl_gene_id, + ext_gene = external_gene_name) + so <- sleuth_prep(s2c, as.formula(opt$full_model), target_mapping = t2g, read_bootstrap_tpm = TRUE, + extra_bootstrap_summary = TRUE) +} else { + so <- sleuth_prep(s2c, as.formula(opt$full_model), read_bootstrap_tpm = TRUE, + extra_bootstrap_summary = TRUE) +} +so <- sleuth_fit(so, as.formula(opt$full_model), "full") +so <- sleuth_fit(so, as.formula(opt$reduced_model), "reduced") +so <- sleuth_lrt(so, "reduced", "full") +sleuth_deploy(so, opt$indir) + +cat("Successfully finished script.\n")