annotate get_deseq_dataset.R @ 5:557362ddaf9e draft default tip

planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 30117fce22f3771c9c0c13e70c3ad14b694de6e2
author iuc
date Fri, 21 Apr 2023 14:08:58 +0000
parents a32619eef8b9
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
1 get_deseq_dataset <- function(sample_table, header, design_formula, tximport, txtype, tx2gene) {
0
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
2
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
3 dir <- ""
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
4
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
5 has_header <- !is.null(header)
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
6 use_txi <- !is.null(tximport)
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
7 if (use_txi) {
1
99e0c629d3a2 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
iuc
parents: 0
diff changeset
8 if (is.null(tx2gene)) stop("A transcript-to-gene map or a GTF/GFF3 file is required for tximport")
99e0c629d3a2 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
iuc
parents: 0
diff changeset
9 if (tolower(file_ext(tx2gene)) == "gff") {
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
10 gff_file <- tx2gene
0
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
11 } else {
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
12 gff_file <- NULL
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
13 tx2gene <- read.table(tx2gene, header = has_header)
0
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
14 }
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
15 }
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
16
5
557362ddaf9e planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 30117fce22f3771c9c0c13e70c3ad14b694de6e2
iuc
parents: 2
diff changeset
17 if (!use_txi && has_header) {
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
18 countfiles <- lapply(as.character(sample_table$filename), read.delim, row.names = 1)
0
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
19 tbl <- do.call("cbind", countfiles)
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
20 colnames(tbl) <- rownames(sample_table) # take sample ids from header
0
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
21
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
22 # check for htseq report lines (from DESeqDataSetFromHTSeqCount function)
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
23 old_special_names <- c(
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
24 "no_feature",
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
25 "ambiguous",
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
26 "too_low_aQual",
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
27 "not_aligned",
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
28 "alignment_not_unique"
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
29 )
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
30 special_rows <- (substr(rownames(tbl), 1, 1) == "_") | rownames(tbl) %in% old_special_names
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
31 tbl <- tbl[!special_rows, , drop = FALSE]
0
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
32
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
33 dds <- DESeqDataSetFromMatrix(
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
34 countData = tbl,
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
35 colData = subset(sample_table, select = -filename),
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
36 design = design_formula
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
37 )
5
557362ddaf9e planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 30117fce22f3771c9c0c13e70c3ad14b694de6e2
iuc
parents: 2
diff changeset
38 } else if (!use_txi && !has_header) {
0
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
39
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
40 # construct the object from HTSeq files
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
41 dds <- DESeqDataSetFromHTSeqCount(
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
42 sampleTable = sample_table,
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
43 directory = dir,
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
44 design = design_formula
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
45 )
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
46 colnames(dds) <- row.names(sample_table)
0
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
47
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
48 } else {
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
49 # construct the object using tximport
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
50 library("tximport")
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
51 txi_files <- as.character(sample_table$filename)
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
52 labs <- row.names(sample_table)
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
53 names(txi_files) <- labs
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
54 if (!is.null(gff_file)) {
1
99e0c629d3a2 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
iuc
parents: 0
diff changeset
55 # first need to make the tx2gene table
99e0c629d3a2 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
iuc
parents: 0
diff changeset
56 # this takes ~2-3 minutes using Bioconductor functions
99e0c629d3a2 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
iuc
parents: 0
diff changeset
57 suppressPackageStartupMessages({
99e0c629d3a2 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
iuc
parents: 0
diff changeset
58 library("GenomicFeatures")
99e0c629d3a2 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
iuc
parents: 0
diff changeset
59 })
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
60 txdb <- makeTxDbFromGFF(gff_file)
1
99e0c629d3a2 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
iuc
parents: 0
diff changeset
61 k <- keys(txdb, keytype = "TXNAME")
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
62 tx2gene <- select(txdb, keys = k, columns = "GENEID", keytype = "TXNAME")
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
63 # Remove 'transcript:' from transcript IDs (when gff_file is a GFF3 from Ensembl and the transcript does not have a Name)
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
64 tx2gene$TXNAME <- sub("^transcript:", "", tx2gene$TXNAME) # nolint
1
99e0c629d3a2 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
iuc
parents: 0
diff changeset
65 }
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
66 try(txi <- tximport(txi_files, type = txtype, tx2gene = tx2gene))
1
99e0c629d3a2 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
iuc
parents: 0
diff changeset
67 if (!exists("txi")) {
99e0c629d3a2 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
iuc
parents: 0
diff changeset
68 # Remove version from transcript IDs in tx2gene...
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
69 tx2gene$TXNAME <- sub("\\.[0-9]+$", "", tx2gene$TXNAME) # nolint
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
70 # ...and in txi_files
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
71 txi <- tximport(txi_files, type = txtype, tx2gene = tx2gene, ignoreTxVersion = TRUE)
1
99e0c629d3a2 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 9ed3d83cc447ee897af867362bf1dd67af8a11c2
iuc
parents: 0
diff changeset
72 }
2
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
73 dds <- DESeqDataSetFromTximport(
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
74 txi,
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
75 subset(sample_table, select = -c(filename)),
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
76 design_formula
a32619eef8b9 "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit 4daa375d022673d2437d609b1865b78c64b04415"
iuc
parents: 1
diff changeset
77 )
0
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
78 }
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
79 return(dds)
5f56d460bce3 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ruvseq commit b95582cea8320d5488056a9576474f79cec53be8
iuc
parents:
diff changeset
80 }