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