Mercurial > repos > iuc > ruvseq

comparison get_deseq_dataset.R @ 1:c24765926774 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/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
author iuc
date Tue, 26 Mar 2019 06:25:38 -0400
parents 61dffb20b6f9
children fed9d0350d72
comparison
equal deleted inserted replaced
0:61dffb20b6f9 1:c24765926774
7 } else { 7 } else {
8 hasHeader <- FALSE 8 hasHeader <- FALSE
9 } 9 }
10 10
11 if (!is.null(tximport)) { 11 if (!is.null(tximport)) {
12 if (is.null(tx2gene)) stop("A transcript-to-gene map or a GTF file is required for tximport") 12 if (is.null(tx2gene)) stop("A transcript-to-gene map or a GTF/GFF3 file is required for tximport")
13 if (tolower(file_ext(opt$tx2gene)) == "gtf") { 13 if (tolower(file_ext(tx2gene)) == "gff") {
14 gtfFile <-tx2gene 14 gffFile <-tx2gene
15 } else { 15 } else {
16 gtfFile <- NULL 16 gffFile <- NULL
17 tx2gene <- read.table(tx2gene, header=FALSE) 17 tx2gene <- read.table(tx2gene, header=hasHeader)
18 } 18 }
19 useTXI <- TRUE 19 useTXI <- TRUE
20 } else { 20 } else {
21 useTXI <- FALSE 21 useTXI <- FALSE
22 } 22 }
43 design = designFormula) 43 design = designFormula)
44 colnames(dds) <- row.names(sampleTable) 44 colnames(dds) <- row.names(sampleTable)
45 45
46 } else { 46 } else {
47 # construct the object using tximport 47 # construct the object using tximport
48 # first need to make the tx2gene table
49 # this takes ~2-3 minutes using Bioconductor functions
50 if (!is.null(gtfFile)) {
51 suppressPackageStartupMessages({
52 library("GenomicFeatures")
53 })
54 txdb <- makeTxDbFromGFF(gtfFile, format="gtf")
55 k <- keys(txdb, keytype = "GENEID")
56 df <- select(txdb, keys = k, keytype = "GENEID", columns = "TXNAME")
57 tx2gene <- df[, 2:1] # tx ID, then gene ID
58 }
59 library("tximport") 48 library("tximport")
60 txiFiles <- as.character(sampleTable$filename) 49 txiFiles <- as.character(sampleTable$filename)
61 labs <- row.names(sampleTable) 50 labs <- row.names(sampleTable)
62 names(txiFiles) <- labs 51 names(txiFiles) <- labs
63 txi <- tximport(txiFiles, type=txtype, tx2gene=tx2gene) 52 if (!is.null(gffFile)) {
53 # first need to make the tx2gene table
54 # this takes ~2-3 minutes using Bioconductor functions
55 suppressPackageStartupMessages({
56 library("GenomicFeatures")
57 })
58 txdb <- makeTxDbFromGFF(gffFile)
59 k <- keys(txdb, keytype = "TXNAME")
60 tx2gene <- select(txdb, keys=k, columns="GENEID", keytype="TXNAME")
61 # Remove 'transcript:' from transcript IDs (when gffFile is a GFF3 from Ensembl and the transcript does not have a Name)
62 tx2gene$TXNAME <- sub('^transcript:', '', tx2gene$TXNAME)
63 }
64 try(txi <- tximport(txiFiles, type=txtype, tx2gene=tx2gene))
65 if (!exists("txi")) {
66 # Remove version from transcript IDs in tx2gene...
67 tx2gene$TXNAME <- sub('\\.[0-9]+$', '', tx2gene$TXNAME)
68 # ...and in txiFiles
69 txi <- tximport(txiFiles, type=txtype, tx2gene=tx2gene, ignoreTxVersion=TRUE)
70 }
64 dds <- DESeqDataSetFromTximport(txi, 71 dds <- DESeqDataSetFromTximport(txi,
65 subset(sampleTable, select=-c(filename)), 72 subset(sampleTable, select=-c(filename)),
66 designFormula) 73 designFormula)
67 } 74 }
68 return(dds) 75 return(dds)