|
1
|
1 # ARGS: 1.inputType -String specifying format of input (fastq or table)
|
|
|
2 # IF inputType is "fastQ":
|
|
|
3 # 2*.fastqPath -One or more strings specifying path to fastq files
|
|
|
4 # 2.annoPath -String specifying path to hairpin annotation table
|
|
|
5 # 3.samplePath -String specifying path to sample annotation table
|
|
|
6 # 4.barStart -Integer specifying starting position of barcode
|
|
|
7 # 5.barEnd -Integer specifying ending position of barcode
|
|
|
8 # 6.hpStart -Integer specifying startins position of hairpin
|
|
|
9 # unique region
|
|
|
10 # 7.hpEnd -Integer specifying ending position of hairpin
|
|
|
11 # unique region
|
|
|
12 # ###
|
|
|
13 # IF inputType is "counts":
|
|
|
14 # 2.countPath -String specifying path to count table
|
|
|
15 # 3.annoPath -String specifying path to hairpin annotation table
|
|
|
16 # 4.samplePath -String specifying path to sample annotation table
|
|
|
17 # ###
|
|
|
18 # 8.cpmReq -Float specifying cpm requirement
|
|
|
19 # 9.sampleReq -Integer specifying cpm requirement
|
|
|
20 # 10.fdrThresh -Float specifying the FDR requirement
|
|
|
21 # 11.lfcThresh -Float specifying the log-fold-change requirement
|
|
|
22 # 12.workMode -String specifying exact test or GLM usage
|
|
|
23 # 13.htmlPath -String specifying path to HTML file
|
|
|
24 # 14.folderPath -STring specifying path to folder for output
|
|
|
25 # IF workMode is "classic" (exact test)
|
|
|
26 # 15.pairData[2] -String specifying first group for exact test
|
|
|
27 # 16.pairData[1] -String specifying second group for exact test
|
|
|
28 # ###
|
|
|
29 # IF workMode is "glm"
|
|
|
30 # 15.contrastData -String specifying contrasts to be made
|
|
|
31 # 16.roastOpt -String specifying usage of gene-wise tests
|
|
|
32 # 17.hairpinReq -String specifying hairpin requirement for gene-
|
|
|
33 # wise test
|
|
|
34 # 18.selectOpt -String specifying type of selection for barcode
|
|
|
35 # plots
|
|
|
36 # 19.selectVals -String specifying members selected for barcode
|
|
|
37 # plots
|
|
|
38 #
|
|
|
39 # OUT: Bar Plot of Counts Per Index
|
|
|
40 # Bar Plot of Counts Per Hairpin
|
|
|
41 # MDS Plot
|
|
|
42 # Smear Plot
|
|
|
43 # Barcode Plots (If Genewise testing was selected)
|
|
|
44 # Top Expression Table
|
|
3
|
45 # Feature Counts Table
|
|
1
|
46 # HTML file linking to the ouputs
|
|
|
47 #
|
|
|
48 # Author: Shian Su - registertonysu@gmail.com - Jan 2014
|
|
|
49
|
|
0
|
50 # Record starting time
|
|
|
51 timeStart <- as.character(Sys.time())
|
|
|
52
|
|
|
53 # Loading and checking required packages
|
|
|
54 library(methods, quietly=TRUE, warn.conflicts=FALSE)
|
|
|
55 library(statmod, quietly=TRUE, warn.conflicts=FALSE)
|
|
|
56 library(splines, quietly=TRUE, warn.conflicts=FALSE)
|
|
|
57 library(edgeR, quietly=TRUE, warn.conflicts=FALSE)
|
|
|
58 library(limma, quietly=TRUE, warn.conflicts=FALSE)
|
|
|
59
|
|
3
|
60 if (packageVersion("edgeR") < "3.5.27") {
|
|
|
61 stop("Please update 'edgeR' to version >= 3.5.23 to run this tool")
|
|
|
62 }
|
|
|
63
|
|
|
64 if (packageVersion("limma")<"3.19.19") {
|
|
|
65 message("Update 'limma' to version >= 3.19.19 to see updated barcode graphs")
|
|
0
|
66 }
|
|
|
67
|
|
|
68 ################################################################################
|
|
|
69 ### Function declarations
|
|
|
70 ################################################################################
|
|
|
71
|
|
3
|
72 # Function to load libaries without messages
|
|
|
73 silentLibrary <- function(...) {
|
|
|
74 list <- c(...)
|
|
|
75 for (package in list){
|
|
|
76 suppressPackageStartupMessages(library(package, character.only=TRUE))
|
|
|
77 }
|
|
|
78 }
|
|
|
79
|
|
0
|
80 # Function to sanitise contrast equations so there are no whitespaces
|
|
|
81 # surrounding the arithmetic operators, leading or trailing whitespace
|
|
|
82 sanitiseEquation <- function(equation) {
|
|
|
83 equation <- gsub(" *[+] *", "+", equation)
|
|
|
84 equation <- gsub(" *[-] *", "-", equation)
|
|
|
85 equation <- gsub(" *[/] *", "/", equation)
|
|
|
86 equation <- gsub(" *[*] *", "*", equation)
|
|
|
87 equation <- gsub("^\\s+|\\s+$", "", equation)
|
|
|
88 return(equation)
|
|
|
89 }
|
|
|
90
|
|
|
91 # Function to sanitise group information
|
|
|
92 sanitiseGroups <- function(string) {
|
|
|
93 string <- gsub(" *[,] *", ",", string)
|
|
|
94 string <- gsub("^\\s+|\\s+$", "", string)
|
|
|
95 return(string)
|
|
|
96 }
|
|
|
97
|
|
|
98 # Function to change periods to whitespace in a string
|
|
|
99 unmake.names <- function(string) {
|
|
|
100 string <- gsub(".", " ", string, fixed=TRUE)
|
|
|
101 return(string)
|
|
|
102 }
|
|
|
103
|
|
|
104 # Function has string input and generates an output path string
|
|
|
105 makeOut <- function(filename) {
|
|
|
106 return(paste0(folderPath, "/", filename))
|
|
|
107 }
|
|
|
108
|
|
|
109 # Function has string input and generates both a pdf and png output strings
|
|
|
110 imgOut <- function(filename) {
|
|
|
111 assign(paste0(filename, "Png"), makeOut(paste0(filename,".png")),
|
|
3
|
112 envir=.GlobalEnv)
|
|
0
|
113 assign(paste0(filename, "Pdf"), makeOut(paste0(filename,".pdf")),
|
|
3
|
114 envir=.GlobalEnv)
|
|
0
|
115 }
|
|
|
116
|
|
|
117 # Create cat function default path set, default seperator empty and appending
|
|
|
118 # true by default (Ripped straight from the cat function with altered argument
|
|
|
119 # defaults)
|
|
3
|
120 cata <- function(..., file=htmlPath, sep="", fill=FALSE, labels=NULL,
|
|
|
121 append=TRUE) {
|
|
0
|
122 if (is.character(file))
|
|
|
123 if (file == "")
|
|
|
124 file <- stdout()
|
|
|
125 else if (substring(file, 1L, 1L) == "|") {
|
|
|
126 file <- pipe(substring(file, 2L), "w")
|
|
|
127 on.exit(close(file))
|
|
|
128 }
|
|
|
129 else {
|
|
|
130 file <- file(file, ifelse(append, "a", "w"))
|
|
|
131 on.exit(close(file))
|
|
|
132 }
|
|
|
133 .Internal(cat(list(...), file, sep, fill, labels, append))
|
|
|
134 }
|
|
|
135
|
|
|
136 # Function to write code for html head and title
|
|
|
137 HtmlHead <- function(title) {
|
|
|
138 cata("<head>\n")
|
|
|
139 cata("<title>", title, "</title>\n")
|
|
|
140 cata("</head>\n")
|
|
|
141 }
|
|
|
142
|
|
|
143 # Function to write code for html links
|
|
|
144 HtmlLink <- function(address, label=address) {
|
|
|
145 cata("<a href=\"", address, "\" target=\"_blank\">", label, "</a><br />\n")
|
|
|
146 }
|
|
|
147
|
|
|
148 # Function to write code for html images
|
|
|
149 HtmlImage <- function(source, label=source, height=600, width=600) {
|
|
|
150 cata("<img src=\"", source, "\" alt=\"", label, "\" height=\"", height)
|
|
|
151 cata("\" width=\"", width, "\"/>\n")
|
|
|
152 }
|
|
|
153
|
|
|
154 # Function to write code for html list items
|
|
|
155 ListItem <- function(...) {
|
|
|
156 cata("<li>", ..., "</li>\n")
|
|
|
157 }
|
|
|
158
|
|
|
159 TableItem <- function(...) {
|
|
|
160 cata("<td>", ..., "</td>\n")
|
|
|
161 }
|
|
|
162
|
|
|
163 TableHeadItem <- function(...) {
|
|
|
164 cata("<th>", ..., "</th>\n")
|
|
|
165 }
|
|
|
166 ################################################################################
|
|
|
167 ### Input Processing
|
|
|
168 ################################################################################
|
|
|
169
|
|
|
170 # Grabbing arguments from command line
|
|
|
171 argv <- commandArgs(TRUE)
|
|
|
172
|
|
|
173 # Remove fastq file paths after collecting from argument vector
|
|
|
174 inputType <- as.character(argv[1])
|
|
|
175 if (inputType=="fastq") {
|
|
|
176 fastqPath <- as.character(gsub("fastq::", "", argv[grepl("fastq::", argv)],
|
|
|
177 fixed=TRUE))
|
|
|
178 argv <- argv[!grepl("fastq::", argv, fixed=TRUE)]
|
|
1
|
179 annoPath <- as.character(argv[2])
|
|
0
|
180 samplePath <- as.character(argv[3])
|
|
|
181 barStart <- as.numeric(argv[4])
|
|
|
182 barEnd <- as.numeric(argv[5])
|
|
|
183 hpStart <- as.numeric(argv[6])
|
|
|
184 hpEnd <- as.numeric(argv[7])
|
|
|
185 } else if (inputType=="counts") {
|
|
|
186 countPath <- as.character(argv[2])
|
|
|
187 annoPath <- as.character(argv[3])
|
|
|
188 samplePath <- as.character(argv[4])
|
|
|
189 }
|
|
|
190
|
|
|
191 cpmReq <- as.numeric(argv[8])
|
|
|
192 sampleReq <- as.numeric(argv[9])
|
|
|
193 fdrThresh <- as.numeric(argv[10])
|
|
|
194 lfcThresh <- as.numeric(argv[11])
|
|
|
195 workMode <- as.character(argv[12])
|
|
|
196 htmlPath <- as.character(argv[13])
|
|
|
197 folderPath <- as.character(argv[14])
|
|
3
|
198
|
|
0
|
199 if (workMode=="classic") {
|
|
|
200 pairData <- character()
|
|
|
201 pairData[2] <- as.character(argv[15])
|
|
|
202 pairData[1] <- as.character(argv[16])
|
|
|
203 } else if (workMode=="glm") {
|
|
|
204 contrastData <- as.character(argv[15])
|
|
|
205 roastOpt <- as.character(argv[16])
|
|
|
206 hairpinReq <- as.numeric(argv[17])
|
|
|
207 selectOpt <- as.character(argv[18])
|
|
|
208 selectVals <- as.character(argv[19])
|
|
|
209 }
|
|
|
210
|
|
|
211 # Read in inputs
|
|
1
|
212
|
|
|
213 samples <- read.table(samplePath, header=TRUE, sep="\t")
|
|
|
214 anno <- read.table(annoPath, header=TRUE, sep="\t")
|
|
|
215 if (inputType=="counts") {
|
|
0
|
216 counts <- read.table(countPath, header=TRUE, sep="\t")
|
|
|
217 }
|
|
1
|
218
|
|
0
|
219 ###################### Check inputs for correctness ############################
|
|
|
220 samples$ID <- make.names(samples$ID)
|
|
|
221
|
|
|
222 if (!any(grepl("group", names(samples)))) {
|
|
|
223 stop("'group' column not specified in sample annotation file")
|
|
|
224 } # Check if grouping variable has been specified
|
|
|
225
|
|
|
226 if (any(table(samples$ID)>1)){
|
|
|
227 tab <- table(samples$ID)
|
|
|
228 offenders <- paste(names(tab[tab>1]), collapse=", ")
|
|
|
229 offenders <- unmake.names(offenders)
|
|
1
|
230 stop("'ID' column of sample annotation must have unique values, values ",
|
|
0
|
231 offenders, " are repeated")
|
|
|
232 } # Check that IDs in sample annotation are unique
|
|
|
233
|
|
|
234 if (inputType=="fastq") {
|
|
1
|
235
|
|
|
236 if (any(table(anno$ID)>1)){
|
|
|
237 tab <- table(anno$ID)
|
|
0
|
238 offenders <- paste(names(tab[tab>1]), collapse=", ")
|
|
1
|
239 stop("'ID' column of hairpin annotation must have unique values, values ",
|
|
0
|
240 offenders, " are repeated")
|
|
|
241 } # Check that IDs in hairpin annotation are unique
|
|
1
|
242
|
|
0
|
243 } else if (inputType=="counts") {
|
|
1
|
244 if (any(is.na(match(samples$ID, colnames(counts))))) {
|
|
|
245 stop("not all samples have groups specified")
|
|
|
246 } # Check that a group has be specifed for each sample
|
|
0
|
247
|
|
|
248 if (any(table(counts$ID)>1)){
|
|
|
249 tab <- table(counts$ID)
|
|
|
250 offenders <- paste(names(tab[tab>1]), collapse=", ")
|
|
1
|
251 stop("'ID' column of count table must have unique values, values ",
|
|
0
|
252 offenders, " are repeated")
|
|
|
253 } # Check that IDs in count table are unique
|
|
|
254 }
|
|
1
|
255 if (workMode=="glm") {
|
|
|
256 if (roastOpt == "yes") {
|
|
|
257 if (is.na(match("Gene", colnames(anno)))) {
|
|
|
258 tempStr <- paste("Gene-wise tests selected but'Gene' column not",
|
|
|
259 "specified in hairpin annotation file")
|
|
|
260 stop(tempStr)
|
|
|
261 }
|
|
|
262 }
|
|
|
263 }
|
|
|
264
|
|
0
|
265 ################################################################################
|
|
|
266
|
|
|
267 # Process arguments
|
|
|
268 if (workMode=="glm") {
|
|
|
269 if (roastOpt=="yes") {
|
|
|
270 wantRoast <- TRUE
|
|
|
271 } else {
|
|
|
272 wantRoast <- FALSE
|
|
|
273 }
|
|
|
274 }
|
|
|
275
|
|
|
276 # Split up contrasts seperated by comma into a vector and replace spaces with
|
|
|
277 # periods
|
|
|
278 if (exists("contrastData")) {
|
|
|
279 contrastData <- unlist(strsplit(contrastData, split=","))
|
|
|
280 contrastData <- sanitiseEquation(contrastData)
|
|
|
281 contrastData <- gsub(" ", ".", contrastData, fixed=TRUE)
|
|
|
282 }
|
|
|
283
|
|
|
284 # Replace spaces with periods in pair data
|
|
|
285 if (exists("pairData")) {
|
|
|
286 pairData <- make.names(pairData)
|
|
|
287 }
|
|
|
288
|
|
|
289 # Generate output folder and paths
|
|
3
|
290 dir.create(folderPath, showWarnings=FALSE)
|
|
0
|
291
|
|
|
292 # Generate links for outputs
|
|
|
293 imgOut("barHairpin")
|
|
|
294 imgOut("barIndex")
|
|
|
295 imgOut("mds")
|
|
|
296 imgOut("bcv")
|
|
|
297 if (workMode == "classic") {
|
|
|
298 smearPng <- makeOut(paste0("smear(", pairData[2], "-", pairData[1],").png"))
|
|
|
299 smearPdf <- makeOut(paste0("smear(", pairData[2], "-", pairData[1],").pdf"))
|
|
|
300 topOut <- makeOut(paste0("toptag(", pairData[2], "-", pairData[1],").tsv"))
|
|
|
301 } else if (workMode=="glm") {
|
|
|
302 smearPng <- character()
|
|
|
303 smearPdf <- character()
|
|
|
304 topOut <- character()
|
|
|
305 roastOut <- character()
|
|
|
306 barcodePng <- character()
|
|
|
307 barcodePdf <- character()
|
|
|
308 for (i in 1:length(contrastData)) {
|
|
|
309 smearPng[i] <- makeOut(paste0("smear(", contrastData[i], ").png"))
|
|
|
310 smearPdf[i] <- makeOut(paste0("smear(", contrastData[i], ").pdf"))
|
|
|
311 topOut[i] <- makeOut(paste0("toptag(", contrastData[i], ").tsv"))
|
|
3
|
312 roastOut[i] <- makeOut(paste0("gene_level(", contrastData[i], ").tsv"))
|
|
0
|
313 barcodePng[i] <- makeOut(paste0("barcode(", contrastData[i], ").png"))
|
|
|
314 barcodePdf[i] <- makeOut(paste0("barcode(", contrastData[i], ").pdf"))
|
|
|
315 }
|
|
|
316 }
|
|
3
|
317 countsOut <- makeOut("counts.tsv")
|
|
|
318 sessionOut <- makeOut("session_info.txt")
|
|
|
319
|
|
0
|
320 # Initialise data for html links and images, table with the link label and
|
|
|
321 # link address
|
|
|
322 linkData <- data.frame(Label=character(), Link=character(),
|
|
|
323 stringsAsFactors=FALSE)
|
|
|
324 imageData <- data.frame(Label=character(), Link=character(),
|
|
|
325 stringsAsFactors=FALSE)
|
|
3
|
326
|
|
|
327 # Initialise vectors for storage of up/down/neutral regulated counts
|
|
|
328 upCount <- numeric()
|
|
|
329 downCount <- numeric()
|
|
|
330 flatCount <- numeric()
|
|
|
331
|
|
0
|
332 ################################################################################
|
|
|
333 ### Data Processing
|
|
|
334 ################################################################################
|
|
|
335
|
|
|
336 # Transform gene selection from string into index values for mroast
|
|
|
337 if (workMode=="glm") {
|
|
|
338 if (selectOpt=="rank") {
|
|
|
339 selectVals <- gsub(" ", "", selectVals, fixed=TRUE)
|
|
|
340 selectVals <- unlist(strsplit(selectVals, ","))
|
|
|
341
|
|
|
342 for (i in 1:length(selectVals)) {
|
|
|
343 if (grepl(":", selectVals[i], fixed=TRUE)) {
|
|
|
344 temp <- unlist(strsplit(selectVals[i], ":"))
|
|
|
345 selectVals <- selectVals[-i]
|
|
|
346 a <- as.numeric(temp[1])
|
|
|
347 b <- as.numeric(temp[2])
|
|
|
348 selectVals <- c(selectVals, a:b)
|
|
|
349 }
|
|
|
350 }
|
|
|
351 selectVals <- as.numeric(unique(selectVals))
|
|
|
352 } else {
|
|
|
353 selectVals <- gsub(" ", "", selectVals, fixed=TRUE)
|
|
3
|
354 selectVals <- unlist(strsplit(selectVals, ","))
|
|
0
|
355 }
|
|
|
356 }
|
|
|
357
|
|
|
358 if (inputType=="fastq") {
|
|
1
|
359 # Use EdgeR hairpin process and capture outputs
|
|
|
360 hpReadout <- capture.output(
|
|
|
361 data <- processHairpinReads(fastqPath, samplePath, annoPath,
|
|
0
|
362 hairpinStart=hpStart, hairpinEnd=hpEnd,
|
|
|
363 verbose=TRUE)
|
|
1
|
364 )
|
|
|
365
|
|
|
366 # Remove function output entries that show processing data or is empty
|
|
|
367 hpReadout <- hpReadout[hpReadout!=""]
|
|
|
368 hpReadout <- hpReadout[!grepl("Processing", hpReadout)]
|
|
|
369 hpReadout <- hpReadout[!grepl("in file", hpReadout)]
|
|
|
370 hpReadout <- gsub(" -- ", "", hpReadout, fixed=TRUE)
|
|
|
371
|
|
|
372 # Make the names of groups syntactically valid (replace spaces with periods)
|
|
|
373 data$samples$group <- make.names(data$samples$group)
|
|
|
374 } else if (inputType=="counts") {
|
|
0
|
375 # Process counts information, set ID column to be row names
|
|
|
376 rownames(counts) <- counts$ID
|
|
|
377 counts <- counts[ , !(colnames(counts)=="ID")]
|
|
|
378 countsRows <- nrow(counts)
|
|
|
379
|
|
|
380 # Process group information
|
|
|
381 factors <- samples$group[match(samples$ID, colnames(counts))]
|
|
|
382 annoRows <- nrow(anno)
|
|
|
383 anno <- anno[match(rownames(counts), anno$ID), ]
|
|
|
384 annoMatched <- sum(!is.na(anno$ID))
|
|
|
385
|
|
|
386 if (any(is.na(anno$ID))) {
|
|
|
387 warningStr <- paste("count table contained more hairpins than",
|
|
|
388 "specified in hairpin annotation file")
|
|
|
389 warning(warningStr)
|
|
|
390 }
|
|
|
391
|
|
|
392 # Filter out rows with zero counts
|
|
|
393 sel <- rowSums(counts)!=0
|
|
|
394 counts <- counts[sel, ]
|
|
|
395 anno <- anno[sel, ]
|
|
|
396
|
|
|
397 # Create DGEList
|
|
|
398 data <- DGEList(counts=counts, lib.size=colSums(counts),
|
|
|
399 norm.factors=rep(1,ncol(counts)), genes=anno, group=factors)
|
|
1
|
400
|
|
0
|
401 # Make the names of groups syntactically valid (replace spaces with periods)
|
|
|
402 data$samples$group <- make.names(data$samples$group)
|
|
|
403 }
|
|
|
404
|
|
|
405 # Filter hairpins with low counts
|
|
1
|
406 preFilterCount <- nrow(data)
|
|
0
|
407 sel <- rowSums(cpm(data$counts) > cpmReq) >= sampleReq
|
|
|
408 data <- data[sel, ]
|
|
1
|
409 postFilterCount <- nrow(data)
|
|
|
410 filteredCount <- preFilterCount-postFilterCount
|
|
0
|
411
|
|
|
412 # Estimate dispersions
|
|
|
413 data <- estimateDisp(data)
|
|
|
414 commonBCV <- sqrt(data$common.dispersion)
|
|
|
415
|
|
|
416 ################################################################################
|
|
|
417 ### Output Processing
|
|
|
418 ################################################################################
|
|
|
419
|
|
|
420 # Plot number of hairpins that could be matched per sample
|
|
|
421 png(barIndexPng, width=600, height=600)
|
|
|
422 barplot(height<-colSums(data$counts), las=2, main="Counts per index",
|
|
|
423 cex.names=1.0, cex.axis=0.8, ylim=c(0, max(height)*1.2))
|
|
|
424 imageData[1, ] <- c("Counts per Index", "barIndex.png")
|
|
|
425 invisible(dev.off())
|
|
|
426
|
|
|
427 pdf(barIndexPdf)
|
|
|
428 barplot(height<-colSums(data$counts), las=2, main="Counts per index",
|
|
|
429 cex.names=1.0, cex.axis=0.8, ylim=c(0, max(height)*1.2))
|
|
|
430 linkData[1, ] <- c("Counts per Index Barplot (.pdf)", "barIndex.pdf")
|
|
|
431 invisible(dev.off())
|
|
|
432
|
|
|
433 # Plot per hairpin totals across all samples
|
|
|
434 png(barHairpinPng, width=600, height=600)
|
|
|
435 if (nrow(data$counts)<50) {
|
|
|
436 barplot(height<-rowSums(data$counts), las=2, main="Counts per hairpin",
|
|
|
437 cex.names=0.8, cex.axis=0.8, ylim=c(0, max(height)*1.2))
|
|
|
438 } else {
|
|
|
439 barplot(height<-rowSums(data$counts), las=2, main="Counts per hairpin",
|
|
|
440 cex.names=0.8, cex.axis=0.8, ylim=c(0, max(height)*1.2),
|
|
|
441 names.arg=FALSE)
|
|
|
442 }
|
|
|
443 imageData <- rbind(imageData, c("Counts per Hairpin", "barHairpin.png"))
|
|
|
444 invisible(dev.off())
|
|
|
445
|
|
|
446 pdf(barHairpinPdf)
|
|
|
447 if (nrow(data$counts)<50) {
|
|
|
448 barplot(height<-rowSums(data$counts), las=2, main="Counts per hairpin",
|
|
|
449 cex.names=0.8, cex.axis=0.8, ylim=c(0, max(height)*1.2))
|
|
|
450 } else {
|
|
|
451 barplot(height<-rowSums(data$counts), las=2, main="Counts per hairpin",
|
|
|
452 cex.names=0.8, cex.axis=0.8, ylim=c(0, max(height)*1.2),
|
|
|
453 names.arg=FALSE)
|
|
|
454 }
|
|
|
455 newEntry <- c("Counts per Hairpin Barplot (.pdf)", "barHairpin.pdf")
|
|
|
456 linkData <- rbind(linkData, newEntry)
|
|
|
457 invisible(dev.off())
|
|
|
458
|
|
|
459 # Make an MDS plot to visualise relationships between replicate samples
|
|
|
460 png(mdsPng, width=600, height=600)
|
|
|
461 plotMDS(data, labels=data$samples$group, col=as.numeric(data$samples$group),
|
|
|
462 main="MDS Plot")
|
|
|
463 imageData <- rbind(imageData, c("MDS Plot", "mds.png"))
|
|
|
464 invisible(dev.off())
|
|
|
465
|
|
|
466 pdf(mdsPdf)
|
|
|
467 plotMDS(data, labels=data$samples$group, col=as.numeric(data$samples$group),
|
|
|
468 main="MDS Plot")
|
|
|
469 newEntry <- c("MDS Plot (.pdf)", "mds.pdf")
|
|
|
470 linkData <- rbind(linkData, newEntry)
|
|
|
471 invisible(dev.off())
|
|
|
472
|
|
3
|
473 # BCV Plot
|
|
|
474 png(bcvPng, width=600, height=600)
|
|
|
475 plotBCV(data, main="BCV Plot")
|
|
|
476 imageData <- rbind(imageData, c("BCV Plot", "bcv.png"))
|
|
|
477 invisible(dev.off())
|
|
|
478
|
|
|
479 pdf(bcvPdf)
|
|
|
480 plotBCV(data, main="BCV Plot")
|
|
|
481 newEntry <- c("BCV Plot (.pdf)", "bcv.pdf")
|
|
|
482 invisible(dev.off())
|
|
|
483
|
|
0
|
484 if (workMode=="classic") {
|
|
|
485 # Assess differential representation using classic exact testing methodology
|
|
|
486 # in edgeR
|
|
|
487 testData <- exactTest(data, pair=pairData)
|
|
|
488
|
|
|
489 top <- topTags(testData, n=Inf)
|
|
|
490 topIDs <- top$table[(top$table$FDR < fdrThresh) &
|
|
|
491 (abs(top$table$logFC) > lfcThresh), 1]
|
|
3
|
492
|
|
0
|
493 write.table(top, file=topOut, row.names=FALSE, sep="\t")
|
|
3
|
494
|
|
0
|
495 linkName <- paste0("Top Tags Table(", pairData[2], "-", pairData[1],
|
|
|
496 ") (.tsv)")
|
|
|
497 linkAddr <- paste0("toptag(", pairData[2], "-", pairData[1], ").tsv")
|
|
|
498 linkData <- rbind(linkData, c(linkName, linkAddr))
|
|
|
499
|
|
3
|
500 upCount[1] <- sum(top$table$FDR < fdrThresh & top$table$logFC > lfcThresh)
|
|
|
501 downCount[1] <- sum(top$table$FDR < fdrThresh &
|
|
|
502 top$table$logFC < -lfcThresh)
|
|
|
503 flatCount[1] <- sum(top$table$FDR > fdrThresh |
|
|
|
504 abs(top$table$logFC) < lfcThresh)
|
|
|
505
|
|
|
506
|
|
|
507
|
|
0
|
508 # Select hairpins with FDR < 0.05 to highlight on plot
|
|
|
509 png(smearPng, width=600, height=600)
|
|
|
510 plotTitle <- gsub(".", " ",
|
|
|
511 paste0("Smear Plot: ", pairData[2], "-", pairData[1]),
|
|
3
|
512 fixed=TRUE)
|
|
0
|
513 plotSmear(testData, de.tags=topIDs,
|
|
|
514 pch=20, cex=1.0, main=plotTitle)
|
|
3
|
515 abline(h=c(-1, 0, 1), col=c("dodgerblue", "yellow", "dodgerblue"), lty=2)
|
|
0
|
516 imgName <- paste0("Smear Plot(", pairData[2], "-", pairData[1], ")")
|
|
|
517 imgAddr <- paste0("smear(", pairData[2], "-", pairData[1],").png")
|
|
|
518 imageData <- rbind(imageData, c(imgName, imgAddr))
|
|
|
519 invisible(dev.off())
|
|
|
520
|
|
|
521 pdf(smearPdf)
|
|
|
522 plotTitle <- gsub(".", " ",
|
|
|
523 paste0("Smear Plot: ", pairData[2], "-", pairData[1]),
|
|
3
|
524 fixed=TRUE)
|
|
0
|
525 plotSmear(testData, de.tags=topIDs,
|
|
|
526 pch=20, cex=1.0, main=plotTitle)
|
|
3
|
527 abline(h=c(-1, 0, 1), col=c("dodgerblue", "yellow", "dodgerblue"), lty=2)
|
|
0
|
528 imgName <- paste0("Smear Plot(", pairData[2], "-", pairData[1], ") (.pdf)")
|
|
|
529 imgAddr <- paste0("smear(", pairData[2], "-", pairData[1], ").pdf")
|
|
|
530 linkData <- rbind(linkData, c(imgName, imgAddr))
|
|
|
531 invisible(dev.off())
|
|
3
|
532
|
|
0
|
533 } else if (workMode=="glm") {
|
|
|
534 # Generating design information
|
|
|
535 factors <- factor(data$sample$group)
|
|
|
536 design <- model.matrix(~0+factors)
|
|
|
537
|
|
|
538 colnames(design) <- gsub("factors", "", colnames(design), fixed=TRUE)
|
|
|
539
|
|
|
540 # Split up contrasts seperated by comma into a vector
|
|
|
541 contrastData <- unlist(strsplit(contrastData, split=","))
|
|
3
|
542
|
|
0
|
543 for (i in 1:length(contrastData)) {
|
|
|
544 # Generate contrasts information
|
|
|
545 contrasts <- makeContrasts(contrasts=contrastData[i], levels=design)
|
|
|
546
|
|
|
547 # Fit negative bionomial GLM
|
|
3
|
548 fit <- glmFit(data, design)
|
|
0
|
549 # Carry out Likelihood ratio test
|
|
3
|
550 testData <- glmLRT(fit, contrast=contrasts)
|
|
0
|
551
|
|
|
552 # Select hairpins with FDR < 0.05 to highlight on plot
|
|
|
553 top <- topTags(testData, n=Inf)
|
|
|
554 topIDs <- top$table[(top$table$FDR < fdrThresh) &
|
|
|
555 (abs(top$table$logFC) > lfcThresh), 1]
|
|
|
556 write.table(top, file=topOut[i], row.names=FALSE, sep="\t")
|
|
|
557
|
|
|
558 linkName <- paste0("Top Tags Table(", contrastData[i], ") (.tsv)")
|
|
|
559 linkAddr <- paste0("toptag(", contrastData[i], ").tsv")
|
|
|
560 linkData <- rbind(linkData, c(linkName, linkAddr))
|
|
|
561
|
|
3
|
562 # Collect counts for differential representation
|
|
|
563 upCount[i] <- sum(top$table$FDR < fdrThresh & top$table$logFC > lfcThresh)
|
|
|
564 downCount[i] <- sum(top$table$FDR < fdrThresh &
|
|
|
565 top$table$logFC < -lfcThresh)
|
|
|
566 flatCount[i] <- sum(top$table$FDR > fdrThresh |
|
|
|
567 abs(top$table$logFC) < lfcThresh)
|
|
|
568
|
|
0
|
569 # Make a plot of logFC versus logCPM
|
|
|
570 png(smearPng[i], height=600, width=600)
|
|
|
571 plotTitle <- paste("Smear Plot:", gsub(".", " ", contrastData[i],
|
|
|
572 fixed=TRUE))
|
|
|
573 plotSmear(testData, de.tags=topIDs, pch=20, cex=0.8, main=plotTitle)
|
|
|
574 abline(h=c(-1, 0, 1), col=c("dodgerblue", "yellow", "dodgerblue"), lty=2)
|
|
|
575
|
|
|
576 imgName <- paste0("Smear Plot(", contrastData[i], ")")
|
|
|
577 imgAddr <- paste0("smear(", contrastData[i], ").png")
|
|
|
578 imageData <- rbind(imageData, c(imgName, imgAddr))
|
|
|
579 invisible(dev.off())
|
|
|
580
|
|
|
581 pdf(smearPdf[i])
|
|
|
582 plotTitle <- paste("Smear Plot:", gsub(".", " ", contrastData[i],
|
|
|
583 fixed=TRUE))
|
|
|
584 plotSmear(testData, de.tags=topIDs, pch=20, cex=0.8, main=plotTitle)
|
|
|
585 abline(h=c(-1, 0, 1), col=c("dodgerblue", "yellow", "dodgerblue"), lty=2)
|
|
|
586
|
|
|
587 linkName <- paste0("Smear Plot(", contrastData[i], ") (.pdf)")
|
|
|
588 linkAddr <- paste0("smear(", contrastData[i], ").pdf")
|
|
|
589 linkData <- rbind(linkData, c(linkName, linkAddr))
|
|
|
590 invisible(dev.off())
|
|
|
591
|
|
|
592 genes <- as.character(data$genes$Gene)
|
|
|
593 unq <- unique(genes)
|
|
|
594 unq <- unq[!is.na(unq)]
|
|
|
595 geneList <- list()
|
|
|
596 for (gene in unq) {
|
|
|
597 if (length(which(genes==gene)) >= hairpinReq) {
|
|
|
598 geneList[[gene]] <- which(genes==gene)
|
|
|
599 }
|
|
|
600 }
|
|
|
601
|
|
|
602 if (wantRoast) {
|
|
|
603 # Input preparaton for roast
|
|
3
|
604 nrot <- 9999
|
|
0
|
605 set.seed(602214129)
|
|
|
606 roastData <- mroast(data, index=geneList, design=design,
|
|
|
607 contrast=contrasts, nrot=nrot)
|
|
|
608 roastData <- cbind(GeneID=rownames(roastData), roastData)
|
|
|
609 write.table(roastData, file=roastOut[i], row.names=FALSE, sep="\t")
|
|
|
610 linkName <- paste0("Gene Level Analysis Table(", contrastData[i],
|
|
|
611 ") (.tsv)")
|
|
3
|
612 linkAddr <- paste0("gene_level(", contrastData[i], ").tsv")
|
|
0
|
613 linkData <- rbind(linkData, c(linkName, linkAddr))
|
|
|
614 if (selectOpt=="rank") {
|
|
|
615 selectedGenes <- rownames(roastData)[selectVals]
|
|
|
616 } else {
|
|
|
617 selectedGenes <- selectVals
|
|
|
618 }
|
|
|
619
|
|
|
620 if (packageVersion("limma")<"3.19.19") {
|
|
|
621 png(barcodePng[i], width=600, height=length(selectedGenes)*150)
|
|
|
622 } else {
|
|
|
623 png(barcodePng[i], width=600, height=length(selectedGenes)*300)
|
|
|
624 }
|
|
|
625 par(mfrow=c(length(selectedGenes), 1))
|
|
|
626 for (gene in selectedGenes) {
|
|
|
627 barcodeplot(testData$table$logFC, index=geneList[[gene]],
|
|
|
628 main=paste("Barcode Plot for", gene, "(logFCs)",
|
|
|
629 gsub(".", " ", contrastData[i])),
|
|
|
630 labels=c("Positive logFC", "Negative logFC"))
|
|
|
631 }
|
|
|
632 imgName <- paste0("Barcode Plot(", contrastData[i], ")")
|
|
|
633 imgAddr <- paste0("barcode(", contrastData[i], ").png")
|
|
|
634 imageData <- rbind(imageData, c(imgName, imgAddr))
|
|
|
635 dev.off()
|
|
|
636 if (packageVersion("limma")<"3.19.19") {
|
|
|
637 pdf(barcodePdf[i], width=8, height=2)
|
|
|
638 } else {
|
|
|
639 pdf(barcodePdf[i], width=8, height=4)
|
|
|
640 }
|
|
|
641 for (gene in selectedGenes) {
|
|
|
642 barcodeplot(testData$table$logFC, index=geneList[[gene]],
|
|
|
643 main=paste("Barcode Plot for", gene, "(logFCs)",
|
|
|
644 gsub(".", " ", contrastData[i])),
|
|
|
645 labels=c("Positive logFC", "Negative logFC"))
|
|
|
646 }
|
|
|
647 linkName <- paste0("Barcode Plot(", contrastData[i], ") (.pdf)")
|
|
|
648 linkAddr <- paste0("barcode(", contrastData[i], ").pdf")
|
|
|
649 linkData <- rbind(linkData, c(linkName, linkAddr))
|
|
|
650 dev.off()
|
|
|
651 }
|
|
|
652 }
|
|
|
653 }
|
|
|
654
|
|
3
|
655 # Generate data frame of the significant differences
|
|
|
656 sigDiff <- data.frame(Up=upCount, Flat=flatCount, Down=downCount)
|
|
|
657 if (workMode == "glm") {
|
|
|
658 row.names(sigDiff) <- contrastData
|
|
|
659 } else if (workMode == "classic") {
|
|
|
660 row.names(sigDiff) <- paste0(pairData[2], "-", pairData[1])
|
|
|
661 }
|
|
|
662
|
|
|
663 # Output table of summarised counts
|
|
|
664 ID <- rownames(data$counts)
|
|
|
665 outputCounts <- cbind(ID, data$counts)
|
|
|
666 write.table(outputCounts, file=countsOut, row.names=FALSE, sep="\t",
|
|
|
667 quote=FALSE)
|
|
|
668 linkName <- "Counts table (.tsv)"
|
|
|
669 linkAddr <- "counts.tsv"
|
|
|
670 linkData <- rbind(linkData, c(linkName, linkAddr))
|
|
|
671
|
|
|
672 # Record session info
|
|
|
673 writeLines(capture.output(sessionInfo()), sessionOut)
|
|
|
674 linkData <- rbind(linkData, c("Session Info", "session_info.txt"))
|
|
|
675
|
|
1
|
676 # Record ending time and calculate total run time
|
|
0
|
677 timeEnd <- as.character(Sys.time())
|
|
1
|
678 timeTaken <- capture.output(round(difftime(timeEnd,timeStart), digits=3))
|
|
|
679 timeTaken <- gsub("Time difference of ", "", timeTaken, fixed=TRUE)
|
|
0
|
680 ################################################################################
|
|
|
681 ### HTML Generation
|
|
|
682 ################################################################################
|
|
|
683 # Clear file
|
|
|
684 cat("", file=htmlPath)
|
|
|
685
|
|
|
686 cata("<html>\n")
|
|
|
687 HtmlHead("EdgeR Output")
|
|
|
688
|
|
|
689 cata("<body>\n")
|
|
|
690 cata("<h3>EdgeR Analysis Output:</h3>\n")
|
|
|
691 cata("<h4>Input Summary:</h4>\n")
|
|
|
692 if (inputType=="fastq") {
|
|
|
693 cata("<ul>\n")
|
|
|
694 ListItem(hpReadout[1])
|
|
|
695 ListItem(hpReadout[2])
|
|
|
696 cata("</ul>\n")
|
|
3
|
697 cata(hpReadout[3], "<br />\n")
|
|
0
|
698 cata("<ul>\n")
|
|
|
699 ListItem(hpReadout[4])
|
|
|
700 ListItem(hpReadout[7])
|
|
|
701 cata("</ul>\n")
|
|
3
|
702 cata(hpReadout[8:11], sep="<br />\n")
|
|
0
|
703 cata("<br />\n")
|
|
|
704 cata("<b>Please check that read percentages are consistent with ")
|
|
|
705 cata("expectations.</b><br >\n")
|
|
|
706 } else if (inputType=="counts") {
|
|
|
707 cata("<ul>\n")
|
|
|
708 ListItem("Number of Samples: ", ncol(data$counts))
|
|
|
709 ListItem("Number of Hairpins: ", countsRows)
|
|
|
710 ListItem("Number of annotations provided: ", annoRows)
|
|
|
711 ListItem("Number of annotations matched to hairpin: ", annoMatched)
|
|
|
712 cata("</ul>\n")
|
|
|
713 }
|
|
|
714
|
|
|
715 cata("The estimated common biological coefficient of variation (BCV) is: ",
|
|
|
716 commonBCV, "<br />\n")
|
|
|
717
|
|
|
718 cata("<h4>Output:</h4>\n")
|
|
3
|
719 cata("PDF copies of JPEGS available in 'Plots' section.<br />\n")
|
|
0
|
720 for (i in 1:nrow(imageData)) {
|
|
|
721 if (grepl("barcode", imageData$Link[i])) {
|
|
|
722 if (packageVersion("limma")<"3.19.19") {
|
|
|
723 HtmlImage(imageData$Link[i], imageData$Label[i],
|
|
|
724 height=length(selectedGenes)*150)
|
|
|
725 } else {
|
|
|
726 HtmlImage(imageData$Link[i], imageData$Label[i],
|
|
|
727 height=length(selectedGenes)*300)
|
|
|
728 }
|
|
|
729 } else {
|
|
|
730 HtmlImage(imageData$Link[i], imageData$Label[i])
|
|
|
731 }
|
|
|
732 }
|
|
3
|
733 cata("<br />\n")
|
|
|
734
|
|
|
735 cata("<h4>Differential Representation Counts:</h4>\n")
|
|
|
736
|
|
|
737 cata("<table border=\"1\" cellpadding=\"4\">\n")
|
|
|
738 cata("<tr>\n")
|
|
|
739 TableItem()
|
|
|
740 for (i in colnames(sigDiff)) {
|
|
|
741 TableHeadItem(i)
|
|
|
742 }
|
|
|
743 cata("</tr>\n")
|
|
|
744 for (i in 1:nrow(sigDiff)) {
|
|
|
745 cata("<tr>\n")
|
|
|
746 TableHeadItem(unmake.names(row.names(sigDiff)[i]))
|
|
|
747 for (j in 1:ncol(sigDiff)) {
|
|
|
748 TableItem(as.character(sigDiff[i, j]))
|
|
|
749 }
|
|
|
750 cata("</tr>\n")
|
|
|
751 }
|
|
|
752 cata("</table>")
|
|
0
|
753
|
|
|
754 cata("<h4>Plots:</h4>\n")
|
|
|
755 for (i in 1:nrow(linkData)) {
|
|
3
|
756 if (grepl(".pdf", linkData$Link[i])) {
|
|
0
|
757 HtmlLink(linkData$Link[i], linkData$Label[i])
|
|
|
758 }
|
|
|
759 }
|
|
|
760
|
|
|
761 cata("<h4>Tables:</h4>\n")
|
|
|
762 for (i in 1:nrow(linkData)) {
|
|
|
763 if (grepl(".tsv", linkData$Link[i])) {
|
|
|
764 HtmlLink(linkData$Link[i], linkData$Label[i])
|
|
|
765 }
|
|
|
766 }
|
|
|
767
|
|
3
|
768 cata("<p>Alt-click links to download file.</p>\n")
|
|
|
769 cata("<p>Click floppy disc icon associated history item to download ")
|
|
|
770 cata("all files.</p>\n")
|
|
|
771 cata("<p>.tsv files can be viewed in Excel or any spreadsheet program.</p>\n")
|
|
0
|
772
|
|
1
|
773 cata("<h4>Additional Information:</h4>\n")
|
|
|
774
|
|
|
775 if (inputType == "fastq") {
|
|
|
776 ListItem("Data was gathered from fastq raw read file(s).")
|
|
|
777 } else if (inputType == "counts") {
|
|
|
778 ListItem("Data was gathered from a table of counts.")
|
|
|
779 }
|
|
|
780
|
|
|
781 if (cpmReq!=0 && sampleReq!=0) {
|
|
3
|
782 tempStr <- paste("Hairpins without more than", cpmReq,
|
|
|
783 "CPM in at least", sampleReq, "samples are insignificant",
|
|
|
784 "and filtered out.")
|
|
1
|
785 ListItem(tempStr)
|
|
|
786 filterProp <- round(filteredCount/preFilterCount*100, digits=2)
|
|
|
787 tempStr <- paste0(filteredCount, " of ", preFilterCount," (", filterProp,
|
|
|
788 "%) hairpins were filtered out for low count-per-million.")
|
|
|
789 ListItem(tempStr)
|
|
|
790 }
|
|
|
791
|
|
|
792 if (workMode == "classic") {
|
|
|
793 ListItem("An exact test was performed on each hairpin.")
|
|
|
794 } else if (workMode == "glm") {
|
|
|
795 ListItem("A generalised linear model was fitted to each hairpin.")
|
|
|
796 }
|
|
|
797
|
|
|
798 cit <- character()
|
|
|
799 link <-character()
|
|
|
800 link[1] <- paste0("<a href=\"",
|
|
|
801 "http://www.bioconductor.org/packages/release/bioc/",
|
|
|
802 "vignettes/limma/inst/doc/usersguide.pdf",
|
|
|
803 "\">", "limma User's Guide", "</a>.")
|
|
|
804 link[2] <- paste0("<a href=\"",
|
|
|
805 "http://www.bioconductor.org/packages/release/bioc/",
|
|
|
806 "vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf",
|
|
|
807 "\">", "edgeR User's Guide", "</a>")
|
|
|
808
|
|
|
809 cit[1] <- paste("Robinson MD, McCarthy DJ and Smyth GK (2010).",
|
|
|
810 "edgeR: a Bioconductor package for differential",
|
|
|
811 "expression analysis of digital gene expression",
|
|
|
812 "data. Bioinformatics 26, 139-140")
|
|
|
813 cit[2] <- paste("Robinson MD and Smyth GK (2007). Moderated statistical tests",
|
|
|
814 "for assessing differences in tag abundance. Bioinformatics",
|
|
|
815 "23, 2881-2887")
|
|
|
816 cit[3] <- paste("Robinson MD and Smyth GK (2008). Small-sample estimation of",
|
|
|
817 "negative binomial dispersion, with applications to SAGE data.",
|
|
|
818 "Biostatistics, 9, 321-332")
|
|
|
819
|
|
|
820 cit[4] <- paste("McCarthy DJ, Chen Y and Smyth GK (2012). Differential",
|
|
|
821 "expression analysis of multifactor RNA-Seq experiments with",
|
|
|
822 "respect to biological variation. Nucleic Acids Research 40,",
|
|
|
823 "4288-4297")
|
|
|
824
|
|
|
825 cata("<h4>Citations</h4>")
|
|
|
826 cata("<ol>\n")
|
|
|
827 ListItem(cit[1])
|
|
|
828 ListItem(cit[2])
|
|
|
829 ListItem(cit[3])
|
|
|
830 ListItem(cit[4])
|
|
|
831 cata("</ol>\n")
|
|
|
832
|
|
3
|
833 cata("<p>Report problems to: su.s@wehi.edu.au</p>\n")
|
|
|
834
|
|
|
835 for (i in 1:nrow(linkData)) {
|
|
|
836 if (grepl("session_info", linkData$Link[i])) {
|
|
|
837 HtmlLink(linkData$Link[i], linkData$Label[i])
|
|
|
838 }
|
|
|
839 }
|
|
|
840
|
|
1
|
841 cata("<table border=\"0\">\n")
|
|
0
|
842 cata("<tr>\n")
|
|
|
843 TableItem("Task started at:"); TableItem(timeStart)
|
|
|
844 cata("</tr>\n")
|
|
|
845 cata("<tr>\n")
|
|
|
846 TableItem("Task ended at:"); TableItem(timeEnd)
|
|
|
847 cata("</tr>\n")
|
|
1
|
848 cata("<tr>\n")
|
|
|
849 TableItem("Task run time:"); TableItem(timeTaken)
|
|
|
850 cata("<tr>\n")
|
|
|
851 cata("</table>\n")
|
|
0
|
852
|
|
|
853 cata("</body>\n")
|
|
|
854 cata("</html>")
|
