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