Mercurial > repos > iuc > tximport

comparison tximport.R @ 0:d53de4a14836 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/tximport commit 66c3b86403faa115751332ea8e0383e26b9ee599"
author iuc
date Wed, 04 Dec 2019 10:56:58 +0000
parents
children 3e2f1aef1aac
comparison
equal deleted inserted replaced
-1:000000000000 0:d53de4a14836
1 # setup R error handling to go to stderr
2 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
3
4 # we need that to not crash galaxy with an UTF8 error on German LC settings.
5 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
6
7 library("getopt")
8 options(stringAsFactors = FALSE, useFancyQuotes = FALSE)
9 args <- commandArgs(trailingOnly = TRUE)
10
11 # get options, using the spec as defined by the enclosed list.
12 # we read the options from the default: commandArgs(TRUE).
13 spec <- matrix(c(
14 "help", "h", 0, "logical",
15 "out_file", "o", 1, "character",
16 "countsFiles", "n", 1, "character",
17 "countsFromAbundance", "r", 1, "character",
18 "format", "v", 1, "character",
19 "gff_file", "H", 0, "character",
20 "tx2gene", "f", 0, "character",
21 "geneIdCol", "l", 0, "character",
22 "txIdCol" , "p", 1, "character",
23 "abundanceCol", "i", 0, "character",
24 "countsCol", "y", 1, "character",
25 "lengthCol", "x", 1, "character"),
26 byrow=TRUE, ncol=4)
27
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 if (is.null(opt$gff_file) & is.null(opt$tx2gene)) {
38 cat("A GFF/GTF file or a tx2gene table is required\n")
39 q(status=1)
40 }
41
42 if (opt$format == 'none'){ #custom format
43 if (is.null(opt$txIdCol) | is.null(opt$abundanceCol) | is.null(opt$countsCol) | is.null(opt$lengthCol)) {
44 cat("If you select a custom format for the input files you need to specify the column names\n")
45 q(status=1)
46 }
47 }
48
49 if (is.null(opt$countsFiles)) {
50 cat("'countsFiles' is required\n")
51 q(status=1)
52 }
53
54
55 # load samples from tab file
56 samples_df <- read.table(opt$countsFiles, sep="\t", header=TRUE)
57 colnames(samples_df) <- c("id","path")
58 rownames(samples_df) <- NULL
59 # Prepare char vector with files and sample names
60 files <- file.path(samples_df[,"path"])
61 names(files) <- samples_df[,"id"]
62
63
64
65 library(tximport)
66
67 ### if the input is a gff/gtf file first need to create the tx2gene table
68 if (!is.null(opt$gff_file)) {
69 suppressPackageStartupMessages({
70 library("GenomicFeatures")
71 })
72 txdb <- makeTxDbFromGFF(opt$gff_file)
73 k <- keys(txdb, keytype = "TXNAME")
74 tx2gene <- select(txdb, keys=k, columns="GENEID", keytype="TXNAME")
75 # Remove 'transcript:' from transcript IDs (when gffFile is a GFF3 from Ensembl and the transcript does not have a Name)
76 tx2gene$TXNAME <- sub('^transcript:', '', tx2gene$TXNAME)
77
78 } else {
79 tx2gene <- read.table(opt$tx2gene,header=FALSE)
80 }
81
82
83
84 ##
85 if (is.null(opt$geneIdCol)) { ## there is a tx2gene table
86 if (opt$format == 'none'){ #predefined format
87 txi_out <- tximport(files, type="none",txIdCol=opt$txIdCol,abundanceCol=opt$abundanceCol,countsCol=opt$countsCol,lengthCol=opt$lengthCol,tx2gene=tx2gene,countsFromAbundance=opt$countsFromAbundance)
88 } else {
89 txi_out <- tximport(files, type=opt$format, tx2gene=tx2gene,countsFromAbundance=opt$countsFromAbundance)
90 }
91 } else { # the gene_ID is a column in the counts table
92 if (opt$format == 'none'){ #predefined format
93 txi_out <- tximport(files, type="none",geneIdCol=opt$geneIdCol,txIdCol=opt$txIdCol,abundanceCol=opt$abundanceCol,countsCol=opt$countsCol,lengthCol=opt$lengthCol,tx2gene=tx2gene,countsFromAbundance=opt$countsFromAbundance)
94 } else {
95 txi_out <- tximport(files, type=opt$format, geneIdCol=opt$geneIdCol,countsFromAbundance=opt$countsFromAbundance)
96 }
97
98 }
99 # write count as table
100 write.table(txi_out$counts, file=opt$out_file, row.names = TRUE, col.names = TRUE, quote = FALSE, sep = "\t")