|
0
|
1 # This tool takes in a matrix of feature counts as well as gene annotations and
|
|
|
2 # outputs a table of top expressions as well as various plots for differential
|
|
|
3 # expression analysis
|
|
|
4 #
|
|
|
5 # ARGS: 1.countPath -Path to RData input containing counts
|
|
|
6 # 2.annoPath -Path to RData input containing gene annotations
|
|
|
7 # 3.htmlPath -Path to html file linking to other outputs
|
|
|
8 # 4.outPath -Path to folder to write all output to
|
|
|
9 # 5.rdaOpt -String specifying if RData should be saved
|
|
|
10 # 6.normOpt -String specifying type of normalisation used
|
|
|
11 # 7.weightOpt -String specifying usage of weights
|
|
|
12 # 8.contrastData -String containing contrasts of interest
|
|
|
13 # 9.cpmReq -Float specifying cpm requirement
|
|
|
14 # 10.sampleReq -Integer specifying cpm requirement
|
|
|
15 # 11.pAdjOpt -String specifying the p-value adjustment method
|
|
|
16 # 12.pValReq -Float specifying the p-value requirement
|
|
|
17 # 13.lfcReq -Float specifying the log-fold-change requirement
|
|
|
18 # 14.factorData -String containing factor names and values
|
|
|
19 #
|
|
|
20 # OUT: Voom Plot
|
|
|
21 # BCV Plot
|
|
|
22 # MA Plot
|
|
|
23 # Top Expression Table
|
|
|
24 # HTML file linking to the ouputs
|
|
|
25 #
|
|
|
26 # Author: Shian Su - registertonysu@gmail.com - Jan 2014
|
|
|
27
|
|
|
28 # Record starting time
|
|
|
29 timeStart <- as.character(Sys.time())
|
|
|
30
|
|
|
31 # Load all required libraries
|
|
|
32 library(methods, quietly=TRUE, warn.conflicts=FALSE)
|
|
|
33 library(statmod, quietly=TRUE, warn.conflicts=FALSE)
|
|
|
34 library(splines, quietly=TRUE, warn.conflicts=FALSE)
|
|
|
35 library(edgeR, quietly=TRUE, warn.conflicts=FALSE)
|
|
|
36 library(limma, quietly=TRUE, warn.conflicts=FALSE)
|
|
|
37 library(scales, quietly=TRUE, warn.conflicts=FALSE)
|
|
|
38
|
|
|
39 if (packageVersion("limma") < "3.20.1") {
|
|
|
40 stop("Please update 'limma' to version >= 3.20.1 to run this tool")
|
|
|
41 }
|
|
|
42
|
|
|
43 ################################################################################
|
|
|
44 ### Function Delcaration
|
|
|
45 ################################################################################
|
|
|
46 # Function to sanitise contrast equations so there are no whitespaces
|
|
|
47 # surrounding the arithmetic operators, leading or trailing whitespace
|
|
|
48 sanitiseEquation <- function(equation) {
|
|
|
49 equation <- gsub(" *[+] *", "+", equation)
|
|
|
50 equation <- gsub(" *[-] *", "-", equation)
|
|
|
51 equation <- gsub(" *[/] *", "/", equation)
|
|
|
52 equation <- gsub(" *[*] *", "*", equation)
|
|
|
53 equation <- gsub("^\\s+|\\s+$", "", equation)
|
|
|
54 return(equation)
|
|
|
55 }
|
|
|
56
|
|
|
57 # Function to sanitise group information
|
|
|
58 sanitiseGroups <- function(string) {
|
|
|
59 string <- gsub(" *[,] *", ",", string)
|
|
|
60 string <- gsub("^\\s+|\\s+$", "", string)
|
|
|
61 return(string)
|
|
|
62 }
|
|
|
63
|
|
|
64 # Function to change periods to whitespace in a string
|
|
|
65 unmake.names <- function(string) {
|
|
|
66 string <- gsub(".", " ", string, fixed=TRUE)
|
|
|
67 return(string)
|
|
|
68 }
|
|
|
69
|
|
|
70 # Generate output folder and paths
|
|
|
71 makeOut <- function(filename) {
|
|
|
72 return(paste0(outPath, "/", filename))
|
|
|
73 }
|
|
|
74
|
|
|
75 # Generating design information
|
|
|
76 pasteListName <- function(string) {
|
|
|
77 return(paste0("factors$", string))
|
|
|
78 }
|
|
|
79
|
|
|
80 # Create cata function: default path set, default seperator empty and appending
|
|
|
81 # true by default (Ripped straight from the cat function with altered argument
|
|
|
82 # defaults)
|
|
|
83 cata <- function(..., file = htmlPath, sep = "", fill = FALSE, labels = NULL,
|
|
|
84 append = TRUE) {
|
|
|
85 if (is.character(file))
|
|
|
86 if (file == "")
|
|
|
87 file <- stdout()
|
|
|
88 else if (substring(file, 1L, 1L) == "|") {
|
|
|
89 file <- pipe(substring(file, 2L), "w")
|
|
|
90 on.exit(close(file))
|
|
|
91 }
|
|
|
92 else {
|
|
|
93 file <- file(file, ifelse(append, "a", "w"))
|
|
|
94 on.exit(close(file))
|
|
|
95 }
|
|
|
96 .Internal(cat(list(...), file, sep, fill, labels, append))
|
|
|
97 }
|
|
|
98
|
|
|
99 # Function to write code for html head and title
|
|
|
100 HtmlHead <- function(title) {
|
|
|
101 cata("<head>\n")
|
|
|
102 cata("<title>", title, "</title>\n")
|
|
|
103 cata("</head>\n")
|
|
|
104 }
|
|
|
105
|
|
|
106 # Function to write code for html links
|
|
|
107 HtmlLink <- function(address, label=address) {
|
|
|
108 cata("<a href=\"", address, "\" target=\"_blank\">", label, "</a><br />\n")
|
|
|
109 }
|
|
|
110
|
|
|
111 # Function to write code for html images
|
|
|
112 HtmlImage <- function(source, label=source, height=600, width=600) {
|
|
|
113 cata("<img src=\"", source, "\" alt=\"", label, "\" height=\"", height)
|
|
|
114 cata("\" width=\"", width, "\"/>\n")
|
|
|
115 }
|
|
|
116
|
|
|
117 # Function to write code for html list items
|
|
|
118 ListItem <- function(...) {
|
|
|
119 cata("<li>", ..., "</li>\n")
|
|
|
120 }
|
|
|
121
|
|
|
122 TableItem <- function(...) {
|
|
|
123 cata("<td>", ..., "</td>\n")
|
|
|
124 }
|
|
|
125
|
|
|
126 TableHeadItem <- function(...) {
|
|
|
127 cata("<th>", ..., "</th>\n")
|
|
|
128 }
|
|
|
129
|
|
|
130 ################################################################################
|
|
|
131 ### Input Processing
|
|
|
132 ################################################################################
|
|
|
133
|
|
|
134 # Collects arguments from command line
|
|
|
135 argv <- commandArgs(TRUE)
|
|
|
136
|
|
|
137 # Grab arguments
|
|
|
138 countPath <- as.character(argv[1])
|
|
|
139 annoPath <- as.character(argv[2])
|
|
|
140 htmlPath <- as.character(argv[3])
|
|
|
141 outPath <- as.character(argv[4])
|
|
|
142 rdaOpt <- as.character(argv[5])
|
|
|
143 normOpt <- as.character(argv[6])
|
|
|
144 weightOpt <- as.character(argv[7])
|
|
|
145 contrastData <- as.character(argv[8])
|
|
|
146 cpmReq <- as.numeric(argv[9])
|
|
|
147 sampleReq <- as.numeric(argv[10])
|
|
|
148 pAdjOpt <- as.character(argv[11])
|
|
|
149 pValReq <- as.numeric(argv[12])
|
|
|
150 lfcReq <- as.numeric(argv[13])
|
|
|
151 factorData <- list()
|
|
|
152 for (i in 14:length(argv)) {
|
|
|
153 newFact <- unlist(strsplit(as.character(argv[i]), split="::"))
|
|
|
154 factorData <- rbind(factorData, newFact)
|
|
|
155 } # Factors have the form: FACT_NAME::LEVEL,LEVEL,LEVEL,LEVEL,...
|
|
|
156
|
|
|
157 # Process arguments
|
|
|
158 if (weightOpt=="yes") {
|
|
|
159 wantWeight <- TRUE
|
|
|
160 } else {
|
|
|
161 wantWeight <- FALSE
|
|
|
162 }
|
|
|
163
|
|
|
164 if (rdaOpt=="yes") {
|
|
|
165 wantRda <- TRUE
|
|
|
166 } else {
|
|
|
167 wantRda <- FALSE
|
|
|
168 }
|
|
|
169
|
|
|
170 if (annoPath=="None") {
|
|
|
171 haveAnno <- FALSE
|
|
|
172 } else {
|
|
|
173 haveAnno <- TRUE
|
|
|
174 }
|
|
|
175
|
|
|
176 # Set the row names to be the name of the factor and delete first row
|
|
|
177 row.names(factorData) <- factorData[, 1]
|
|
|
178 factorData <- factorData[, -1]
|
|
|
179 factorData <- sapply(factorData, sanitiseGroups)
|
|
|
180 factorData <- sapply(factorData, strsplit, split=",")
|
|
|
181 factorData <- sapply(factorData, make.names)
|
|
|
182
|
|
|
183 # Transform factor data into data frame of R factor objects
|
|
|
184 factors <- data.frame(factorData)
|
|
|
185
|
|
|
186 #Create output directory
|
|
|
187 dir.create(outPath, showWarnings=FALSE)
|
|
|
188
|
|
|
189 # Split up contrasts seperated by comma into a vector then sanitise
|
|
|
190 contrastData <- unlist(strsplit(contrastData, split=","))
|
|
|
191 contrastData <- sanitiseEquation(contrastData)
|
|
|
192 contrastData <- gsub(" ", ".", contrastData, fixed=TRUE)
|
|
|
193
|
|
|
194 bcvOutPdf <- makeOut("bcvplot.pdf")
|
|
|
195 bcvOutPng <- makeOut("bcvplot.png")
|
|
|
196 mdsOutPdf <- makeOut("mdsplot.pdf")
|
|
|
197 mdsOutPng <- makeOut("mdsplot.png")
|
|
|
198 voomOutPdf <- makeOut("voomplot.pdf")
|
|
|
199 voomOutPng <- makeOut("voomplot.png")
|
|
|
200 maOutPdf <- character() # Initialise character vector
|
|
|
201 maOutPng <- character()
|
|
|
202 topOut <- character()
|
|
|
203 for (i in 1:length(contrastData)) {
|
|
|
204 maOutPdf[i] <- makeOut(paste0("maplot(", contrastData[i], ").pdf"))
|
|
|
205 maOutPng[i] <- makeOut(paste0("maplot(", contrastData[i], ").png"))
|
|
|
206 topOut[i] <- makeOut(paste0("toptab(", contrastData[i], ").tsv"))
|
|
|
207 } # Save output paths for each contrast as vectors
|
|
|
208 rdaOut <- makeOut("RData.rda")
|
|
|
209 sessionOut <- makeOut("session_info.txt")
|
|
|
210
|
|
|
211 # Initialise data for html links and images, data frame with columns Label and
|
|
|
212 # Link
|
|
|
213 linkData <- data.frame(Label=character(), Link=character(),
|
|
|
214 stringsAsFactors=FALSE)
|
|
|
215 imageData <- data.frame(Label=character(), Link=character(),
|
|
|
216 stringsAsFactors=FALSE)
|
|
|
217
|
|
|
218 # Initialise vectors for storage of up/down/neutral regulated counts
|
|
|
219 upCount <- numeric()
|
|
|
220 downCount <- numeric()
|
|
|
221 flatCount <- numeric()
|
|
|
222
|
|
|
223 # Read in counts and geneanno data
|
|
|
224 counts <- read.table(countPath, header=TRUE, sep="\t")
|
|
|
225 row.names(counts) <- counts$GeneID
|
|
|
226 counts <- counts[ , !(colnames(counts)=="GeneID")]
|
|
|
227 countsRows <- nrow(counts)
|
|
|
228 if (haveAnno) {
|
|
|
229 geneanno <- read.table(annoPath, header=TRUE, sep="\t")
|
|
|
230 }
|
|
|
231
|
|
|
232 ################################################################################
|
|
|
233 ### Data Processing
|
|
|
234 ################################################################################
|
|
|
235
|
|
|
236 # Extract counts and annotation data
|
|
|
237 data <- list()
|
|
|
238 data$counts <- counts
|
|
|
239 if (haveAnno) {
|
|
|
240 data$genes <- geneanno
|
|
|
241 } else {
|
|
|
242 data$genes <- data.frame(GeneID=row.names(counts))
|
|
|
243 }
|
|
|
244
|
|
|
245 # Filter out genes that do not have a required cpm in a required number of
|
|
|
246 # samples
|
|
|
247 preFilterCount <- nrow(data$counts)
|
|
|
248 sel <- rowSums(cpm(data$counts) > cpmReq) >= sampleReq
|
|
|
249 data$counts <- data$counts[sel, ]
|
|
|
250 data$genes <- data$genes[sel, ]
|
|
|
251 postFilterCount <- nrow(data$counts)
|
|
|
252 filteredCount <- preFilterCount-postFilterCount
|
|
|
253
|
|
|
254 # Creating naming data
|
|
|
255 samplenames <- colnames(data$counts)
|
|
|
256 sampleanno <- data.frame("sampleID"=samplenames, factors)
|
|
|
257
|
|
|
258 # Generating the DGEList object "data"
|
|
|
259 data$samples <- sampleanno
|
|
|
260 data$samples$lib.size <- colSums(data$counts)
|
|
|
261 data$samples$norm.factors <- 1
|
|
|
262 row.names(data$samples) <- colnames(data$counts)
|
|
|
263 data <- new("DGEList", data)
|
|
|
264
|
|
|
265 factorList <- sapply(names(factors), pasteListName)
|
|
|
266 formula <- "~0"
|
|
|
267 for (i in 1:length(factorList)) {
|
|
|
268 formula <- paste(formula, factorList[i], sep="+")
|
|
|
269 }
|
|
|
270 formula <- formula(formula)
|
|
|
271 design <- model.matrix(formula)
|
|
|
272 for (i in 1:length(factorList)) {
|
|
|
273 colnames(design) <- gsub(factorList[i], "", colnames(design), fixed=TRUE)
|
|
|
274 }
|
|
|
275
|
|
|
276 # Calculating normalising factor, estimating dispersion
|
|
|
277 data <- calcNormFactors(data, method=normOpt)
|
|
|
278 #data <- estimateDisp(data, design=design, robust=TRUE)
|
|
|
279
|
|
|
280 # Generate contrasts information
|
|
|
281 contrasts <- makeContrasts(contrasts=contrastData, levels=design)
|
|
|
282
|
|
|
283 # Name rows of factors according to their sample
|
|
|
284 row.names(factors) <- names(data$counts)
|
|
|
285
|
|
|
286 ################################################################################
|
|
|
287 ### Data Output
|
|
|
288 ################################################################################
|
|
|
289
|
|
|
290 # BCV Plot
|
|
|
291 #png(bcvOutPng, width=600, height=600)
|
|
|
292 #plotBCV(data, main="BCV Plot")
|
|
|
293 #imageData[1, ] <- c("BCV Plot", "bcvplot.png")
|
|
|
294 #invisible(dev.off())
|
|
|
295
|
|
|
296 #pdf(bcvOutPdf)
|
|
|
297 #plotBCV(data, main="BCV Plot")
|
|
|
298 #invisible(dev.off())
|
|
|
299
|
|
|
300 if (wantWeight) {
|
|
|
301 # Creating voom data object and plot
|
|
|
302 png(voomOutPng, width=1000, height=600)
|
|
|
303 vData <- voomWithQualityWeights(data, design=design, plot=TRUE)
|
|
|
304 imageData[1, ] <- c("Voom Plot", "voomplot.png")
|
|
|
305 invisible(dev.off())
|
|
|
306
|
|
|
307 pdf(voomOutPdf, width=14)
|
|
|
308 vData <- voomWithQualityWeights(data, design=design, plot=TRUE)
|
|
|
309 linkData[1, ] <- c("Voom Plot (.pdf)", "voomplot.pdf")
|
|
|
310 invisible(dev.off())
|
|
|
311
|
|
|
312 # Generating fit data and top table with weights
|
|
|
313 wts <- vData$weights
|
|
|
314 voomFit <- lmFit(vData, design, weights=wts)
|
|
|
315
|
|
|
316 } else {
|
|
|
317 # Creating voom data object and plot
|
|
|
318 png(voomOutPng, width=600, height=600)
|
|
|
319 vData <- voom(data, design=design, plot=TRUE)
|
|
|
320 imageData[1, ] <- c("Voom Plot", "voomplot.png")
|
|
|
321 invisible(dev.off())
|
|
|
322
|
|
|
323 pdf(voomOutPdf)
|
|
|
324 vData <- voom(data, design=design, plot=TRUE)
|
|
|
325 linkData[1, ] <- c("Voom Plot (.pdf)", "voomplot.pdf")
|
|
|
326 invisible(dev.off())
|
|
|
327
|
|
|
328 # Generate voom fit
|
|
|
329 voomFit <- lmFit(vData, design)
|
|
|
330
|
|
|
331 }
|
|
|
332
|
|
|
333 # Fit linear model and estimate dispersion with eBayes
|
|
|
334 voomFit <- contrasts.fit(voomFit, contrasts)
|
|
|
335 voomFit <- eBayes(voomFit)
|
|
|
336
|
|
|
337 # Plot MDS
|
|
|
338 labels <- names(counts)
|
|
|
339 png(mdsOutPng, width=600, height=600)
|
|
|
340 # Currently only using a single factor
|
|
|
341 plotMDS(vData, labels=labels, col=as.numeric(factors[, 1]), cex=0.8)
|
|
|
342 imgName <- "Voom Plot"
|
|
|
343 imgAddr <- "mdsplot.png"
|
|
|
344 imageData <- rbind(imageData, c(imgName, imgAddr))
|
|
|
345 invisible(dev.off())
|
|
|
346
|
|
|
347 pdf(mdsOutPdf)
|
|
|
348 plotMDS(vData, labels=labels, cex=0.5)
|
|
|
349 linkName <- paste0("MDS Plot (.pdf)")
|
|
|
350 linkAddr <- paste0("mdsplot.pdf")
|
|
|
351 linkData <- rbind(linkData, c(linkName, linkAddr))
|
|
|
352 invisible(dev.off())
|
|
|
353
|
|
|
354
|
|
|
355 for (i in 1:length(contrastData)) {
|
|
|
356
|
|
|
357 status = decideTests(voomFit[, i], adjust.method=pAdjOpt, p.value=pValReq,
|
|
|
358 lfc=lfcReq)
|
|
|
359
|
|
|
360 sumStatus <- summary(status)
|
|
|
361
|
|
|
362 # Collect counts for differential expression
|
|
|
363 upCount[i] <- sumStatus["1",]
|
|
|
364 downCount[i] <- sumStatus["-1",]
|
|
|
365 flatCount[i] <- sumStatus["0",]
|
|
|
366
|
|
|
367 # Write top expressions table
|
|
|
368 top <- topTable(voomFit, coef=i, number=Inf, sort.by="P")
|
|
|
369 write.table(top, file=topOut[i], row.names=FALSE, sep="\t")
|
|
|
370
|
|
|
371 linkName <- paste0("Top Differential Expressions(", contrastData[i],
|
|
|
372 ") (.tsv)")
|
|
|
373 linkAddr <- paste0("toptab(", contrastData[i], ").tsv")
|
|
|
374 linkData <- rbind(linkData, c(linkName, linkAddr))
|
|
|
375
|
|
|
376 # Plot MA (log ratios vs mean average) using limma package on weighted data
|
|
|
377 pdf(maOutPdf[i])
|
|
|
378 limma::plotMA(voomFit, status=status, coef=i,
|
|
|
379 main=paste("MA Plot:", unmake.names(contrastData[i])),
|
|
|
380 col=alpha(c("firebrick", "blue"), 0.4), values=c("1", "-1"),
|
|
|
381 xlab="Average Expression", ylab="logFC")
|
|
|
382
|
|
|
383 abline(h=0, col="grey", lty=2)
|
|
|
384
|
|
|
385 linkName <- paste0("MA Plot(", contrastData[i], ")", " (.pdf)")
|
|
|
386 linkAddr <- paste0("maplot(", contrastData[i], ").pdf")
|
|
|
387 linkData <- rbind(linkData, c(linkName, linkAddr))
|
|
|
388 invisible(dev.off())
|
|
|
389
|
|
|
390 png(maOutPng[i], height=600, width=600)
|
|
|
391 limma::plotMA(voomFit, status=status, coef=i,
|
|
|
392 main=paste("MA Plot:", unmake.names(contrastData[i])),
|
|
|
393 col=alpha(c("firebrick", "blue"), 0.4), values=c("1", "-1"),
|
|
|
394 xlab="Average Expression", ylab="logFC")
|
|
|
395
|
|
|
396 abline(h=0, col="grey", lty=2)
|
|
|
397
|
|
|
398 imgName <- paste0("MA Plot(", contrastData[i], ")")
|
|
|
399 imgAddr <- paste0("maplot(", contrastData[i], ").png")
|
|
|
400 imageData <- rbind(imageData, c(imgName, imgAddr))
|
|
|
401 invisible(dev.off())
|
|
|
402 }
|
|
|
403 sigDiff <- data.frame(Up=upCount, Flat=flatCount, Down=downCount)
|
|
|
404 row.names(sigDiff) <- contrastData
|
|
|
405
|
|
|
406 # Save relevant items as rda object
|
|
|
407 if (wantRda) {
|
|
|
408 if (wantWeight) {
|
|
|
409 save(data, status, vData, labels, factors, wts, voomFit, top, contrasts,
|
|
|
410 design,
|
|
|
411 file=rdaOut, ascii=TRUE)
|
|
|
412 } else {
|
|
|
413 save(data, status, vData, labels, factors, voomFit, top, contrasts, design,
|
|
|
414 file=rdaOut, ascii=TRUE)
|
|
|
415 }
|
|
|
416 linkData <- rbind(linkData, c("RData (.rda)", "RData.rda"))
|
|
|
417 }
|
|
|
418
|
|
|
419 # Record session info
|
|
|
420 writeLines(capture.output(sessionInfo()), sessionOut)
|
|
|
421 linkData <- rbind(linkData, c("Session Info", "session_info.txt"))
|
|
|
422
|
|
|
423 # Record ending time and calculate total run time
|
|
|
424 timeEnd <- as.character(Sys.time())
|
|
|
425 timeTaken <- capture.output(round(difftime(timeEnd,timeStart), digits=3))
|
|
|
426 timeTaken <- gsub("Time difference of ", "", timeTaken, fixed=TRUE)
|
|
|
427 ################################################################################
|
|
|
428 ### HTML Generation
|
|
|
429 ################################################################################
|
|
|
430
|
|
|
431 # Clear file
|
|
|
432 cat("", file=htmlPath)
|
|
|
433
|
|
|
434 cata("<html>\n")
|
|
|
435
|
|
|
436 cata("<body>\n")
|
|
|
437 cata("<h3>Limma Voom Analysis Output:</h3>\n")
|
|
|
438 cata("PDF copies of JPEGS available in 'Plots' section.<br />\n")
|
|
|
439 if (wantWeight) {
|
|
|
440 HtmlImage(imageData$Link[1], imageData$Label[1], width=1000)
|
|
|
441 } else {
|
|
|
442 HtmlImage(imageData$Link[1], imageData$Label[1])
|
|
|
443 }
|
|
|
444
|
|
|
445 for (i in 2:nrow(imageData)) {
|
|
|
446 HtmlImage(imageData$Link[i], imageData$Label[i])
|
|
|
447 }
|
|
|
448
|
|
|
449 cata("<h4>Differential Expression Counts:</h4>\n")
|
|
|
450
|
|
|
451 cata("<table border=\"1\" cellpadding=\"4\">\n")
|
|
|
452 cata("<tr>\n")
|
|
|
453 TableItem()
|
|
|
454 for (i in colnames(sigDiff)) {
|
|
|
455 TableHeadItem(i)
|
|
|
456 }
|
|
|
457 cata("</tr>\n")
|
|
|
458 for (i in 1:nrow(sigDiff)) {
|
|
|
459 cata("<tr>\n")
|
|
|
460 TableHeadItem(unmake.names(row.names(sigDiff)[i]))
|
|
|
461 for (j in 1:ncol(sigDiff)) {
|
|
|
462 TableItem(as.character(sigDiff[i, j]))
|
|
|
463 }
|
|
|
464 cata("</tr>\n")
|
|
|
465 }
|
|
|
466 cata("</table>")
|
|
|
467
|
|
|
468 cata("<h4>Plots:</h4>\n")
|
|
|
469 for (i in 1:nrow(linkData)) {
|
|
|
470 if (grepl(".pdf", linkData$Link[i])) {
|
|
|
471 HtmlLink(linkData$Link[i], linkData$Label[i])
|
|
|
472 }
|
|
|
473 }
|
|
|
474
|
|
|
475 cata("<h4>Tables:</h4>\n")
|
|
|
476 for (i in 1:nrow(linkData)) {
|
|
|
477 if (grepl(".tsv", linkData$Link[i])) {
|
|
|
478 HtmlLink(linkData$Link[i], linkData$Label[i])
|
|
|
479 }
|
|
|
480 }
|
|
|
481
|
|
|
482 if (wantRda) {
|
|
|
483 cata("<h4>R Data Object:</h4>\n")
|
|
|
484 for (i in 1:nrow(linkData)) {
|
|
|
485 if (grepl(".rda", linkData$Link[i])) {
|
|
|
486 HtmlLink(linkData$Link[i], linkData$Label[i])
|
|
|
487 }
|
|
|
488 }
|
|
|
489 }
|
|
|
490
|
|
|
491 cata("<p>Alt-click links to download file.</p>\n")
|
|
|
492 cata("<p>Click floppy disc icon associated history item to download ")
|
|
|
493 cata("all files.</p>\n")
|
|
|
494 cata("<p>.tsv files can be viewed in Excel or any spreadsheet program.</p>\n")
|
|
|
495
|
|
|
496 cata("<h4>Additional Information</h4>\n")
|
|
|
497 cata("<ul>\n")
|
|
|
498 if (cpmReq!=0 && sampleReq!=0) {
|
|
|
499 tempStr <- paste("Genes without more than", cpmReq,
|
|
|
500 "CPM in at least", sampleReq, "samples are insignificant",
|
|
|
501 "and filtered out.")
|
|
|
502 ListItem(tempStr)
|
|
|
503 filterProp <- round(filteredCount/preFilterCount*100, digits=2)
|
|
|
504 tempStr <- paste0(filteredCount, " of ", preFilterCount," (", filterProp,
|
|
|
505 "%) genes were filtered out for low expression.")
|
|
|
506 ListItem(tempStr)
|
|
|
507 }
|
|
|
508 ListItem(normOpt, " was the method used to normalise library sizes.")
|
|
|
509 if (wantWeight) {
|
|
|
510 ListItem("Weights were applied to samples.")
|
|
|
511 } else {
|
|
|
512 ListItem("Weights were not applied to samples.")
|
|
|
513 }
|
|
|
514 if (pAdjOpt!="none") {
|
|
|
515 if (pAdjOpt=="BH" || pAdjOpt=="BY") {
|
|
|
516 tempStr <- paste0("MA-Plot highlighted genes are significant at FDR ",
|
|
|
517 "of ", pValReq," and exhibit log2-fold-change of at ",
|
|
|
518 "least ", lfcReq, ".")
|
|
|
519 ListItem(tempStr)
|
|
|
520 } else if (pAdjOpt=="holm") {
|
|
|
521 tempStr <- paste0("MA-Plot highlighted genes are significant at adjusted ",
|
|
|
522 "p-value of ", pValReq," by the Holm(1979) ",
|
|
|
523 "method, and exhibit log2-fold-change of at least ",
|
|
|
524 lfcReq, ".")
|
|
|
525 ListItem(tempStr)
|
|
|
526 }
|
|
|
527 } else {
|
|
|
528 tempStr <- paste0("MA-Plot highlighted genes are significant at p-value ",
|
|
|
529 "of ", pValReq," and exhibit log2-fold-change of at ",
|
|
|
530 "least ", lfcReq, ".")
|
|
|
531 ListItem(tempStr)
|
|
|
532 }
|
|
|
533 cata("</ul>\n")
|
|
|
534
|
|
|
535 cata("<h4>Summary of experimental data:</h4>\n")
|
|
|
536
|
|
|
537 cata("<p>*CHECK THAT SAMPLES ARE ASSOCIATED WITH CORRECT GROUP*</p>\n")
|
|
|
538
|
|
|
539 cata("<table border=\"1\" cellpadding=\"3\">\n")
|
|
|
540 cata("<tr>\n")
|
|
|
541 TableItem()
|
|
|
542 for (i in names(factors)) {
|
|
|
543 TableHeadItem(i)
|
|
|
544 }
|
|
|
545 cata("</tr>\n")
|
|
|
546
|
|
|
547 for (i in 1:nrow(factors)) {
|
|
|
548 cata("<tr>\n")
|
|
|
549 TableHeadItem(row.names(factors)[i])
|
|
|
550 for (j in ncol(factors)) {
|
|
|
551 TableItem(as.character(unmake.names(factors[i, j])))
|
|
|
552 }
|
|
|
553 cata("</tr>\n")
|
|
|
554 }
|
|
|
555 cata("</table>")
|
|
|
556
|
|
|
557 cit <- character()
|
|
|
558 link <- character()
|
|
|
559 link[1] <- paste0("<a href=\"",
|
|
|
560 "http://www.bioconductor.org/packages/release/bioc/",
|
|
|
561 "vignettes/limma/inst/doc/usersguide.pdf",
|
|
|
562 "\">", "limma User's Guide", "</a>.")
|
|
|
563
|
|
|
564 link[2] <- paste0("<a href=\"",
|
|
|
565 "http://www.bioconductor.org/packages/release/bioc/",
|
|
|
566 "vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf",
|
|
|
567 "\">", "edgeR User's Guide", "</a>")
|
|
|
568
|
|
|
569 cit[1] <- paste("Please cite the paper below for the limma software itself.",
|
|
|
570 "Please also try to cite the appropriate methodology articles",
|
|
|
571 "that describe the statistical methods implemented in limma,",
|
|
|
572 "depending on which limma functions you are using. The",
|
|
|
573 "methodology articles are listed in Section 2.1 of the",
|
|
|
574 link[1],
|
|
|
575 "Cite no. 3 only if sample weights were used.")
|
|
|
576 cit[2] <- paste("Smyth, GK (2005). Limma: linear models for microarray data.",
|
|
|
577 "In: 'Bioinformatics and Computational Biology Solutions using",
|
|
|
578 "R and Bioconductor'. R. Gentleman, V. Carey, S. doit,.",
|
|
|
579 "Irizarry, W. Huber (eds), Springer, New York, pages 397-420.")
|
|
|
580 cit[3] <- paste("Please cite the first paper for the software itself and the",
|
|
|
581 "other papers for the various original statistical methods",
|
|
|
582 "implemented in edgeR. See Section 1.2 in the", link[2],
|
|
|
583 "for more detail.")
|
|
|
584 cit[4] <- paste("Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a",
|
|
|
585 "Bioconductor package for differential expression analysis",
|
|
|
586 "of digital gene expression data. Bioinformatics 26, 139-140")
|
|
|
587 cit[5] <- paste("Robinson MD and Smyth GK (2007). Moderated statistical tests",
|
|
|
588 "for assessing differences in tag abundance. Bioinformatics",
|
|
|
589 "23, 2881-2887")
|
|
|
590 cit[6] <- paste("Robinson MD and Smyth GK (2008). Small-sample estimation of",
|
|
|
591 "negative binomial dispersion, with applications to SAGE data.",
|
|
|
592 "Biostatistics, 9, 321-332")
|
|
|
593 cit[7] <- paste("McCarthy DJ, Chen Y and Smyth GK (2012). Differential",
|
|
|
594 "expression analysis of multifactor RNA-Seq experiments with",
|
|
|
595 "respect to biological variation. Nucleic Acids Research 40,",
|
|
|
596 "4288-4297")
|
|
|
597 cit[8] <- paste("Law, CW, Chen, Y, Shi, W, and Smyth, GK (2014). Voom:",
|
|
|
598 "precision weights unlock linear model analysis tools for",
|
|
|
599 "RNA-seq read counts. Genome Biology 15, R29.")
|
|
|
600 cit[9] <- paste("Ritchie, M. E., Diyagama, D., Neilson, J., van Laar,",
|
|
|
601 "R., Dobrovic, A., Holloway, A., and Smyth, G. K. (2006).",
|
|
|
602 "Empirical array quality weights for microarray data.",
|
|
|
603 "BMC Bioinformatics 7, Article 261.")
|
|
|
604 cata("<h3>Citations</h3>\n")
|
|
|
605
|
|
|
606 cata("<h4>limma</h4>\n")
|
|
|
607 cata(cit[1], "\n")
|
|
|
608 cata("<ol>\n")
|
|
|
609 ListItem(cit[2])
|
|
|
610 ListItem(cit[8])
|
|
|
611 ListItem(cit[9])
|
|
|
612 cata("</ol>\n")
|
|
|
613
|
|
|
614 cata("<h4>edgeR</h4>\n")
|
|
|
615 cata(cit[3], "\n")
|
|
|
616 cata("<ol>\n")
|
|
|
617 ListItem(cit[4])
|
|
|
618 ListItem(cit[5])
|
|
|
619 ListItem(cit[6])
|
|
|
620 ListItem(cit[7])
|
|
|
621 cata("</ol>\n")
|
|
|
622
|
|
|
623 cata("<p>Report problems to: su.s@wehi.edu.au</p>\n")
|
|
|
624
|
|
|
625 for (i in 1:nrow(linkData)) {
|
|
|
626 if (grepl("session_info", linkData$Link[i])) {
|
|
|
627 HtmlLink(linkData$Link[i], linkData$Label[i])
|
|
|
628 }
|
|
|
629 }
|
|
|
630
|
|
|
631 cata("<table border=\"0\">\n")
|
|
|
632 cata("<tr>\n")
|
|
|
633 TableItem("Task started at:"); TableItem(timeStart)
|
|
|
634 cata("</tr>\n")
|
|
|
635 cata("<tr>\n")
|
|
|
636 TableItem("Task ended at:"); TableItem(timeEnd)
|
|
|
637 cata("</tr>\n")
|
|
|
638 cata("<tr>\n")
|
|
|
639 TableItem("Task run time:"); TableItem(timeTaken)
|
|
|
640 cata("<tr>\n")
|
|
|
641 cata("</table>\n")
|
|
|
642
|
|
|
643 cata("</body>\n")
|
|
|
644 cata("</html>") |
