|
2
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 # A command-line interface to sleuth for use with Galaxy
|
|
|
4 # This script modified from
|
|
|
5 # https://github.com/pachterlab/bears_analyses/blob/master/sleuth.R
|
|
|
6
|
|
|
7 # setup R error handling to go to stderr
|
|
|
8 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
|
|
|
9
|
|
|
10 # we need that to not crash galaxy with an UTF8 error on German LC settings.
|
|
|
11 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
|
|
|
12
|
|
|
13 library("getopt")
|
|
|
14 library("tools")
|
|
|
15 options(stringAsFactors = FALSE, useFancyQuotes = FALSE)
|
|
|
16 args <- commandArgs(trailingOnly = TRUE)
|
|
|
17
|
|
|
18 # get options, using the spec as defined by the enclosed list.
|
|
|
19 spec <- matrix(c(
|
|
|
20 "quiet", "q", 0, "logical",
|
|
|
21 "help", "h", 0, "logical",
|
|
|
22 "indir", "i", 1, "character",
|
|
|
23 "metadata", "m", 1, "character",
|
|
|
24 "full_model", "f", 1, "character",
|
|
|
25 "reduced_model", "r", 1, "character",
|
|
|
26 "gene_anno_name", "a", 2, "character"
|
|
|
27 ), byrow=TRUE, ncol=4)
|
|
|
28 opt <- getopt(spec)
|
|
|
29
|
|
|
30 # if help was asked for print a friendly message
|
|
|
31 # and exit with a non-zero error code
|
|
|
32 if (!is.null(opt$help)) {
|
|
|
33 cat(getopt(spec, usage=TRUE))
|
|
|
34 q(status=1)
|
|
|
35 }
|
|
|
36
|
|
|
37 # enforce the following required arguments
|
|
|
38 if (is.null(opt$indir)) {
|
|
|
39 cat("'indir' is required\n")
|
|
|
40 q(status=1)
|
|
|
41 }
|
|
|
42 if (is.null(opt$metadata)) {
|
|
|
43 cat("'metadata' is required\n")
|
|
|
44 q(status=1)
|
|
|
45 }
|
|
|
46 if (is.null(opt$full_model)) {
|
|
|
47 cat("'full_model' is required\n")
|
|
|
48 q(status=1)
|
|
|
49 }
|
|
|
50 if (is.null(opt$reduced_model)) {
|
|
|
51 cat("'reduced_model' is required\n")
|
|
|
52 q(status=1)
|
|
|
53 }
|
|
|
54
|
|
|
55 verbose <- if (is.null(opt$quiet)) {
|
|
|
56 TRUE
|
|
|
57 } else {
|
|
|
58 FALSE
|
|
|
59 }
|
|
|
60
|
|
|
61 suppressPackageStartupMessages({
|
|
|
62 library("sleuth")
|
|
|
63 library("biomaRt")
|
|
|
64 })
|
|
|
65
|
|
|
66 s2c <- read.table(file.path(opt$metadata), header = TRUE, stringsAsFactors=FALSE)
|
|
|
67 run_dirs <- s2c$sample
|
|
|
68 kal_dirs <- c()
|
|
|
69
|
|
|
70 for (dir in run_dirs) {
|
|
|
71 kal_dirs <- c(kal_dirs, file.path(opt$indir, dir, "kallisto"))
|
|
|
72 }
|
|
|
73
|
|
|
74 s2c <- dplyr::mutate(s2c, path = kal_dirs)
|
|
|
75
|
|
|
76 if (!is.null(opt$gene_anno_name)) {
|
|
|
77 mart <- biomaRt::useMart(biomart = "ensembl", dataset = opt$gene_anno_name)
|
|
|
78 t2g <- biomaRt::getBM(attributes = c("ensembl_transcript_id", "ensembl_gene_id", "external_gene_name"), mart = mart)
|
|
|
79 t2g <- dplyr::rename(t2g, target_id = ensembl_transcript_id, ens_gene = ensembl_gene_id, ext_gene = external_gene_name)
|
|
|
80 so <- sleuth_prep(s2c, as.formula(opt$full_model), target_mapping = t2g, read_bootstrap_tpm=TRUE, extra_bootstrap_summary=TRUE)
|
|
|
81 } else {
|
|
|
82 so <- sleuth_prep(s2c, as.formula(opt$full_model), read_bootstrap_tpm=TRUE, extra_bootstrap_summary=TRUE)
|
|
|
83 }
|
|
|
84 so <- sleuth_fit(so, as.formula(opt$full_model), "full")
|
|
|
85 so <- sleuth_fit(so, as.formula(opt$reduced_model), "reduced")
|
|
|
86 so <- sleuth_lrt(so, "reduced", "full")
|
|
|
87 sleuth_deploy(so, opt$indir)
|
