Mercurial > repos > shians > shrnaseq
annotate hairpinTool.R @ 4:f8af57d6f60b
Fixed bug causing fastq input to break
- Moved misplaced error check causing fastq inputs to not function
| author | shian_su <registertonysu@gmail.com> |
|---|---|
| date | Mon, 24 Feb 2014 14:41:39 +1100 |
| parents | 076ca575208f |
| children | 3d04308a99f9 |
| rev | line source |
|---|---|
| 2 | 1 # Record starting time |
| 2 timeStart <- as.character(Sys.time()) | |
| 3 | |
| 4 # Loading and checking required packages | |
| 5 library(methods, quietly=TRUE, warn.conflicts=FALSE) | |
| 6 library(statmod, quietly=TRUE, warn.conflicts=FALSE) | |
| 7 library(splines, quietly=TRUE, warn.conflicts=FALSE) | |
| 8 library(edgeR, quietly=TRUE, warn.conflicts=FALSE) | |
| 9 library(limma, quietly=TRUE, warn.conflicts=FALSE) | |
| 10 | |
| 11 if (packageVersion("edgeR") < "3.5.23") { | |
| 12 message("Please update 'edgeR' to version >= 3.5.23 to run this script") | |
| 13 } | |
| 14 | |
| 15 ################################################################################ | |
| 16 ### Function declarations | |
| 17 ################################################################################ | |
| 18 | |
| 19 # Function to sanitise contrast equations so there are no whitespaces | |
| 20 # surrounding the arithmetic operators, leading or trailing whitespace | |
| 21 sanitiseEquation <- function(equation) { | |
| 22 equation <- gsub(" *[+] *", "+", equation) | |
| 23 equation <- gsub(" *[-] *", "-", equation) | |
| 24 equation <- gsub(" *[/] *", "/", equation) | |
| 25 equation <- gsub(" *[*] *", "*", equation) | |
| 26 equation <- gsub("^\\s+|\\s+$", "", equation) | |
| 27 return(equation) | |
| 28 } | |
| 29 | |
| 30 # Function to sanitise group information | |
| 31 sanitiseGroups <- function(string) { | |
| 32 string <- gsub(" *[,] *", ",", string) | |
| 33 string <- gsub("^\\s+|\\s+$", "", string) | |
| 34 return(string) | |
| 35 } | |
| 36 | |
| 37 # Function to change periods to whitespace in a string | |
| 38 unmake.names <- function(string) { | |
| 39 string <- gsub(".", " ", string, fixed=TRUE) | |
| 40 return(string) | |
| 41 } | |
| 42 | |
| 43 # Function has string input and generates an output path string | |
| 44 makeOut <- function(filename) { | |
| 45 return(paste0(folderPath, "/", filename)) | |
| 46 } | |
| 47 | |
| 48 # Function has string input and generates both a pdf and png output strings | |
| 49 imgOut <- function(filename) { | |
| 50 assign(paste0(filename, "Png"), makeOut(paste0(filename,".png")), | |
| 51 envir = .GlobalEnv) | |
| 52 assign(paste0(filename, "Pdf"), makeOut(paste0(filename,".pdf")), | |
| 53 envir = .GlobalEnv) | |
| 54 } | |
| 55 | |
| 56 # Create cat function default path set, default seperator empty and appending | |
| 57 # true by default (Ripped straight from the cat function with altered argument | |
| 58 # defaults) | |
| 59 cata <- function(..., file = htmlPath, sep = "", fill = FALSE, labels = NULL, | |
| 60 append = TRUE) { | |
| 61 if (is.character(file)) | |
| 62 if (file == "") | |
| 63 file <- stdout() | |
| 64 else if (substring(file, 1L, 1L) == "|") { | |
| 65 file <- pipe(substring(file, 2L), "w") | |
| 66 on.exit(close(file)) | |
| 67 } | |
| 68 else { | |
| 69 file <- file(file, ifelse(append, "a", "w")) | |
| 70 on.exit(close(file)) | |
| 71 } | |
| 72 .Internal(cat(list(...), file, sep, fill, labels, append)) | |
| 73 } | |
| 74 | |
| 75 # Function to write code for html head and title | |
| 76 HtmlHead <- function(title) { | |
| 77 cata("<head>\n") | |
| 78 cata("<title>", title, "</title>\n") | |
| 79 cata("</head>\n") | |
| 80 } | |
| 81 | |
| 82 # Function to write code for html links | |
| 83 HtmlLink <- function(address, label=address) { | |
| 84 cata("<a href=\"", address, "\" target=\"_blank\">", label, "</a><br />\n") | |
| 85 } | |
| 86 | |
| 87 # Function to write code for html images | |
| 88 HtmlImage <- function(source, label=source, height=600, width=600) { | |
| 89 cata("<img src=\"", source, "\" alt=\"", label, "\" height=\"", height) | |
| 90 cata("\" width=\"", width, "\"/>\n") | |
| 91 } | |
| 92 | |
| 93 # Function to write code for html list items | |
| 94 ListItem <- function(...) { | |
| 95 cata("<li>", ..., "</li>\n") | |
| 96 } | |
| 97 | |
| 98 TableItem <- function(...) { | |
| 99 cata("<td>", ..., "</td>\n") | |
| 100 } | |
| 101 | |
| 102 TableHeadItem <- function(...) { | |
| 103 cata("<th>", ..., "</th>\n") | |
| 104 } | |
| 105 ################################################################################ | |
| 106 ### Input Processing | |
| 107 ################################################################################ | |
| 108 | |
| 109 # Grabbing arguments from command line | |
| 110 argv <- commandArgs(TRUE) | |
| 111 | |
| 112 # Remove fastq file paths after collecting from argument vector | |
| 113 inputType <- as.character(argv[1]) | |
| 114 if (inputType=="fastq") { | |
| 115 fastqPath <- as.character(gsub("fastq::", "", argv[grepl("fastq::", argv)], | |
| 116 fixed=TRUE)) | |
| 117 argv <- argv[!grepl("fastq::", argv, fixed=TRUE)] | |
| 118 hairpinPath <- as.character(argv[2]) | |
| 119 samplePath <- as.character(argv[3]) | |
| 120 barStart <- as.numeric(argv[4]) | |
| 121 barEnd <- as.numeric(argv[5]) | |
| 122 hpStart <- as.numeric(argv[6]) | |
| 123 hpEnd <- as.numeric(argv[7]) | |
| 124 } else if (inputType=="counts") { | |
| 125 countPath <- as.character(argv[2]) | |
| 126 annoPath <- as.character(argv[3]) | |
| 127 samplePath <- as.character(argv[4]) | |
| 128 } | |
| 129 | |
| 130 cpmReq <- as.numeric(argv[8]) | |
| 131 sampleReq <- as.numeric(argv[9]) | |
| 132 fdrThresh <- as.numeric(argv[10]) | |
| 133 lfcThresh <- as.numeric(argv[11]) | |
| 134 workMode <- as.character(argv[12]) | |
| 135 htmlPath <- as.character(argv[13]) | |
| 136 folderPath <- as.character(argv[14]) | |
| 137 if (workMode=="classic") { | |
| 138 pairData <- character() | |
| 139 pairData[2] <- as.character(argv[15]) | |
| 140 pairData[1] <- as.character(argv[16]) | |
| 141 } else if (workMode=="glm") { | |
| 142 contrastData <- as.character(argv[15]) | |
| 143 roastOpt <- as.character(argv[16]) | |
| 144 hairpinReq <- as.numeric(argv[17]) | |
| 145 selectOpt <- as.character(argv[18]) | |
| 146 selectVals <- as.character(argv[19]) | |
| 147 } | |
| 148 | |
| 149 # Read in inputs | |
| 150 if (inputType=="fastq") { | |
| 151 samples <- read.table(samplePath, header=TRUE, sep="\t") | |
| 152 hairpins <- read.table(hairpinPath, header=TRUE, sep="\t") | |
| 153 } else if (inputType=="counts") { | |
| 154 samples <- read.table(samplePath, header=TRUE, sep="\t") | |
| 155 counts <- read.table(countPath, header=TRUE, sep="\t") | |
| 156 anno <- read.table(annoPath, header=TRUE, sep="\t") | |
| 157 } | |
| 158 ###################### Check inputs for correctness ############################ | |
| 159 samples$ID <- make.names(samples$ID) | |
| 160 | |
| 161 if (!any(grepl("group", names(samples)))) { | |
| 162 stop("'group' column not specified in sample annotation file") | |
| 163 } # Check if grouping variable has been specified | |
| 164 | |
| 165 if (any(table(samples$ID)>1)){ | |
| 166 tab <- table(samples$ID) | |
| 167 offenders <- paste(names(tab[tab>1]), collapse=", ") | |
| 168 offenders <- unmake.names(offenders) | |
| 169 stop("ID column of sample annotation must have unique values, values ", | |
| 170 offenders, " are repeated") | |
| 171 } # Check that IDs in sample annotation are unique | |
| 172 | |
| 173 if (inputType=="fastq") { | |
| 174 | |
| 175 if (any(table(hairpins$ID)>1)){ | |
| 176 tab <- table(hairpins$ID) | |
| 177 offenders <- paste(names(tab[tab>1]), collapse=", ") | |
| 178 stop("ID column of hairpin annotation must have unique values, values ", | |
| 179 offenders, " are repeated") | |
| 180 } # Check that IDs in hairpin annotation are unique | |
|
4
f8af57d6f60b
Fixed bug causing fastq input to break
shian_su <registertonysu@gmail.com>
parents:
2
diff
changeset
|
181 |
| 2 | 182 } else if (inputType=="counts") { |
|
4
f8af57d6f60b
Fixed bug causing fastq input to break
shian_su <registertonysu@gmail.com>
parents:
2
diff
changeset
|
183 if (any(is.na(match(samples$ID, colnames(counts))))) { |
|
f8af57d6f60b
Fixed bug causing fastq input to break
shian_su <registertonysu@gmail.com>
parents:
2
diff
changeset
|
184 stop("not all samples have groups specified") |
|
f8af57d6f60b
Fixed bug causing fastq input to break
shian_su <registertonysu@gmail.com>
parents:
2
diff
changeset
|
185 } # Check that a group has be specifed for each sample |
| 2 | 186 |
| 187 if (any(table(counts$ID)>1)){ | |
| 188 tab <- table(counts$ID) | |
| 189 offenders <- paste(names(tab[tab>1]), collapse=", ") | |
| 190 stop("ID column of count table must have unique values, values ", | |
| 191 offenders, " are repeated") | |
| 192 } # Check that IDs in count table are unique | |
| 193 } | |
| 194 ################################################################################ | |
| 195 | |
| 196 # Process arguments | |
| 197 if (workMode=="glm") { | |
| 198 if (roastOpt=="yes") { | |
| 199 wantRoast <- TRUE | |
| 200 } else { | |
| 201 wantRoast <- FALSE | |
| 202 } | |
| 203 } | |
| 204 | |
| 205 # Split up contrasts seperated by comma into a vector and replace spaces with | |
| 206 # periods | |
| 207 if (exists("contrastData")) { | |
| 208 contrastData <- unlist(strsplit(contrastData, split=",")) | |
| 209 contrastData <- sanitiseEquation(contrastData) | |
| 210 contrastData <- gsub(" ", ".", contrastData, fixed=TRUE) | |
| 211 } | |
| 212 | |
| 213 # Replace spaces with periods in pair data | |
| 214 if (exists("pairData")) { | |
| 215 pairData <- make.names(pairData) | |
| 216 } | |
| 217 | |
| 218 # Generate output folder and paths | |
| 219 dir.create(folderPath) | |
| 220 | |
| 221 # Generate links for outputs | |
| 222 imgOut("barHairpin") | |
| 223 imgOut("barIndex") | |
| 224 imgOut("mds") | |
| 225 imgOut("bcv") | |
| 226 if (workMode == "classic") { | |
| 227 smearPng <- makeOut(paste0("smear(", pairData[2], "-", pairData[1],").png")) | |
| 228 smearPdf <- makeOut(paste0("smear(", pairData[2], "-", pairData[1],").pdf")) | |
| 229 topOut <- makeOut(paste0("toptag(", pairData[2], "-", pairData[1],").tsv")) | |
| 230 } else if (workMode=="glm") { | |
| 231 smearPng <- character() | |
| 232 smearPdf <- character() | |
| 233 topOut <- character() | |
| 234 roastOut <- character() | |
| 235 barcodePng <- character() | |
| 236 barcodePdf <- character() | |
| 237 for (i in 1:length(contrastData)) { | |
| 238 smearPng[i] <- makeOut(paste0("smear(", contrastData[i], ").png")) | |
| 239 smearPdf[i] <- makeOut(paste0("smear(", contrastData[i], ").pdf")) | |
| 240 topOut[i] <- makeOut(paste0("toptag(", contrastData[i], ").tsv")) | |
| 241 roastOut[i] <- makeOut(paste0("roast(", contrastData[i], ").tsv")) | |
| 242 barcodePng[i] <- makeOut(paste0("barcode(", contrastData[i], ").png")) | |
| 243 barcodePdf[i] <- makeOut(paste0("barcode(", contrastData[i], ").pdf")) | |
| 244 } | |
| 245 } | |
| 246 # Initialise data for html links and images, table with the link label and | |
| 247 # link address | |
| 248 linkData <- data.frame(Label=character(), Link=character(), | |
| 249 stringsAsFactors=FALSE) | |
| 250 imageData <- data.frame(Label=character(), Link=character(), | |
| 251 stringsAsFactors=FALSE) | |
| 252 ################################################################################ | |
| 253 ### Data Processing | |
| 254 ################################################################################ | |
| 255 | |
| 256 # Transform gene selection from string into index values for mroast | |
| 257 if (workMode=="glm") { | |
| 258 if (selectOpt=="rank") { | |
| 259 selectVals <- gsub(" ", "", selectVals, fixed=TRUE) | |
| 260 selectVals <- unlist(strsplit(selectVals, ",")) | |
| 261 | |
| 262 for (i in 1:length(selectVals)) { | |
| 263 if (grepl(":", selectVals[i], fixed=TRUE)) { | |
| 264 temp <- unlist(strsplit(selectVals[i], ":")) | |
| 265 selectVals <- selectVals[-i] | |
| 266 a <- as.numeric(temp[1]) | |
| 267 b <- as.numeric(temp[2]) | |
| 268 selectVals <- c(selectVals, a:b) | |
| 269 } | |
| 270 } | |
| 271 selectVals <- as.numeric(unique(selectVals)) | |
| 272 } else { | |
| 273 selectVals <- gsub(" ", "", selectVals, fixed=TRUE) | |
| 274 selectVals <- unlist(strsplit(selectVals, " ")) | |
| 275 } | |
| 276 } | |
| 277 | |
| 278 if (inputType=="fastq") { | |
| 279 # Use EdgeR hairpin process and capture outputs | |
| 280 hpReadout <- capture.output( | |
| 281 data <- processHairpinReads(fastqPath, samplePath, hairpinPath, | |
| 282 hairpinStart=hpStart, hairpinEnd=hpEnd, | |
| 283 verbose=TRUE) | |
| 284 ) | |
| 285 | |
| 286 # Remove function output entries that show processing data or is empty | |
| 287 hpReadout <- hpReadout[hpReadout!=""] | |
| 288 hpReadout <- hpReadout[!grepl("Processing", hpReadout)] | |
| 289 hpReadout <- hpReadout[!grepl("in file", hpReadout)] | |
| 290 hpReadout <- gsub(" -- ", "", hpReadout, fixed=TRUE) | |
| 291 | |
| 292 | |
| 293 # Make the names of groups syntactically valid (replace spaces with periods) | |
| 294 data$samples$group <- make.names(data$samples$group) | |
|
4
f8af57d6f60b
Fixed bug causing fastq input to break
shian_su <registertonysu@gmail.com>
parents:
2
diff
changeset
|
295 } else if (inputType=="counts") { |
| 2 | 296 # Process counts information, set ID column to be row names |
| 297 rownames(counts) <- counts$ID | |
| 298 counts <- counts[ , !(colnames(counts)=="ID")] | |
| 299 countsRows <- nrow(counts) | |
| 300 | |
| 301 # Process group information | |
| 302 factors <- samples$group[match(samples$ID, colnames(counts))] | |
| 303 annoRows <- nrow(anno) | |
| 304 anno <- anno[match(rownames(counts), anno$ID), ] | |
| 305 annoMatched <- sum(!is.na(anno$ID)) | |
| 306 | |
| 307 if (any(is.na(anno$ID))) { | |
| 308 warningStr <- paste("count table contained more hairpins than", | |
| 309 "specified in hairpin annotation file") | |
| 310 warning(warningStr) | |
| 311 } | |
| 312 | |
| 313 # Filter out rows with zero counts | |
| 314 sel <- rowSums(counts)!=0 | |
| 315 counts <- counts[sel, ] | |
| 316 anno <- anno[sel, ] | |
| 317 | |
| 318 # Create DGEList | |
| 319 data <- DGEList(counts=counts, lib.size=colSums(counts), | |
| 320 norm.factors=rep(1,ncol(counts)), genes=anno, group=factors) | |
| 321 | |
| 322 # Make the names of groups syntactically valid (replace spaces with periods) | |
| 323 data$samples$group <- make.names(data$samples$group) | |
| 324 } | |
| 325 | |
| 326 # Filter hairpins with low counts | |
| 327 sel <- rowSums(cpm(data$counts) > cpmReq) >= sampleReq | |
| 328 data <- data[sel, ] | |
| 329 | |
| 330 # Estimate dispersions | |
| 331 data <- estimateDisp(data) | |
| 332 commonBCV <- sqrt(data$common.dispersion) | |
| 333 | |
| 334 ################################################################################ | |
| 335 ### Output Processing | |
| 336 ################################################################################ | |
| 337 | |
| 338 # Plot number of hairpins that could be matched per sample | |
| 339 png(barIndexPng, width=600, height=600) | |
| 340 barplot(height<-colSums(data$counts), las=2, main="Counts per index", | |
| 341 cex.names=1.0, cex.axis=0.8, ylim=c(0, max(height)*1.2)) | |
| 342 imageData[1, ] <- c("Counts per Index", "barIndex.png") | |
| 343 invisible(dev.off()) | |
| 344 | |
| 345 pdf(barIndexPdf) | |
| 346 barplot(height<-colSums(data$counts), las=2, main="Counts per index", | |
| 347 cex.names=1.0, cex.axis=0.8, ylim=c(0, max(height)*1.2)) | |
| 348 linkData[1, ] <- c("Counts per Index Barplot (.pdf)", "barIndex.pdf") | |
| 349 invisible(dev.off()) | |
| 350 | |
| 351 # Plot per hairpin totals across all samples | |
| 352 png(barHairpinPng, width=600, height=600) | |
| 353 if (nrow(data$counts)<50) { | |
| 354 barplot(height<-rowSums(data$counts), las=2, main="Counts per hairpin", | |
| 355 cex.names=0.8, cex.axis=0.8, ylim=c(0, max(height)*1.2)) | |
| 356 } else { | |
| 357 barplot(height<-rowSums(data$counts), las=2, main="Counts per hairpin", | |
| 358 cex.names=0.8, cex.axis=0.8, ylim=c(0, max(height)*1.2), | |
| 359 names.arg=FALSE) | |
| 360 } | |
| 361 imageData <- rbind(imageData, c("Counts per Hairpin", "barHairpin.png")) | |
| 362 invisible(dev.off()) | |
| 363 | |
| 364 pdf(barHairpinPdf) | |
| 365 if (nrow(data$counts)<50) { | |
| 366 barplot(height<-rowSums(data$counts), las=2, main="Counts per hairpin", | |
| 367 cex.names=0.8, cex.axis=0.8, ylim=c(0, max(height)*1.2)) | |
| 368 } else { | |
| 369 barplot(height<-rowSums(data$counts), las=2, main="Counts per hairpin", | |
| 370 cex.names=0.8, cex.axis=0.8, ylim=c(0, max(height)*1.2), | |
| 371 names.arg=FALSE) | |
| 372 } | |
| 373 newEntry <- c("Counts per Hairpin Barplot (.pdf)", "barHairpin.pdf") | |
| 374 linkData <- rbind(linkData, newEntry) | |
| 375 invisible(dev.off()) | |
| 376 | |
| 377 # Make an MDS plot to visualise relationships between replicate samples | |
| 378 png(mdsPng, width=600, height=600) | |
| 379 plotMDS(data, labels=data$samples$group, col=as.numeric(data$samples$group), | |
| 380 main="MDS Plot") | |
| 381 imageData <- rbind(imageData, c("MDS Plot", "mds.png")) | |
| 382 invisible(dev.off()) | |
| 383 | |
| 384 pdf(mdsPdf) | |
| 385 plotMDS(data, labels=data$samples$group, col=as.numeric(data$samples$group), | |
| 386 main="MDS Plot") | |
| 387 newEntry <- c("MDS Plot (.pdf)", "mds.pdf") | |
| 388 linkData <- rbind(linkData, newEntry) | |
| 389 invisible(dev.off()) | |
| 390 | |
| 391 if (workMode=="classic") { | |
| 392 # Assess differential representation using classic exact testing methodology | |
| 393 # in edgeR | |
| 394 testData <- exactTest(data, pair=pairData) | |
| 395 | |
| 396 top <- topTags(testData, n=Inf) | |
| 397 topIDs <- top$table[(top$table$FDR < fdrThresh) & | |
| 398 (abs(top$table$logFC) > lfcThresh), 1] | |
| 399 write.table(top, file=topOut, row.names=FALSE, sep="\t") | |
| 400 linkName <- paste0("Top Tags Table(", pairData[2], "-", pairData[1], | |
| 401 ") (.tsv)") | |
| 402 linkAddr <- paste0("toptag(", pairData[2], "-", pairData[1], ").tsv") | |
| 403 linkData <- rbind(linkData, c(linkName, linkAddr)) | |
| 404 | |
| 405 # Select hairpins with FDR < 0.05 to highlight on plot | |
| 406 png(smearPng, width=600, height=600) | |
| 407 plotTitle <- gsub(".", " ", | |
| 408 paste0("Smear Plot: ", pairData[2], "-", pairData[1]), | |
| 409 fixed = TRUE) | |
| 410 plotSmear(testData, de.tags=topIDs, | |
| 411 pch=20, cex=1.0, main=plotTitle) | |
| 412 abline(h = c(-1, 0, 1), col = c("dodgerblue", "yellow", "dodgerblue"), lty=2) | |
| 413 imgName <- paste0("Smear Plot(", pairData[2], "-", pairData[1], ")") | |
| 414 imgAddr <- paste0("smear(", pairData[2], "-", pairData[1],").png") | |
| 415 imageData <- rbind(imageData, c(imgName, imgAddr)) | |
| 416 invisible(dev.off()) | |
| 417 | |
| 418 pdf(smearPdf) | |
| 419 plotTitle <- gsub(".", " ", | |
| 420 paste0("Smear Plot: ", pairData[2], "-", pairData[1]), | |
| 421 fixed = TRUE) | |
| 422 plotSmear(testData, de.tags=topIDs, | |
| 423 pch=20, cex=1.0, main=plotTitle) | |
| 424 abline(h = c(-1, 0, 1), col = c("dodgerblue", "yellow", "dodgerblue"), lty=2) | |
| 425 imgName <- paste0("Smear Plot(", pairData[2], "-", pairData[1], ") (.pdf)") | |
| 426 imgAddr <- paste0("smear(", pairData[2], "-", pairData[1], ").pdf") | |
| 427 linkData <- rbind(linkData, c(imgName, imgAddr)) | |
| 428 invisible(dev.off()) | |
| 429 } else if (workMode=="glm") { | |
| 430 # Generating design information | |
| 431 factors <- factor(data$sample$group) | |
| 432 design <- model.matrix(~0+factors) | |
| 433 | |
| 434 colnames(design) <- gsub("factors", "", colnames(design), fixed=TRUE) | |
| 435 | |
| 436 # Split up contrasts seperated by comma into a vector | |
| 437 contrastData <- unlist(strsplit(contrastData, split=",")) | |
| 438 for (i in 1:length(contrastData)) { | |
| 439 # Generate contrasts information | |
| 440 contrasts <- makeContrasts(contrasts=contrastData[i], levels=design) | |
| 441 | |
| 442 # Fit negative bionomial GLM | |
| 443 fit = glmFit(data, design) | |
| 444 # Carry out Likelihood ratio test | |
| 445 testData = glmLRT(fit, contrast=contrasts) | |
| 446 | |
| 447 # Select hairpins with FDR < 0.05 to highlight on plot | |
| 448 top <- topTags(testData, n=Inf) | |
| 449 topIDs <- top$table[(top$table$FDR < fdrThresh) & | |
| 450 (abs(top$table$logFC) > lfcThresh), 1] | |
| 451 write.table(top, file=topOut[i], row.names=FALSE, sep="\t") | |
| 452 | |
| 453 linkName <- paste0("Top Tags Table(", contrastData[i], ") (.tsv)") | |
| 454 linkAddr <- paste0("toptag(", contrastData[i], ").tsv") | |
| 455 linkData <- rbind(linkData, c(linkName, linkAddr)) | |
| 456 | |
| 457 # Make a plot of logFC versus logCPM | |
| 458 png(smearPng[i], height=600, width=600) | |
| 459 plotTitle <- paste("Smear Plot:", gsub(".", " ", contrastData[i], | |
| 460 fixed=TRUE)) | |
| 461 plotSmear(testData, de.tags=topIDs, pch=20, cex=0.8, main=plotTitle) | |
| 462 abline(h=c(-1, 0, 1), col=c("dodgerblue", "yellow", "dodgerblue"), lty=2) | |
| 463 | |
| 464 imgName <- paste0("Smear Plot(", contrastData[i], ")") | |
| 465 imgAddr <- paste0("smear(", contrastData[i], ").png") | |
| 466 imageData <- rbind(imageData, c(imgName, imgAddr)) | |
| 467 invisible(dev.off()) | |
| 468 | |
| 469 pdf(smearPdf[i]) | |
| 470 plotTitle <- paste("Smear Plot:", gsub(".", " ", contrastData[i], | |
| 471 fixed=TRUE)) | |
| 472 plotSmear(testData, de.tags=topIDs, pch=20, cex=0.8, main=plotTitle) | |
| 473 abline(h=c(-1, 0, 1), col=c("dodgerblue", "yellow", "dodgerblue"), lty=2) | |
| 474 | |
| 475 linkName <- paste0("Smear Plot(", contrastData[i], ") (.pdf)") | |
| 476 linkAddr <- paste0("smear(", contrastData[i], ").pdf") | |
| 477 linkData <- rbind(linkData, c(linkName, linkAddr)) | |
| 478 invisible(dev.off()) | |
| 479 | |
| 480 genes <- as.character(data$genes$Gene) | |
| 481 unq <- unique(genes) | |
| 482 unq <- unq[!is.na(unq)] | |
| 483 geneList <- list() | |
| 484 for (gene in unq) { | |
| 485 if (length(which(genes==gene)) >= hairpinReq) { | |
| 486 geneList[[gene]] <- which(genes==gene) | |
| 487 } | |
| 488 } | |
| 489 | |
| 490 if (wantRoast) { | |
| 491 # Input preparaton for roast | |
| 492 nrot = 9999 | |
| 493 set.seed(602214129) | |
| 494 roastData <- mroast(data, index=geneList, design=design, | |
| 495 contrast=contrasts, nrot=nrot) | |
| 496 roastData <- cbind(GeneID=rownames(roastData), roastData) | |
| 497 write.table(roastData, file=roastOut[i], row.names=FALSE, sep="\t") | |
| 498 linkName <- paste0("Gene Level Analysis Table(", contrastData[i], | |
| 499 ") (.tsv)") | |
| 500 linkAddr <- paste0("roast(", contrastData[i], ").tsv") | |
| 501 linkData <- rbind(linkData, c(linkName, linkAddr)) | |
| 502 if (selectOpt=="rank") { | |
| 503 selectedGenes <- rownames(roastData)[selectVals] | |
| 504 } else { | |
| 505 selectedGenes <- selectVals | |
| 506 } | |
| 507 | |
| 508 if (packageVersion("limma")<"3.19.19") { | |
| 509 png(barcodePng[i], width=600, height=length(selectedGenes)*150) | |
| 510 } else { | |
| 511 png(barcodePng[i], width=600, height=length(selectedGenes)*300) | |
| 512 } | |
| 513 par(mfrow=c(length(selectedGenes), 1)) | |
| 514 for (gene in selectedGenes) { | |
| 515 barcodeplot(testData$table$logFC, index=geneList[[gene]], | |
| 516 main=paste("Barcode Plot for", gene, "(logFCs)", | |
| 517 gsub(".", " ", contrastData[i])), | |
| 518 labels=c("Positive logFC", "Negative logFC")) | |
| 519 } | |
| 520 imgName <- paste0("Barcode Plot(", contrastData[i], ")") | |
| 521 imgAddr <- paste0("barcode(", contrastData[i], ").png") | |
| 522 imageData <- rbind(imageData, c(imgName, imgAddr)) | |
| 523 dev.off() | |
| 524 if (packageVersion("limma")<"3.19.19") { | |
| 525 pdf(barcodePdf[i], width=8, height=2) | |
| 526 } else { | |
| 527 pdf(barcodePdf[i], width=8, height=4) | |
| 528 } | |
| 529 for (gene in selectedGenes) { | |
| 530 barcodeplot(testData$table$logFC, index=geneList[[gene]], | |
| 531 main=paste("Barcode Plot for", gene, "(logFCs)", | |
| 532 gsub(".", " ", contrastData[i])), | |
| 533 labels=c("Positive logFC", "Negative logFC")) | |
| 534 } | |
| 535 linkName <- paste0("Barcode Plot(", contrastData[i], ") (.pdf)") | |
| 536 linkAddr <- paste0("barcode(", contrastData[i], ").pdf") | |
| 537 linkData <- rbind(linkData, c(linkName, linkAddr)) | |
| 538 dev.off() | |
| 539 } | |
| 540 } | |
| 541 } | |
| 542 | |
| 543 # Record ending time | |
| 544 timeEnd <- as.character(Sys.time()) | |
| 545 ################################################################################ | |
| 546 ### HTML Generation | |
| 547 ################################################################################ | |
| 548 # Clear file | |
| 549 cat("", file=htmlPath) | |
| 550 | |
| 551 cata("<html>\n") | |
| 552 HtmlHead("EdgeR Output") | |
| 553 | |
| 554 cata("<body>\n") | |
| 555 cata("<h3>EdgeR Analysis Output:</h3>\n") | |
| 556 cata("<h4>Input Summary:</h4>\n") | |
| 557 if (inputType=="fastq") { | |
| 558 cata("<ul>\n") | |
| 559 ListItem(hpReadout[1]) | |
| 560 ListItem(hpReadout[2]) | |
| 561 cata("</ul>\n") | |
| 562 cata(hpReadout[3], "<br/>\n") | |
| 563 cata("<ul>\n") | |
| 564 ListItem(hpReadout[4]) | |
| 565 ListItem(hpReadout[7]) | |
| 566 cata("</ul>\n") | |
| 567 cata(hpReadout[8:11], sep="<br/>\n") | |
| 568 cata("<br />\n") | |
| 569 cata("<b>Please check that read percentages are consistent with ") | |
| 570 cata("expectations.</b><br >\n") | |
| 571 } else if (inputType=="counts") { | |
| 572 cata("<ul>\n") | |
| 573 ListItem("Number of Samples: ", ncol(data$counts)) | |
| 574 ListItem("Number of Hairpins: ", countsRows) | |
| 575 ListItem("Number of annotations provided: ", annoRows) | |
| 576 ListItem("Number of annotations matched to hairpin: ", annoMatched) | |
| 577 cata("</ul>\n") | |
| 578 } | |
| 579 | |
| 580 cata("The estimated common biological coefficient of variation (BCV) is: ", | |
| 581 commonBCV, "<br />\n") | |
| 582 | |
| 583 cata("<h4>Output:</h4>\n") | |
| 584 cata("All images displayed have PDF copy at the bottom of the page, these can ") | |
| 585 cata("exported in a pdf viewer to high resolution image format. <br/>\n") | |
| 586 for (i in 1:nrow(imageData)) { | |
| 587 if (grepl("barcode", imageData$Link[i])) { | |
| 588 if (packageVersion("limma")<"3.19.19") { | |
| 589 HtmlImage(imageData$Link[i], imageData$Label[i], | |
| 590 height=length(selectedGenes)*150) | |
| 591 } else { | |
| 592 HtmlImage(imageData$Link[i], imageData$Label[i], | |
| 593 height=length(selectedGenes)*300) | |
| 594 } | |
| 595 } else { | |
| 596 HtmlImage(imageData$Link[i], imageData$Label[i]) | |
| 597 } | |
| 598 } | |
| 599 cata("<br/>\n") | |
| 600 | |
| 601 cata("<h4>Plots:</h4>\n") | |
| 602 for (i in 1:nrow(linkData)) { | |
| 603 if (!grepl(".tsv", linkData$Link[i])) { | |
| 604 HtmlLink(linkData$Link[i], linkData$Label[i]) | |
| 605 } | |
| 606 } | |
| 607 | |
| 608 cata("<h4>Tables:</h4>\n") | |
| 609 for (i in 1:nrow(linkData)) { | |
| 610 if (grepl(".tsv", linkData$Link[i])) { | |
| 611 HtmlLink(linkData$Link[i], linkData$Label[i]) | |
| 612 } | |
| 613 } | |
| 614 | |
| 615 cata("<p>alt-click any of the links to download the file, or click the name ") | |
| 616 cata("of this task in the galaxy history panel and click on the floppy ") | |
| 617 cata("disk icon to download all files in a zip archive.</p>\n") | |
| 618 cata("<p>.tsv files are tab seperated files that can be viewed using Excel ") | |
| 619 cata("or other spreadsheet programs</p>\n") | |
| 620 cata("<table border=\"0\">\n") | |
| 621 | |
| 622 cata("<tr>\n") | |
| 623 TableItem("Task started at:"); TableItem(timeStart) | |
| 624 cata("</tr>\n") | |
| 625 cata("<tr>\n") | |
| 626 TableItem("Task ended at:"); TableItem(timeEnd) | |
| 627 cata("</tr>\n") | |
| 628 | |
| 629 cata("</body>\n") | |
| 630 cata("</html>") |