Mercurial > repos > iuc > edger

comparison edger.R @ 16:ae2aad0a6d50 draft default tip

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/edger commit 4b17306763bc8eb92c8967c7b4b57abcc514e670
author iuc
date Wed, 04 Sep 2024 15:50:01 +0000
parents 5bf899c13979
children
comparison
equal deleted inserted replaced
15:5bf899c13979 16:ae2aad0a6d50
40 # Record starting time 40 # Record starting time
41 time_start <- as.character(Sys.time()) 41 time_start <- as.character(Sys.time())
42 42
43 # setup R error handling to go to stderr 43 # setup R error handling to go to stderr
44 options(show.error.messages = FALSE, error = function() { 44 options(show.error.messages = FALSE, error = function() {
45 cat(geterrmessage(), file = stderr()) 45 cat(geterrmessage(), file = stderr())
46 q("no", 1, FALSE) 46 q("no", 1, FALSE)
47 }) 47 })
48 48
49 # we need that to not crash galaxy with an UTF8 error on German LC settings. 49 # we need that to not crash galaxy with an UTF8 error on German LC settings.
50 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") 50 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
51 51
62 ### Function Delcaration 62 ### Function Delcaration
63 ################################################################################ 63 ################################################################################
64 # Function to sanitise contrast equations so there are no whitespaces 64 # Function to sanitise contrast equations so there are no whitespaces
65 # surrounding the arithmetic operators, leading or trailing whitespace 65 # surrounding the arithmetic operators, leading or trailing whitespace
66 sanitise_equation <- function(equation) { 66 sanitise_equation <- function(equation) {
67 equation <- gsub(" *[+] *", "+", equation) 67 equation <- gsub(" *[+] *", "+", equation)
68 equation <- gsub(" *[-] *", "-", equation) 68 equation <- gsub(" *[-] *", "-", equation)
69 equation <- gsub(" *[/] *", "/", equation) 69 equation <- gsub(" *[/] *", "/", equation)
70 equation <- gsub(" *[*] *", "*", equation) 70 equation <- gsub(" *[*] *", "*", equation)
71 equation <- gsub("^\\s+|\\s+$", "", equation) 71 equation <- gsub("^\\s+|\\s+$", "", equation)
72 return(equation) 72 return(equation)
73 } 73 }
74 74
75 # Function to sanitise group information 75 # Function to sanitise group information
76 sanitise_groups <- function(string) { 76 sanitise_groups <- function(string) {
77 string <- gsub(" *[,] *", ",", string) 77 string <- gsub(" *[,] *", ",", string)
78 string <- gsub("^\\s+|\\s+$", "", string) 78 string <- gsub("^\\s+|\\s+$", "", string)
79 return(string) 79 return(string)
80 } 80 }
81 81
82 # Function to change periods to whitespace in a string 82 # Function to change periods to whitespace in a string
83 unmake_names <- function(string) { 83 unmake_names <- function(string) {
84 string <- gsub(".", " ", string, fixed = TRUE) 84 string <- gsub(".", " ", string, fixed = TRUE)
85 return(string) 85 return(string)
86 } 86 }
87 87
88 # Sanitise file base names coming from factors or contrasts 88 # Sanitise file base names coming from factors or contrasts
89 sanitise_basename <- function(string) { 89 sanitise_basename <- function(string) {
90 string <- gsub("[/^]", "_", string) 90 string <- gsub("[/^]", "_", string)
91 return(string) 91 return(string)
92 } 92 }
93 93
94 # Generate output folder and paths 94 # Generate output folder and paths
95 make_out <- function(filename) { 95 make_out <- function(filename) {
96 return(paste0(out_path, "/", filename)) 96 return(paste0(out_path, "/", filename))
97 } 97 }
98 98
99 # Generating design information 99 # Generating design information
100 paste_listname <- function(string) { 100 paste_listname <- function(string) {
101 return(paste0("factors$", string)) 101 return(paste0("factors$", string))
102 } 102 }
103 103
104 # Create cata function: default path set, default seperator empty and appending 104 # Create cata function: default path set, default seperator empty and appending
105 # true by default (Ripped straight from the cat function with altered argument 105 # true by default (Ripped straight from the cat function with altered argument
106 # defaults) 106 # defaults)
107 cata <- function(..., file = opt$htmlPath, sep = "", fill = FALSE, labels = NULL, 107 cata <- function(..., file = opt$htmlPath, sep = "", fill = FALSE, labels = NULL,
108 append = TRUE) { 108 append = TRUE) {
109 if (is.character(file)) { 109 if (is.character(file)) {
110 if (file == "") { 110 if (file == "") {
111 file <- stdout() 111 file <- stdout()
112 } else if (substring(file, 1L, 1L) == "|") { 112 } else if (substring(file, 1L, 1L) == "|") {
113 file <- pipe(substring(file, 2L), "w") 113 file <- pipe(substring(file, 2L), "w")
114 on.exit(close(file)) 114 on.exit(close(file))
115 } else { 115 } else {
116 file <- file(file, ifelse(append, "a", "w")) 116 file <- file(file, ifelse(append, "a", "w"))
117 on.exit(close(file)) 117 on.exit(close(file))
118 } 118 }
119 } 119 }
120 .Internal(cat(list(...), file, sep, fill, labels, append)) 120 .Internal(cat(list(...), file, sep, fill, labels, append))
121 } 121 }
122 122
123 # Function to write code for html head and title 123 # Function to write code for html head and title
124 html_head <- function(title) { 124 html_head <- function(title) {
125 cata("<head>\n") 125 cata("<head>\n")
126 cata("<title>", title, "</title>\n") 126 cata("<title>", title, "</title>\n")
127 cata("</head>\n") 127 cata("</head>\n")
128 } 128 }
129 129
130 # Function to write code for html links 130 # Function to write code for html links
131 html_link <- function(address, label = address) { 131 html_link <- function(address, label = address) {
132 cata("<a href=\"", address, "\" target=\"_blank\">", label, "</a><br />\n") 132 cata("<a href=\"", address, "\" target=\"_blank\">", label, "</a><br />\n")
133 } 133 }
134 134
135 # Function to write code for html images 135 # Function to write code for html images
136 html_image <- function(source, label = source, height = 600, width = 600) { 136 html_image <- function(source, label = source, height = 600, width = 600) {
137 cata("<img src=\"", source, "\" alt=\"", label, "\" height=\"", height) 137 cata("<img src=\"", source, "\" alt=\"", label, "\" height=\"", height)
138 cata("\" width=\"", width, "\"/>\n") 138 cata("\" width=\"", width, "\"/>\n")
139 } 139 }
140 140
141 # Function to write code for html list items 141 # Function to write code for html list items
142 list_item <- function(...) { 142 list_item <- function(...) {
143 cata("<li>", ..., "</li>\n") 143 cata("<li>", ..., "</li>\n")
144 } 144 }
145 145
146 table_item <- function(...) { 146 table_item <- function(...) {
147 cata("<td>", ..., "</td>\n") 147 cata("<td>", ..., "</td>\n")
148 } 148 }
149 149
150 table_head_item <- function(...) { 150 table_head_item <- function(...) {
151 cata("<th>", ..., "</th>\n") 151 cata("<th>", ..., "</th>\n")
152 } 152 }
153 153
154 ################################################################################ 154 ################################################################################
155 ### Input Processing 155 ### Input Processing
156 ################################################################################ 156 ################################################################################
158 # Collect arguments from command line 158 # Collect arguments from command line
159 args <- commandArgs(trailingOnly = TRUE) 159 args <- commandArgs(trailingOnly = TRUE)
160 160
161 # Get options, using the spec as defined by the enclosed list. 161 # Get options, using the spec as defined by the enclosed list.
162 # Read the options from the default: commandArgs(TRUE). 162 # Read the options from the default: commandArgs(TRUE).
163 spec <- matrix(c( 163 spec <- matrix(
164 "htmlPath", "R", 1, "character", 164 c(
165 "outPath", "o", 1, "character", 165 "htmlPath", "R", 1, "character",
166 "filesPath", "j", 2, "character", 166 "outPath", "o", 1, "character",
167 "matrixPath", "m", 2, "character", 167 "filesPath", "j", 2, "character",
168 "factFile", "f", 2, "character", 168 "matrixPath", "m", 2, "character",
169 "formula", "F", 2, "character", 169 "factFile", "f", 2, "character",
170 "factInput", "i", 2, "character", 170 "formula", "F", 2, "character",
171 "annoPath", "a", 2, "character", 171 "factInput", "i", 2, "character",
172 "contrastData", "C", 1, "character", 172 "annoPath", "a", 2, "character",
173 "cpmReq", "c", 1, "double", 173 "contrastData", "C", 1, "character",
174 "totReq", "y", 0, "logical", 174 "cpmReq", "c", 1, "double",
175 "cntReq", "z", 1, "integer", 175 "totReq", "y", 0, "logical",
176 "sampleReq", "s", 1, "integer", 176 "cntReq", "z", 1, "integer",
177 "normCounts", "x", 0, "logical", 177 "sampleReq", "s", 1, "integer",
178 "rdaOpt", "r", 0, "logical", 178 "normCounts", "x", 0, "logical",
179 "lfcReq", "l", 1, "double", 179 "rdaOpt", "r", 0, "logical",
180 "pValReq", "p", 1, "double", 180 "lfcReq", "l", 1, "double",
181 "pAdjOpt", "d", 1, "character", 181 "pValReq", "p", 1, "double",
182 "normOpt", "n", 1, "character", 182 "pAdjOpt", "d", 1, "character",
183 "robOpt", "b", 0, "logical", 183 "normOpt", "n", 1, "character",
184 "lrtOpt", "t", 0, "logical" 184 "robOpt", "b", 0, "logical",
185 ), 185 "lrtOpt", "t", 0, "logical"
186 byrow = TRUE, ncol = 4 186 ),
187 byrow = TRUE, ncol = 4
187 ) 188 )
188 opt <- getopt(spec) 189 opt <- getopt(spec)
189 190
190 191
191 if (is.null(opt$matrixPath) && is.null(opt$filesPath)) { 192 if (is.null(opt$matrixPath) && is.null(opt$filesPath)) {
192 cat("A counts matrix (or a set of counts files) is required.\n") 193 cat("A counts matrix (or a set of counts files) is required.\n")
193 q(status = 1) 194 q(status = 1)
194 } 195 }
195 196
196 if (is.null(opt$cpmReq)) { 197 if (is.null(opt$cpmReq)) {
197 filt_cpm <- FALSE 198 filt_cpm <- FALSE
198 } else { 199 } else {
199 filt_cpm <- TRUE 200 filt_cpm <- TRUE
200 } 201 }
201 202
202 if (is.null(opt$cntReq) || is.null(opt$sampleReq)) { 203 if (is.null(opt$cntReq) || is.null(opt$sampleReq)) {
203 filt_smpcount <- FALSE 204 filt_smpcount <- FALSE
204 } else { 205 } else {
205 filt_smpcount <- TRUE 206 filt_smpcount <- TRUE
206 } 207 }
207 208
208 if (is.null(opt$totReq)) { 209 if (is.null(opt$totReq)) {
209 filt_totcount <- FALSE 210 filt_totcount <- FALSE
210 } else { 211 } else {
211 filt_totcount <- TRUE 212 filt_totcount <- TRUE
212 } 213 }
213 214
214 if (is.null(opt$lrtOpt)) { 215 if (is.null(opt$lrtOpt)) {
215 want_lrt <- FALSE 216 want_lrt <- FALSE
216 } else { 217 } else {
217 want_lrt <- TRUE 218 want_lrt <- TRUE
218 } 219 }
219 220
220 if (is.null(opt$rdaOpt)) { 221 if (is.null(opt$rdaOpt)) {
221 want_rda <- FALSE 222 want_rda <- FALSE
222 } else { 223 } else {
223 want_rda <- TRUE 224 want_rda <- TRUE
224 } 225 }
225 226
226 if (is.null(opt$annoPath)) { 227 if (is.null(opt$annoPath)) {
227 have_anno <- FALSE 228 have_anno <- FALSE
228 } else { 229 } else {
229 have_anno <- TRUE 230 have_anno <- TRUE
230 } 231 }
231 232
232 if (is.null(opt$normCounts)) { 233 if (is.null(opt$normCounts)) {
233 want_norm <- FALSE 234 want_norm <- FALSE
234 } else { 235 } else {
235 want_norm <- TRUE 236 want_norm <- TRUE
236 } 237 }
237 238
238 if (is.null(opt$robOpt)) { 239 if (is.null(opt$robOpt)) {
239 want_robust <- FALSE 240 want_robust <- FALSE
240 } else { 241 } else {
241 want_robust <- TRUE 242 want_robust <- TRUE
242 } 243 }
243 244
244 245
245 if (!is.null(opt$filesPath)) { 246 if (!is.null(opt$filesPath)) {
246 # Process the separate count files (adapted from DESeq2 wrapper) 247 # Process the separate count files (adapted from DESeq2 wrapper)
247 library("rjson") 248 library("rjson")
248 parser <- newJSONParser() 249 parser <- newJSONParser()
249 parser$addData(opt$filesPath) 250 parser$addData(opt$filesPath)
250 factor_list <- parser$getObject() 251 factor_list <- parser$getObject()
251 factors <- sapply(factor_list, function(x) x[[1]]) 252 factors <- sapply(factor_list, function(x) x[[1]])
252 filenames_in <- unname(unlist(factor_list[[1]][[2]])) 253 filenames_in <- unname(unlist(factor_list[[1]][[2]]))
253 sampletable <- data.frame( 254 sampletable <- data.frame(
254 sample = basename(filenames_in), 255 sample = basename(filenames_in),
255 filename = filenames_in, 256 filename = filenames_in,
256 row.names = filenames_in, 257 row.names = filenames_in,
257 stringsAsFactors = FALSE 258 stringsAsFactors = FALSE
258 ) 259 )
259 for (factor in factor_list) { 260 for (factor in factor_list) {
260 factorname <- factor[[1]] 261 factorname <- factor[[1]]
261 sampletable[[factorname]] <- character(nrow(sampletable)) 262 sampletable[[factorname]] <- character(nrow(sampletable))
262 lvls <- sapply(factor[[2]], function(x) names(x)) 263 lvls <- sapply(factor[[2]], function(x) names(x))
263 for (i in seq_along(factor[[2]])) { 264 for (i in seq_along(factor[[2]])) {
264 files <- factor[[2]][[i]][[1]] 265 files <- factor[[2]][[i]][[1]]
265 sampletable[files, factorname] <- lvls[i] 266 sampletable[files, factorname] <- lvls[i]
266 } 267 }
267 sampletable[[factorname]] <- factor(sampletable[[factorname]], levels = lvls) 268 sampletable[[factorname]] <- factor(sampletable[[factorname]], levels = lvls)
268 } 269 }
269 rownames(sampletable) <- sampletable$sample 270 rownames(sampletable) <- sampletable$sample
270 rem <- c("sample", "filename") 271 rem <- c("sample", "filename")
271 factors <- sampletable[, !(names(sampletable) %in% rem), drop = FALSE] 272 factors <- sampletable[, !(names(sampletable) %in% rem), drop = FALSE]
272 273
273 # read in count files and create single table 274 # read in count files and create single table
274 countfiles <- lapply(sampletable$filename, function(x) { 275 countfiles <- lapply(sampletable$filename, function(x) {
275 read.delim(x, row.names = 1) 276 read.delim(x, row.names = 1)
276 }) 277 })
277 counts <- do.call("cbind", countfiles) 278 counts <- do.call("cbind", countfiles)
278 } else { 279 } else {
279 # Process the single count matrix 280 # Process the single count matrix
280 counts <- read.table(opt$matrixPath, header = TRUE, sep = "\t", strip.white = TRUE, stringsAsFactors = FALSE) 281 counts <- read.table(opt$matrixPath, header = TRUE, sep = "\t", strip.white = TRUE, stringsAsFactors = FALSE)
281 row.names(counts) <- counts[, 1] 282 row.names(counts) <- counts[, 1]
282 counts <- counts[, -1] 283 counts <- counts[, -1]
283 countsrows <- nrow(counts) 284 countsrows <- nrow(counts)
284 285
285 # Process factors 286 # Process factors
286 if (is.null(opt$factInput)) { 287 if (is.null(opt$factInput)) {
287 factordata <- read.table(opt$factFile, header = TRUE, sep = "\t", strip.white = TRUE, stringsAsFactors = TRUE) 288 factordata <- read.table(opt$factFile, header = TRUE, sep = "\t", strip.white = TRUE, stringsAsFactors = TRUE)
288 # check samples names match 289 # check samples names match
289 if (!any(factordata[, 1] %in% colnames(counts))) { 290 if (!any(factordata[, 1] %in% colnames(counts))) {
290 stop("Sample IDs in factors file and count matrix don't match") 291 stop("Sample IDs in factors file and count matrix don't match")
291 } 292 }
292 # order samples as in counts matrix 293 # order samples as in counts matrix
293 factordata <- factordata[match(colnames(counts), factordata[, 1]), ] 294 factordata <- factordata[match(colnames(counts), factordata[, 1]), ]
294 factors <- data.frame(sapply(factordata[, -1, drop = FALSE], make.names)) 295 factors <- data.frame(sapply(factordata[, -1, drop = FALSE], make.names))
295 } else { 296 } else {
296 factors <- unlist(strsplit(opt$factInput, "|", fixed = TRUE)) 297 factors <- unlist(strsplit(opt$factInput, "|", fixed = TRUE))
297 factordata <- list() 298 factordata <- list()
298 for (fact in factors) { 299 for (fact in factors) {
299 newfact <- unlist(strsplit(fact, split = "::")) 300 newfact <- unlist(strsplit(fact, split = "::"))
300 factordata <- rbind(factordata, newfact) 301 factordata <- rbind(factordata, newfact)
301 } # Factors have the form: FACT_NAME::LEVEL,LEVEL,LEVEL,LEVEL,... The first factor is the Primary Factor. 302 } # Factors have the form: FACT_NAME::LEVEL,LEVEL,LEVEL,LEVEL,... The first factor is the Primary Factor.
302 303
303 # Set the row names to be the name of the factor and delete first row 304 # Set the row names to be the name of the factor and delete first row
304 row.names(factordata) <- factordata[, 1] 305 row.names(factordata) <- factordata[, 1]
305 factordata <- factordata[, -1] 306 factordata <- factordata[, -1]
306 factordata <- sapply(factordata, sanitise_groups) 307 factordata <- sapply(factordata, sanitise_groups)
307 factordata <- sapply(factordata, strsplit, split = ",") 308 factordata <- sapply(factordata, strsplit, split = ",")
308 factordata <- sapply(factordata, make.names) 309 factordata <- sapply(factordata, make.names)
309 # Transform factor data into data frame of R factor objects 310 # Transform factor data into data frame of R factor objects
310 factors <- data.frame(factordata) 311 factors <- data.frame(factordata)
311 } 312 }
312 } 313 }
313 314
314 # if annotation file provided 315 # if annotation file provided
315 if (have_anno) { 316 if (have_anno) {
316 geneanno <- read.table(opt$annoPath, header = TRUE, sep = "\t", quote = "", strip.white = TRUE, stringsAsFactors = FALSE) 317 geneanno <- read.table(opt$annoPath, header = TRUE, sep = "\t", quote = "", strip.white = TRUE, stringsAsFactors = FALSE)
317 } 318 }
318 319
319 # Create output directory 320 # Create output directory
320 out_path <- opt$outPath 321 out_path <- opt$outPath
321 dir.create(out_path, showWarnings = FALSE) 322 dir.create(out_path, showWarnings = FALSE)
322 323
323 # Check if contrastData is a file or not 324 # Check if contrastData is a file or not
324 if (file.exists(opt$contrastData)) { 325 if (file.exists(opt$contrastData)) {
325 contrast_data <- unlist(read.table(opt$contrastData, sep = "\t", header = TRUE)[[1]]) 326 contrast_data <- unlist(read.table(opt$contrastData, sep = "\t", header = TRUE)[[1]])
326 } else { 327 } else {
327 # Split up contrasts separated by comma into a vector then sanitise 328 # Split up contrasts separated by comma into a vector then sanitise
328 contrast_data <- unlist(strsplit(opt$contrastData, split = ",")) 329 contrast_data <- unlist(strsplit(opt$contrastData, split = ","))
329 } 330 }
330 contrast_data <- sanitise_equation(contrast_data) 331 contrast_data <- sanitise_equation(contrast_data)
331 contrast_data <- gsub(" ", ".", contrast_data, fixed = TRUE) 332 contrast_data <- gsub(" ", ".", contrast_data, fixed = TRUE)
332 333
333 bcv_pdf <- make_out("bcvplot.pdf") 334 bcv_pdf <- make_out("bcvplot.pdf")
335 ql_pdf <- make_out("qlplot.pdf") 336 ql_pdf <- make_out("qlplot.pdf")
336 ql_png <- make_out("qlplot.png") 337 ql_png <- make_out("qlplot.png")
337 mds_pdf <- character() # Initialise character vector 338 mds_pdf <- character() # Initialise character vector
338 mds_png <- character() 339 mds_png <- character()
339 for (i in seq_len(ncol(factors))) { 340 for (i in seq_len(ncol(factors))) {
340 mds_pdf[i] <- make_out(paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".pdf")) 341 mds_pdf[i] <- make_out(paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".pdf"))
341 mds_png[i] <- make_out(paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".png")) 342 mds_png[i] <- make_out(paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".png"))
342 } 343 }
343 md_pdf <- character() 344 md_pdf <- character()
344 md_png <- character() 345 md_png <- character()
345 top_out <- character() 346 top_out <- character()
346 for (i in seq_along(contrast_data)) { 347 for (i in seq_along(contrast_data)) {
347 md_pdf[i] <- make_out(paste0("mdplot_", sanitise_basename(contrast_data[i]), ".pdf")) 348 md_pdf[i] <- make_out(paste0("mdplot_", sanitise_basename(contrast_data[i]), ".pdf"))
348 md_png[i] <- make_out(paste0("mdplot_", sanitise_basename(contrast_data[i]), ".png")) 349 md_png[i] <- make_out(paste0("mdplot_", sanitise_basename(contrast_data[i]), ".png"))
349 top_out[i] <- make_out(paste0("edgeR_", sanitise_basename(contrast_data[i]), ".tsv")) 350 top_out[i] <- make_out(paste0("edgeR_", sanitise_basename(contrast_data[i]), ".tsv"))
350 } # Save output paths for each contrast as vectors 351 } # Save output paths for each contrast as vectors
351 norm_out <- make_out("edgeR_normcounts.tsv") 352 norm_out <- make_out("edgeR_normcounts.tsv")
352 rda_out <- make_out("edgeR_analysis.RData") 353 rda_out <- make_out("edgeR_analysis.RData")
353 session_out <- make_out("session_info.txt") 354 session_out <- make_out("session_info.txt")
354 355
368 369
369 # Extract counts and annotation data 370 # Extract counts and annotation data
370 data <- list() 371 data <- list()
371 data$counts <- counts 372 data$counts <- counts
372 if (have_anno) { 373 if (have_anno) {
373 # order annotation by genes in counts (assumes gene ids are in 1st column of geneanno) 374 # order annotation by genes in counts (assumes gene ids are in 1st column of geneanno)
374 annoord <- geneanno[match(row.names(counts), geneanno[, 1]), ] 375 annoord <- geneanno[match(row.names(counts), geneanno[, 1]), ]
375 data$genes <- annoord 376 data$genes <- annoord
376 } else { 377 } else {
377 data$genes <- data.frame(GeneID = row.names(counts)) 378 data$genes <- data.frame(GeneID = row.names(counts))
378 } 379 }
379 380
380 # If filter crieteria set, filter out genes that do not have a required cpm/counts in a required number of 381 # If filter crieteria set, filter out genes that do not have a required cpm/counts in a required number of
381 # samples. Default is no filtering 382 # samples. Default is no filtering
382 prefilter_count <- nrow(data$counts) 383 prefilter_count <- nrow(data$counts)
383 384
384 if (filt_cpm || filt_smpcount || filt_totcount) { 385 if (filt_cpm || filt_smpcount || filt_totcount) {
385 if (filt_totcount) { 386 if (filt_totcount) {
386 keep <- rowSums(data$counts) >= opt$cntReq 387 keep <- rowSums(data$counts) >= opt$cntReq
387 } else if (filt_smpcount) { 388 } else if (filt_smpcount) {
388 keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq 389 keep <- rowSums(data$counts >= opt$cntReq) >= opt$sampleReq
389 } else if (filt_cpm) { 390 } else if (filt_cpm) {
390 keep <- rowSums(cpm(data$counts) >= opt$cpmReq) >= opt$sampleReq 391 keep <- rowSums(cpm(data$counts) >= opt$cpmReq) >= opt$sampleReq
391 } 392 }
392 393
393 data$counts <- data$counts[keep, ] 394 data$counts <- data$counts[keep, ]
394 data$genes <- data$genes[keep, , drop = FALSE] 395 data$genes <- data$genes[keep, , drop = FALSE]
395 } 396 }
396 397
397 postfilter_count <- nrow(data$counts) 398 postfilter_count <- nrow(data$counts)
398 filtered_count <- prefilter_count - postfilter_count 399 filtered_count <- prefilter_count - postfilter_count
399 400
409 data$samples <- factors 410 data$samples <- factors
410 data$genes <- genes 411 data$genes <- genes
411 412
412 413
413 if (!is.null(opt$formula)) { 414 if (!is.null(opt$formula)) {
414 formula <- opt$formula 415 formula <- opt$formula
415 # sanitisation can be getting rid of the "~" 416 # sanitisation can be getting rid of the "~"
416 if (!startsWith(formula, "~")) { 417 if (!startsWith(formula, "~")) {
417 formula <- paste0("~", formula) 418 formula <- paste0("~", formula)
418 } 419 }
419 } else { 420 } else {
420 formula <- "~0" 421 formula <- "~0"
421 for (i in seq_along(factor_list)) { 422 for (i in seq_along(factor_list)) {
422 formula <- paste(formula, factor_list[i], sep = "+") 423 formula <- paste(formula, factor_list[i], sep = "+")
423 } 424 }
424 } 425 }
425 426
426 formula <- formula(formula) 427 formula <- formula(formula)
427 design <- model.matrix(formula, factors) 428 design <- model.matrix(formula, factors)
428 429
429 for (i in seq_along(factor_list)) { 430 for (i in seq_along(factor_list)) {
430 colnames(design) <- gsub(factor_list[i], "", colnames(design), fixed = TRUE) 431 colnames(design) <- gsub(factor_list[i], "", colnames(design), fixed = TRUE)
431 } 432 }
432 433
433 # Calculating normalising factor, estimating dispersion 434 # Calculating normalising factor, estimating dispersion
434 data <- calcNormFactors(data, method = opt$normOpt) 435 data <- calcNormFactors(data, method = opt$normOpt)
435 436
436 if (want_robust) { 437 if (want_robust) {
437 data <- estimateDisp(data, design = design, robust = TRUE) 438 data <- estimateDisp(data, design = design, robust = TRUE)
438 } else { 439 } else {
439 data <- estimateDisp(data, design = design) 440 data <- estimateDisp(data, design = design)
440 } 441 }
441 442
442 # Generate contrasts information 443 # Generate contrasts information
443 contrasts <- makeContrasts(contrasts = contrast_data, levels = design) 444 contrasts <- makeContrasts(contrasts = contrast_data, levels = design)
444 445
464 link_data[1, ] <- c(link_name, link_addr) 465 link_data[1, ] <- c(link_name, link_addr)
465 invisible(dev.off()) 466 invisible(dev.off())
466 467
467 # If additional factors create additional MDS plots coloured by factor 468 # If additional factors create additional MDS plots coloured by factor
468 if (ncol(factors) > 1) { 469 if (ncol(factors) > 1) {
469 for (i in 2:ncol(factors)) { 470 for (i in 2:ncol(factors)) {
470 png(mds_png[i], width = 600, height = 600) 471 png(mds_png[i], width = 600, height = 600)
471 plotMDS(data, labels = labels, col = as.numeric(factors[, i]), cex = 0.8, main = paste("MDS Plot:", names(factors)[i])) 472 plotMDS(data, labels = labels, col = as.numeric(factors[, i]), cex = 0.8, main = paste("MDS Plot:", names(factors)[i]))
472 img_name <- paste0("MDS Plot_", sanitise_basename(names(factors)[i]), ".png") 473 img_name <- paste0("MDS Plot_", sanitise_basename(names(factors)[i]), ".png")
473 img_addr <- paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".png") 474 img_addr <- paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".png")
474 image_data <- rbind(image_data, c(img_name, img_addr)) 475 image_data <- rbind(image_data, c(img_name, img_addr))
475 invisible(dev.off()) 476 invisible(dev.off())
476 477
477 pdf(mds_pdf[i]) 478 pdf(mds_pdf[i])
478 plotMDS(data, labels = labels, col = as.numeric(factors[, i]), cex = 0.8, main = paste("MDS Plot:", names(factors)[i])) 479 plotMDS(data, labels = labels, col = as.numeric(factors[, i]), cex = 0.8, main = paste("MDS Plot:", names(factors)[i]))
479 link_name <- paste0("MDS Plot_", sanitise_basename(names(factors)[i]), ".pdf") 480 link_name <- paste0("MDS Plot_", sanitise_basename(names(factors)[i]), ".pdf")
480 link_addr <- paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".pdf") 481 link_addr <- paste0("mdsplot_", sanitise_basename(names(factors)[i]), ".pdf")
481 link_data <- rbind(link_data, c(link_name, link_addr)) 482 link_data <- rbind(link_data, c(link_name, link_addr))
482 invisible(dev.off()) 483 invisible(dev.off())
483 } 484 }
484 } 485 }
485 486
486 # BCV Plot 487 # BCV Plot
487 png(bcv_png, width = 600, height = 600) 488 png(bcv_png, width = 600, height = 600)
488 plotBCV(data, main = "BCV Plot") 489 plotBCV(data, main = "BCV Plot")
498 link_data <- rbind(link_data, c(link_name, link_addr)) 499 link_data <- rbind(link_data, c(link_name, link_addr))
499 invisible(dev.off()) 500 invisible(dev.off())
500 501
501 # Generate fit 502 # Generate fit
502 if (want_lrt) { 503 if (want_lrt) {
503 fit <- glmFit(data, design) 504 fit <- glmFit(data, design)
504 } else { 505 } else {
505 if (want_robust) { 506 if (want_robust) {
506 fit <- glmQLFit(data, design, robust = TRUE) 507 fit <- glmQLFit(data, design, robust = TRUE)
507 } else { 508 } else {
508 fit <- glmQLFit(data, design) 509 fit <- glmQLFit(data, design)
509 } 510 }
510 511
511 # Plot QL dispersions 512 # Plot QL dispersions
512 png(ql_png, width = 600, height = 600) 513 png(ql_png, width = 600, height = 600)
513 plotQLDisp(fit, main = "QL Plot") 514 plotQLDisp(fit, main = "QL Plot")
514 img_name <- "QL Plot" 515 img_name <- "QL Plot"
515 img_addr <- "qlplot.png" 516 img_addr <- "qlplot.png"
516 image_data <- rbind(image_data, c(img_name, img_addr)) 517 image_data <- rbind(image_data, c(img_name, img_addr))
517 invisible(dev.off()) 518 invisible(dev.off())
518 519
519 pdf(ql_pdf) 520 pdf(ql_pdf)
520 plotQLDisp(fit, main = "QL Plot") 521 plotQLDisp(fit, main = "QL Plot")
521 link_name <- "QL Plot.pdf" 522 link_name <- "QL Plot.pdf"
522 link_addr <- "qlplot.pdf" 523 link_addr <- "qlplot.pdf"
523 link_data <- rbind(link_data, c(link_name, link_addr)) 524 link_data <- rbind(link_data, c(link_name, link_addr))
524 invisible(dev.off()) 525 invisible(dev.off())
525 } 526 }
526 527
527 # Save normalised counts (log2cpm) 528 # Save normalised counts (log2cpm)
528 if (want_norm) { 529 if (want_norm) {
529 normalised_counts <- cpm(data, normalized.lib.sizes = TRUE, log = TRUE) 530 normalised_counts <- cpm(data, normalized.lib.sizes = TRUE, log = TRUE)
530 normalised_counts <- data.frame(data$genes, normalised_counts) 531 normalised_counts <- data.frame(data$genes, normalised_counts)
531 write.table(normalised_counts, file = norm_out, row.names = FALSE, sep = "\t", quote = FALSE) 532 write.table(normalised_counts, file = norm_out, row.names = FALSE, sep = "\t", quote = FALSE)
532 link_data <- rbind(link_data, c("edgeR_normcounts.tsv", "edgeR_normcounts.tsv")) 533 link_data <- rbind(link_data, c("edgeR_normcounts.tsv", "edgeR_normcounts.tsv"))
533 } 534 }
534 535
535 536
536 for (i in seq_along(contrast_data)) { 537 for (i in seq_along(contrast_data)) {
537 if (want_lrt) { 538 if (want_lrt) {
538 res <- glmLRT(fit, contrast = contrasts[, i]) 539 res <- glmLRT(fit, contrast = contrasts[, i])
539 } else { 540 } else {
540 res <- glmQLFTest(fit, contrast = contrasts[, i]) 541 res <- glmQLFTest(fit, contrast = contrasts[, i])
541 } 542 }
542 543
543 status <- decideTestsDGE(res, 544 status <- decideTestsDGE(res,
544 adjust.method = opt$pAdjOpt, p.value = opt$pValReq, 545 adjust.method = opt$pAdjOpt, p.value = opt$pValReq,
545 lfc = opt$lfcReq 546 lfc = opt$lfcReq
546 ) 547 )
547 sum_status <- summary(status) 548 sum_status <- summary(status)
548 549
549 # Collect counts for differential expression 550 # Collect counts for differential expression
550 up_count[i] <- sum_status["Up", ] 551 up_count[i] <- sum_status["Up", ]
551 down_count[i] <- sum_status["Down", ] 552 down_count[i] <- sum_status["Down", ]
552 flat_count[i] <- sum_status["NotSig", ] 553 flat_count[i] <- sum_status["NotSig", ]
553 554
554 # Write top expressions table 555 # Write top expressions table
555 top <- topTags(res, adjust.method = opt$pAdjOpt, n = Inf, sort.by = "PValue") 556 top <- topTags(res, adjust.method = opt$pAdjOpt, n = Inf, sort.by = "PValue")
556 write.table(top, file = top_out[i], row.names = FALSE, sep = "\t", quote = FALSE) 557 write.table(top, file = top_out[i], row.names = FALSE, sep = "\t", quote = FALSE)
557 558
558 link_name <- paste0("edgeR_", sanitise_basename(contrast_data[i]), ".tsv") 559 link_name <- paste0("edgeR_", sanitise_basename(contrast_data[i]), ".tsv")
559 link_addr <- paste0("edgeR_", sanitise_basename(contrast_data[i]), ".tsv") 560 link_addr <- paste0("edgeR_", sanitise_basename(contrast_data[i]), ".tsv")
560 link_data <- rbind(link_data, c(link_name, link_addr)) 561 link_data <- rbind(link_data, c(link_name, link_addr))
561 562
562 # Plot MD (log ratios vs mean difference) using limma package 563 # Plot MD (log ratios vs mean difference) using limma package
563 pdf(md_pdf[i]) 564 pdf(md_pdf[i])
564 limma::plotMD(res, 565 limma::plotMD(res,
565 status = status, 566 status = status,
566 main = paste("MD Plot:", unmake_names(contrast_data[i])), 567 main = paste("MD Plot:", unmake_names(contrast_data[i])),
567 hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1), 568 hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1),
568 xlab = "Average Expression", ylab = "logFC" 569 xlab = "Average Expression", ylab = "logFC"
569 ) 570 )
570 571
571 abline(h = 0, col = "grey", lty = 2) 572 abline(h = 0, col = "grey", lty = 2)
572 573
573 link_name <- paste0("MD Plot_", sanitise_basename(contrast_data[i]), ".pdf") 574 link_name <- paste0("MD Plot_", sanitise_basename(contrast_data[i]), ".pdf")
574 link_addr <- paste0("mdplot_", sanitise_basename(contrast_data[i]), ".pdf") 575 link_addr <- paste0("mdplot_", sanitise_basename(contrast_data[i]), ".pdf")
575 link_data <- rbind(link_data, c(link_name, link_addr)) 576 link_data <- rbind(link_data, c(link_name, link_addr))
576 invisible(dev.off()) 577 invisible(dev.off())
577 578
578 png(md_png[i], height = 600, width = 600) 579 png(md_png[i], height = 600, width = 600)
579 limma::plotMD(res, 580 limma::plotMD(res,
580 status = status, 581 status = status,
581 main = paste("MD Plot:", unmake_names(contrast_data[i])), 582 main = paste("MD Plot:", unmake_names(contrast_data[i])),
582 hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1), 583 hl.col = alpha(c("firebrick", "blue"), 0.4), values = c(1, -1),
583 xlab = "Average Expression", ylab = "logFC" 584 xlab = "Average Expression", ylab = "logFC"
584 ) 585 )
585 586
586 abline(h = 0, col = "grey", lty = 2) 587 abline(h = 0, col = "grey", lty = 2)
587 588
588 img_name <- paste0("MD Plot_", sanitise_basename(contrast_data[i]), ".png") 589 img_name <- paste0("MD Plot_", sanitise_basename(contrast_data[i]), ".png")
589 img_addr <- paste0("mdplot_", sanitise_basename(contrast_data[i]), ".png") 590 img_addr <- paste0("mdplot_", sanitise_basename(contrast_data[i]), ".png")
590 image_data <- rbind(image_data, c(img_name, img_addr)) 591 image_data <- rbind(image_data, c(img_name, img_addr))
591 invisible(dev.off()) 592 invisible(dev.off())
592 } 593 }
593 sig_diff <- data.frame(Up = up_count, Flat = flat_count, Down = down_count) 594 sig_diff <- data.frame(Up = up_count, Flat = flat_count, Down = down_count)
594 row.names(sig_diff) <- contrast_data 595 row.names(sig_diff) <- contrast_data
595 596
596 # Save relevant items as rda object 597 # Save relevant items as rda object
597 if (want_rda) { 598 if (want_rda) {
598 if (want_norm) { 599 if (want_norm) {
599 save(counts, data, status, normalised_counts, labels, factors, fit, res, top, contrasts, design, 600 save(counts, data, status, normalised_counts, labels, factors, fit, res, top, contrasts, design,
600 file = rda_out, ascii = TRUE 601 file = rda_out, ascii = TRUE
601 ) 602 )
602 } else { 603 } else {
603 save(counts, data, status, labels, factors, fit, res, top, contrasts, design, 604 save(counts, data, status, labels, factors, fit, res, top, contrasts, design,
604 file = rda_out, ascii = TRUE 605 file = rda_out, ascii = TRUE
605 ) 606 )
606 } 607 }
607 link_data <- rbind(link_data, c("edgeR_analysis.RData", "edgeR_analysis.RData")) 608 link_data <- rbind(link_data, c("edgeR_analysis.RData", "edgeR_analysis.RData"))
608 } 609 }
609 610
610 # Record session info 611 # Record session info
611 writeLines(capture.output(sessionInfo()), session_out) 612 writeLines(capture.output(sessionInfo()), session_out)
612 link_data <- rbind(link_data, c("Session Info", "session_info.txt")) 613 link_data <- rbind(link_data, c("Session Info", "session_info.txt"))
630 cata("Links to PDF copies of plots are in 'Plots' section below.<br />\n") 631 cata("Links to PDF copies of plots are in 'Plots' section below.<br />\n")
631 632
632 html_image(image_data$Link[1], image_data$Label[1]) 633 html_image(image_data$Link[1], image_data$Label[1])
633 634
634 for (i in 2:nrow(image_data)) { 635 for (i in 2:nrow(image_data)) {
635 html_image(image_data$Link[i], image_data$Label[i]) 636 html_image(image_data$Link[i], image_data$Label[i])
636 } 637 }
637 638
638 cata("<h4>Differential Expression Counts:</h4>\n") 639 cata("<h4>Differential Expression Counts:</h4>\n")
639 640
640 cata("<table border=\"1\" cellpadding=\"4\">\n") 641 cata("<table border=\"1\" cellpadding=\"4\">\n")
641 cata("<tr>\n") 642 cata("<tr>\n")
642 table_item() 643 table_item()
643 for (i in colnames(sig_diff)) { 644 for (i in colnames(sig_diff)) {
644 table_head_item(i) 645 table_head_item(i)
645 } 646 }
646 cata("</tr>\n") 647 cata("</tr>\n")
647 for (i in seq_len(nrow(sig_diff))) { 648 for (i in seq_len(nrow(sig_diff))) {
648 cata("<tr>\n") 649 cata("<tr>\n")
649 table_head_item(unmake_names(row.names(sig_diff)[i])) 650 table_head_item(unmake_names(row.names(sig_diff)[i]))
650 for (j in seq_len(ncol(sig_diff))) { 651 for (j in seq_len(ncol(sig_diff))) {
651 table_item(as.character(sig_diff[i, j])) 652 table_item(as.character(sig_diff[i, j]))
652 } 653 }
653 cata("</tr>\n") 654 cata("</tr>\n")
654 } 655 }
655 cata("</table>") 656 cata("</table>")
656 657
657 cata("<h4>Plots:</h4>\n") 658 cata("<h4>Plots:</h4>\n")
658 for (i in seq_len(nrow(link_data))) { 659 for (i in seq_len(nrow(link_data))) {
659 if (grepl(".pdf", link_data$Link[i])) { 660 if (grepl(".pdf", link_data$Link[i])) {
660 html_link(link_data$Link[i], link_data$Label[i]) 661 html_link(link_data$Link[i], link_data$Label[i])
661 } 662 }
662 } 663 }
663 664
664 cata("<h4>Tables:</h4>\n") 665 cata("<h4>Tables:</h4>\n")
665 for (i in seq_len(nrow(link_data))) { 666 for (i in seq_len(nrow(link_data))) {
666 if (grepl(".tsv", link_data$Link[i])) { 667 if (grepl(".tsv", link_data$Link[i])) {
667 html_link(link_data$Link[i], link_data$Label[i]) 668 html_link(link_data$Link[i], link_data$Label[i])
668 } 669 }
669 } 670 }
670 671
671 if (want_rda) { 672 if (want_rda) {
672 cata("<h4>R Data Objects:</h4>\n") 673 cata("<h4>R Data Objects:</h4>\n")
673 for (i in seq_len(nrow(link_data))) { 674 for (i in seq_len(nrow(link_data))) {
674 if (grepl(".RData", link_data$Link[i])) { 675 if (grepl(".RData", link_data$Link[i])) {
675 html_link(link_data$Link[i], link_data$Label[i]) 676 html_link(link_data$Link[i], link_data$Label[i])
676 } 677 }
677 } 678 }
678 } 679 }
679 680
680 cata("<p>Alt-click links to download file.</p>\n") 681 cata("<p>Alt-click links to download file.</p>\n")
681 cata("<p>Click floppy disc icon associated history item to download ") 682 cata("<p>Click floppy disc icon associated history item to download ")
682 cata("all files.</p>\n") 683 cata("all files.</p>\n")
684 685
685 cata("<h4>Additional Information</h4>\n") 686 cata("<h4>Additional Information</h4>\n")
686 cata("<ul>\n") 687 cata("<ul>\n")
687 688
688 if (filt_cpm || filt_smpcount || filt_totcount) { 689 if (filt_cpm || filt_smpcount || filt_totcount) {
689 if (filt_cpm) { 690 if (filt_cpm) {
690 temp_str <- paste( 691 temp_str <- paste(
691 "Genes without more than", opt$cpmReq, 692 "Genes without more than", opt$cpmReq,
692 "CPM in at least", opt$sampleReq, "samples are insignificant", 693 "CPM in at least", opt$sampleReq, "samples are insignificant",
693 "and filtered out." 694 "and filtered out."
695 )
696 } else if (filt_smpcount) {
697 temp_str <- paste(
698 "Genes without more than", opt$cntReq,
699 "counts in at least", opt$sampleReq, "samples are insignificant",
700 "and filtered out."
701 )
702 } else if (filt_totcount) {
703 temp_str <- paste(
704 "Genes without more than", opt$cntReq,
705 "counts, after summing counts for all samples, are insignificant",
706 "and filtered out."
707 )
708 }
709
710 list_item(temp_str)
711 filter_prop <- round(filtered_count / prefilter_count * 100, digits = 2)
712 temp_str <- paste0(
713 filtered_count, " of ", prefilter_count, " (", filter_prop,
714 "%) genes were filtered out for low expression."
694 ) 715 )
695 } else if (filt_smpcount) { 716 list_item(temp_str)
696 temp_str <- paste(
697 "Genes without more than", opt$cntReq,
698 "counts in at least", opt$sampleReq, "samples are insignificant",
699 "and filtered out."
700 )
701 } else if (filt_totcount) {
702 temp_str <- paste(
703 "Genes without more than", opt$cntReq,
704 "counts, after summing counts for all samples, are insignificant",
705 "and filtered out."
706 )
707 }
708
709 list_item(temp_str)
710 filter_prop <- round(filtered_count / prefilter_count * 100, digits = 2)
711 temp_str <- paste0(
712 filtered_count, " of ", prefilter_count, " (", filter_prop,
713 "%) genes were filtered out for low expression."
714 )
715 list_item(temp_str)
716 } 717 }
717 list_item(opt$normOpt, " was the method used to normalise library sizes.") 718 list_item(opt$normOpt, " was the method used to normalise library sizes.")
718 if (want_lrt) { 719 if (want_lrt) {
719 list_item("The edgeR likelihood ratio test was used.") 720 list_item("The edgeR likelihood ratio test was used.")
720 } else { 721 } else {
721 if (want_robust) { 722 if (want_robust) {
722 list_item("The edgeR quasi-likelihood test was used with robust settings (robust=TRUE with estimateDisp and glmQLFit).") 723 list_item("The edgeR quasi-likelihood test was used with robust settings (robust=TRUE with estimateDisp and glmQLFit).")
723 } else { 724 } else {
724 list_item("The edgeR quasi-likelihood test was used.") 725 list_item("The edgeR quasi-likelihood test was used.")
725 } 726 }
726 } 727 }
727 if (opt$pAdjOpt != "none") { 728 if (opt$pAdjOpt != "none") {
728 if (opt$pAdjOpt == "BH" || opt$pAdjOpt == "BY") { 729 if (opt$pAdjOpt == "BH" || opt$pAdjOpt == "BY") {
730 temp_str <- paste0(
731 "MD-Plot highlighted genes are significant at FDR ",
732 "of ", opt$pValReq, " and exhibit log2-fold-change of at ",
733 "least ", opt$lfcReq, "."
734 )
735 list_item(temp_str)
736 } else if (opt$pAdjOpt == "holm") {
737 temp_str <- paste0(
738 "MD-Plot highlighted genes are significant at adjusted ",
739 "p-value of ", opt$pValReq, " by the Holm(1979) ",
740 "method, and exhibit log2-fold-change of at least ",
741 opt$lfcReq, "."
742 )
743 list_item(temp_str)
744 }
745 } else {
729 temp_str <- paste0( 746 temp_str <- paste0(
730 "MD-Plot highlighted genes are significant at FDR ", 747 "MD-Plot highlighted genes are significant at p-value ",
731 "of ", opt$pValReq, " and exhibit log2-fold-change of at ", 748 "of ", opt$pValReq, " and exhibit log2-fold-change of at ",
732 "least ", opt$lfcReq, "." 749 "least ", opt$lfcReq, "."
733 ) 750 )
734 list_item(temp_str) 751 list_item(temp_str)
735 } else if (opt$pAdjOpt == "holm") {
736 temp_str <- paste0(
737 "MD-Plot highlighted genes are significant at adjusted ",
738 "p-value of ", opt$pValReq, " by the Holm(1979) ",
739 "method, and exhibit log2-fold-change of at least ",
740 opt$lfcReq, "."
741 )
742 list_item(temp_str)
743 }
744 } else {
745 temp_str <- paste0(
746 "MD-Plot highlighted genes are significant at p-value ",
747 "of ", opt$pValReq, " and exhibit log2-fold-change of at ",
748 "least ", opt$lfcReq, "."
749 )
750 list_item(temp_str)
751 } 752 }
752 cata("</ul>\n") 753 cata("</ul>\n")
753 754
754 cata("<h4>Summary of experimental data:</h4>\n") 755 cata("<h4>Summary of experimental data:</h4>\n")
755 756
759 cata("<tr>\n") 760 cata("<tr>\n")
760 table_head_item("SampleID") 761 table_head_item("SampleID")
761 table_head_item(names(factors)[1], " (Primary Factor)") 762 table_head_item(names(factors)[1], " (Primary Factor)")
762 763
763 if (ncol(factors) > 1) { 764 if (ncol(factors) > 1) {
764 for (i in names(factors)[2:length(names(factors))]) { 765 for (i in names(factors)[2:length(names(factors))]) {
765 table_head_item(i) 766 table_head_item(i)
766 } 767 }
767 cata("</tr>\n") 768 cata("</tr>\n")
768 } 769 }
769 770
770 for (i in seq_len(nrow((factors)))) { 771 for (i in seq_len(nrow((factors)))) {
771 cata("<tr>\n") 772 cata("<tr>\n")
772 table_head_item(row.names(factors)[i]) 773 table_head_item(row.names(factors)[i])
773 for (j in seq_len(ncol(factors))) { 774 for (j in seq_len(ncol(factors))) {
774 table_item(as.character(unmake_names(factors[i, j]))) 775 table_item(as.character(unmake_names(factors[i, j])))
775 } 776 }
776 cata("</tr>\n") 777 cata("</tr>\n")
777 } 778 }
778 cata("</table>") 779 cata("</table>")
779 780
780 for (i in seq_len(nrow(link_data))) { 781 for (i in seq_len(nrow(link_data))) {
781 if (grepl("session_info", link_data$Link[i])) { 782 if (grepl("session_info", link_data$Link[i])) {
782 html_link(link_data$Link[i], link_data$Label[i]) 783 html_link(link_data$Link[i], link_data$Label[i])
783 } 784 }
784 } 785 }
785 786
786 cata("<table border=\"0\">\n") 787 cata("<table border=\"0\">\n")
787 cata("<tr>\n") 788 cata("<tr>\n")
788 table_item("Task started at:") 789 table_item("Task started at:")