|
6
|
1 #!/usr/bin/Rscript
|
|
|
2
|
|
|
3 ## --------------------
|
|
|
4 ## prints all arguments as msg
|
|
|
5 ## --------------------
|
|
|
6 catMsg <- function( msg=c() ){
|
|
|
7 cat( MAIN_NAME, paste( msg, collapse="" ), "\n", sep='')
|
|
|
8 }
|
|
|
9 ## --------------------
|
|
|
10 ## return the location of this script
|
|
|
11 ## --------------------
|
|
|
12 getScriptPath <- function(){
|
|
|
13 cmd.args <- commandArgs()
|
|
|
14 m <- regexpr("(?<=^--file=).+", cmd.args, perl=TRUE)
|
|
|
15 script.dir <- dirname(regmatches(cmd.args, m))
|
|
|
16 if( length(script.dir) == 0 ) stop("[ERR] Can't determine script dir: please call the script with Rscript\n")
|
|
|
17 if( length(script.dir) > 1 ) stop("[ERR] Can't determine script dir: more than one '--file' argument detected\n")
|
|
|
18 return(script.dir)
|
|
|
19 }
|
|
|
20 ## --------------------
|
|
|
21 ## Some html creation functions
|
|
|
22 ## --------------------
|
|
|
23 htmlTableRow <- function( string_array=c() ){
|
|
|
24 td_cells <- ''
|
|
|
25 for ( i in string_array ){
|
|
|
26 td_cells <- paste( td_cells, '<td>', i, '</td>', sep='' )
|
|
|
27 }
|
|
|
28 return( paste( "<tr>", td_cells, "</tr>") )
|
|
|
29 }
|
|
|
30 htmlLink <- function( path, desc="LINK" ){
|
|
|
31 return( paste( '<a href="', path, '">', desc, "</a>", sep='') )
|
|
|
32 }
|
|
|
33 ## --------------------
|
|
|
34 ## constructs a list with input bam file info
|
|
|
35 ## --------------------
|
|
|
36 makeBamFileList <- function( paths, names ){
|
|
|
37 tmp <- list()
|
|
|
38 l1 <- length(paths)
|
|
|
39 l2 <- length(names)
|
|
|
40 if ( l1 != l2 ) stop( "Unequal amount of bam-paths (", l1, ") and -names (", l2, ") in makeBamFileList!!!\n" )
|
|
|
41 if ( l1 == 0 ){ return(tmp) } # empty list in debug mode
|
|
|
42
|
|
|
43 for ( i in 1:length(paths) ){
|
|
|
44 path <- paths[i]
|
|
|
45 name <- names[i]
|
|
|
46 file <- basename(path)
|
|
|
47
|
|
|
48 tmp[[ file ]] <- name
|
|
|
49 tmp[[ 'all_paths' ]] <- c( tmp[[ 'all_paths' ]], path )
|
|
|
50 tmp[[ 'all_files' ]] <- c( tmp[[ 'all_files' ]], file )
|
|
|
51 tmp[[ 'all_names' ]] <- c( tmp[[ 'all_names' ]], name )
|
|
|
52 }
|
|
|
53 return( tmp )
|
|
|
54 }
|
|
|
55
|
|
|
56 ## --------------------
|
|
|
57 ## copied code for extracting the regions by segment call status
|
|
|
58 ## --------------------
|
|
|
59 fuseRegions <- function( obj, minRatio=0 ) {
|
|
|
60 if ( ncol(obj) > 1 ) stop('Please specify which sample...')
|
|
|
61
|
|
|
62 data <- data.frame( obj@featureData@data[,1:3], copynumber(obj), segmented(obj), check.names=FALSE, stringsAsFactors=FALSE)
|
|
|
63 colnames( data ) <- c( "chr", "start", "end", "log2", "segmentval" )
|
|
|
64
|
|
|
65 fused.data <- data.frame()
|
|
|
66 curr.bin <- 1
|
|
|
67 for ( chr in unique( data$chr ) ) {
|
|
|
68 chr.data <- data[ data$chr == chr, ]
|
|
|
69 prev.bin <- curr.bin
|
|
|
70 prev.log2 <- chr.data[ 1, 'log2' ]
|
|
|
71 prev.segm <- chr.data[ 1, 'segmentval' ]
|
|
|
72 start <- chr.data[ 1, 'start' ]
|
|
|
73
|
|
|
74 if ( nrow(chr.data) > 1) {
|
|
|
75 for ( i in 2:nrow(chr.data) ) {
|
|
|
76 curr.bin <- curr.bin + 1
|
|
|
77 curr.segm <- chr.data[ i, 'segmentval']
|
|
|
78
|
|
|
79 if ( curr.segm != prev.segm ) {
|
|
|
80 fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], segmentval=round(prev.segm, digits=DECIMALS) ) )
|
|
|
81 prev.segm <- curr.segm
|
|
|
82 prev.bin <- curr.bin
|
|
|
83 start <- chr.data[ i, 'start']
|
|
|
84 }
|
|
|
85 }
|
|
|
86 fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], segmentval=round(prev.segm, digits=DECIMALS) ) )
|
|
|
87 }else{
|
|
|
88 fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], segmentval=round(prev.segm, digits=DECIMALS) ) )
|
|
|
89 }
|
|
|
90 }
|
|
|
91 ## remove regions with low amplitude
|
|
|
92 fused.data <- fused.data[ abs(fused.data$segmentval) >= minRatio, ]
|
|
|
93 fused.data
|
|
|
94 }
|
|
|
95
|
|
|
96 ## DESC: takes the output of fuse.regions and outputs a txt file per sample
|
|
|
97 outputRegionsFromList <- function ( regionsList, outputBasename, outputDir="./", binSize, storeList ){
|
|
|
98 if ( missing(regionsList) ) stop( 'Please provide regionsList...' )
|
|
|
99 if ( missing(outputBasename) ) stop( 'Please provide outputBasename...' )
|
|
|
100 if ( !is.list(regionsList) ) stop( 'Input not a list...?' )
|
|
|
101 if ( length(regionsList) < 1 ) stop( 'List seems empty...?' )
|
|
|
102 if ( file.exists( outputDir ) ) catMsg( c(" Using dir ", outputDir, " for output") )
|
|
|
103 else dir.create( outputDir )
|
|
|
104
|
|
|
105 ## have to set R output options otherwise scientific method is used at some point
|
|
|
106 options( "scipen"=100 )
|
|
|
107
|
|
|
108 sampleCount <- length( regionsList )
|
|
|
109 sampleNames <- names( regionsList )
|
|
|
110 bedgraphColumns <- c( 'chr', 'start', 'end', 'segmentval' )
|
|
|
111
|
|
|
112 catMsg( c( " There are ", sampleCount, " samples found in input list") )
|
|
|
113
|
|
|
114 for ( sample in sampleNames ){
|
|
|
115 catMsg( c(" Working on sample ", sample ) )
|
|
|
116 regionCount <- nrow( regionsList[[sample]] )
|
|
|
117
|
|
|
118 outSampleBase <- paste( outputBasename, '_', sample, '_', binSize, 'kbp', sep='')
|
|
|
119 outBedgraphFile <- paste( outSampleBase, '.bedGraph', sep="" )
|
|
|
120 outBedgraphPath <- paste( outputDir, '/', outBedgraphFile, sep="" )
|
|
|
121
|
|
|
122 ## ---------- BEDGRAPH ----------
|
|
|
123 txt <- paste( "track type=bedGraph color=0,100,0 altColor=255,0,0 name=", sample," description=segmented_regions_from_QDNAseq_",binSize,"kbp\n", sep="")
|
|
|
124 sink( outBedgraphPath )
|
|
|
125 cat( txt )
|
|
|
126 sink()
|
|
|
127 write.table( regionsList[[sample]][,bedgraphColumns], outBedgraphPath, quote=F, sep="\t", row.names=F, append=T, col.names=F)
|
|
|
128 #outFiles[[sample]] <- c( outBedgraphFile )
|
|
|
129 storeList[[ paste( binSize, sample, 'bedgraph', sep="_")]] <- outBedgraphFile
|
|
|
130 }
|
|
|
131 return(storeList)
|
|
|
132 }
|
|
|
133
|
|
|
134
|
|
|
135 ## ==================================================
|
|
|
136 ## Start of analysis
|
|
|
137 ## ==================================================
|
|
|
138 MAIN_NAME <- '[INFO] '
|
|
|
139 TOOL_PATH <- getScriptPath()
|
|
|
140 CSS_FILE <- paste( TOOL_PATH, '/static/css/QDNAseq.css', sep="" )
|
|
|
141 DECIMALS <- 3
|
|
|
142 WEB_LINK <- 'http://www.bioconductor.org/packages/release/bioc/html/QDNAseq.html'
|
|
|
143 PURE_CSS <- 'http://yui.yahooapis.com/pure/0.5.0/pure-min.css'
|
|
|
144
|
|
|
145 catMsg( "Starting QDNAseq wrapper" )
|
|
|
146 catMsg( "Loading R libraries" )
|
|
7
|
147 catMsg( R.version.string )
|
|
|
148
|
|
|
149 catMsg( "Loading Limma" )
|
|
|
150 library( "limma" )
|
|
|
151 catMsg( "Loading DNAcopy" )
|
|
|
152 library( "DNAcopy" )
|
|
|
153
|
|
|
154 ## supress msg to allow R to finish with non-error msg
|
|
6
|
155 suppressWarnings( suppressMessages( library( QDNAseq, quietly = TRUE ) ) )
|
|
|
156 suppressWarnings( suppressMessages( library( CGHcall, quietly = TRUE ) ) )
|
|
|
157
|
|
|
158 ## only one param: the tmp config file
|
|
|
159 cmdLineArgs <- commandArgs(TRUE)
|
|
|
160 config <- cmdLineArgs[1]
|
|
|
161
|
|
|
162 ## sourcing the config file will load all input params
|
|
|
163 ## many variables are imported via sourced "config"
|
|
|
164 source( config )
|
|
|
165
|
|
|
166 ## desparate tries to make png text scale well, damn you R...!
|
|
|
167 PLOT_RES <- min( PLOT_WIDTH, PLOT_HEIGHT ) / 6.3
|
|
|
168 PAR_SET <- list( pch=22 )
|
|
|
169
|
|
|
170 systemUser <- system("whoami",T)
|
|
|
171 qdnaseqVersion <- packageDescription( "QDNAseq" )$Version
|
|
|
172 rVersion <- R.version.string
|
|
|
173 startTime <- Sys.time()
|
|
|
174 analysisStart <- as.character( startTime )
|
|
|
175 catMsg( c("QDNAseq version: ", qdnaseqVersion) )
|
|
|
176 catMsg( c( rVersion ) )
|
|
|
177
|
|
|
178 ## get the comma separated list of chromosomes to exclude
|
|
|
179 excludeChrs <- unlist( strsplit( excludeChrsString, ",") )
|
|
|
180 binSizes <- as.numeric( unlist( strsplit( binSizesString, ",") ) )
|
|
|
181
|
|
|
182
|
|
|
183 ## ------------------------
|
|
|
184 ## DEBUG
|
|
|
185 if ( debug ){
|
|
|
186 catMsg( c("Analysis run by user: ", systemUser ) )
|
|
|
187 catMsg( c("DEBUG SessionInfo: " ) )
|
|
|
188 print( sessionInfo() )
|
|
|
189 }
|
|
|
190 ## /DEBUG
|
|
|
191 ## ------------------------
|
|
|
192
|
|
|
193 ## prepare output dir
|
|
|
194 if ( !file.exists( outputPath) ){
|
|
|
195 dir.create( outputPath )
|
|
|
196 }
|
|
|
197
|
|
|
198 ## copy source config file to output dir to include it in output zip
|
|
|
199 if ( inGalaxy ){
|
|
|
200 file.copy( config, paste(outputPath, 'galaxyConfigFile.R', sep='/') )
|
|
|
201 }
|
|
|
202
|
|
|
203 ## setup bam filelist for easy retrieval later
|
|
|
204 fileList <- makeBamFileList( bamsPaths, bamsNames )
|
|
|
205 bamCount <- length( fileList[[ 'all_paths' ]] )
|
|
|
206
|
|
|
207 gzipOutputName <- paste( 'QDNAseqResults_', outputName, '.zip', sep='' )
|
|
|
208 gzipOutputPath <- paste( outputPath, '/', gzipOutputName, sep='')
|
|
|
209 htmlOutputName <- 'index.html'
|
|
|
210 htmlOutputPath <- paste( outputPath, '/', htmlOutputName, sep='')
|
|
|
211
|
|
|
212 plotted_images <- list() # to keep track of images for later linking
|
|
|
213 regions <- list() # will contain the segments
|
|
|
214 outputFiles <- list()
|
|
|
215
|
|
|
216 ## ------------------------
|
|
|
217 ## in case of debug just use inbuilt LGG data for speedup
|
|
|
218 if ( debug ){
|
|
|
219 binSizes <- c(15)
|
|
|
220 bamsPaths <- c( "BUILD_IN_DATA")
|
|
|
221 bamsNames <- c( "LGG150")
|
|
|
222 fileList <- makeBamFileList( bamsPaths, bamsNames )
|
|
|
223 bamCount <- length( fileList[[ 'all_paths' ]] )
|
|
|
224 }
|
|
|
225
|
|
|
226 for ( binSize in binSizes ){
|
|
|
227
|
|
|
228 ## ------------------------
|
|
|
229 ## construct output file-names and -paths
|
|
|
230 ## ------------------------
|
|
|
231 robjReadCoName <- paste( binSize, 'kbp_QDNAseqReadCounts.rds', sep='')
|
|
|
232 robjCopyNrName <- paste( binSize, 'kbp_QDNAseqCopyNumbers.rds', sep='')
|
|
|
233 igvCopyNrName <- paste( binSize, 'kbp_QDNAseqCopyNumbers.igv', sep='')
|
|
|
234 robjSegmntName <- paste( binSize, 'kbp_QDNAseqCopyNumbersSegmented.rds', sep='')
|
|
|
235 regiOutputName <- paste( binSize, 'kbp_QDNAseqRegions.rds', sep='')
|
|
|
236 noiseImgName <- paste( binSize, 'kbp_QDNAseqNoiseplot.png', sep='')
|
|
|
237
|
|
|
238 robjReadCoPath <- paste( outputPath, '/', robjReadCoName, sep='')
|
|
|
239 robjCopyNrPath <- paste( outputPath, '/', robjCopyNrName, sep='')
|
|
|
240 robjSegmntPath <- paste( outputPath, '/', robjSegmntName, sep='')
|
|
|
241 robjRegionPath <- paste( outputPath, '/', regiOutputName, sep='')
|
|
|
242 igvCopyNrPath <- paste( outputPath, '/', igvCopyNrName, sep='')
|
|
|
243 noiseImgPath <- paste( outputPath, '/', noiseImgName, sep='')
|
|
|
244
|
|
|
245 binAnnFile <- paste( TOOL_PATH, '/static/binannotation/', binSize, 'kbp_binAnnotations.rds', sep="" )
|
|
|
246 if ( file.exists(binAnnFile) ){
|
|
|
247 binAnnotations <- readRDS( binAnnFile )
|
|
|
248 catMsg( c("Using local binAnnotations file" ) )
|
|
|
249 }else{
|
|
|
250 binAnnotations <- getBinAnnotations( binSize=binSize, type=experimentType )
|
|
|
251 }
|
|
|
252
|
|
|
253 ## in case of debug just use inbuilt LGG data for speedup
|
|
|
254 if ( debug ){
|
|
|
255 data(LGG150)
|
|
|
256 readCounts <- LGG150
|
|
|
257 }else{
|
|
|
258 ## provide bamnames because in galaxy everyting is called "dataset_###"
|
|
|
259 readCounts <- binReadCounts( binAnnotations, bamfiles=fileList[[ 'all_paths' ]], bamnames=fileList[[ 'all_names' ]] )
|
|
|
260 }
|
|
|
261
|
|
|
262 readCountsFiltered <- applyFilters( readCounts, residual=TRUE, blacklist=filterBlacklistedBins, mappability=mappabilityCutoff, chromosomes=excludeChrs )
|
|
|
263 readCountsFiltered <- estimateCorrection( readCountsFiltered )
|
|
|
264 copyNumbers <- correctBins( readCountsFiltered )
|
|
|
265 copyNumbersNormalized <- normalizeBins( copyNumbers )
|
|
|
266 copyNumbersSmooth <- smoothOutlierBins( copyNumbersNormalized )
|
|
|
267 sampleNames <- readCountsFiltered@phenoData@data$name
|
|
|
268
|
|
|
269 ## save objects to output dir
|
|
|
270 saveRDS( readCountsFiltered, robjReadCoPath );
|
|
|
271 saveRDS( copyNumbersSmooth, robjCopyNrPath );
|
|
|
272 exportBins( copyNumbersSmooth, file=igvCopyNrPath, format="igv" )
|
|
|
273
|
|
|
274 ## also save objects for galaxy history output if requested
|
|
|
275 if ( doOutputCopynumbersIgv ){
|
|
|
276 #@ a bit hacky galaxy way to allow an unknown number of output files based on param selection
|
|
|
277 #@ see: https://wiki.galaxyproject.org/Admin/Tools/Multiple%20Output%20Files
|
|
|
278 historyName <- paste(binSize, 'kbp-IGV', sep="")
|
|
|
279 igvFile <- paste( newFilePath, "/primary_", outputId, "_", historyName, "_visible_txt", sep="" )
|
|
|
280 exportBins( copyNumbersSmooth, file=igvFile, format="igv" )
|
|
|
281 catMsg( c("Exported igv file to history for ", binSize, "kbp bin") )
|
|
|
282 }
|
|
|
283
|
|
|
284 ## proceed with calling if requested
|
|
|
285 if ( doSegment ){
|
|
|
286 copyNumbersSegmented <- segmentBins( copyNumbersSmooth, undo.splits=undoSplits, undo.SD=undoSD )
|
|
|
287 copyNumbersSegmented <- normalizeSegmentedBins( copyNumbersSegmented )
|
|
|
288 cgh <- makeCgh( copyNumbersSegmented )
|
|
|
289 saveRDS( copyNumbersSegmented, robjSegmntPath );
|
|
|
290 }
|
|
|
291
|
|
|
292 ## ------------------------
|
|
|
293 ## create output files
|
|
|
294 ## ------------------------
|
|
|
295 png( noiseImgPath, width=PLOT_HEIGHT, height=PLOT_HEIGHT, res=PLOT_RES );
|
|
|
296 par( PAR_SET )
|
|
|
297 noisePlot( readCountsFiltered, main=paste( "Noise Plot ", binSize, "kbp", sep=''), col="darkgreen" )
|
|
|
298 dev.off()
|
|
|
299
|
|
|
300 binSize <- as.character( binSize ) # to avoid R using it as array index... (*#$^@ you R!)
|
|
|
301 binSizeString <- paste( binSize, 'kbp', sep='')
|
|
|
302
|
|
|
303 for (i in 1:length(sampleNames) ){
|
|
|
304
|
|
|
305 sample <- sampleNames[i]
|
|
|
306 usedReads <- readCountsFiltered@phenoData@data$used.reads[i]
|
|
|
307 catMsg( c("Creating plots for sample: ", sample, " (", binSizeString, ")" ) )
|
|
|
308
|
|
|
309 type <- 'CopyNumbers'
|
|
|
310 img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png', sep='')
|
|
|
311 img_file_path <- paste( outputPath, '/', img_file, sep='' )
|
|
|
312
|
|
|
313 ## COPYNUMBER PLOT
|
|
|
314 png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES );
|
|
|
315 par( PAR_SET )
|
|
|
316 plot( copyNumbersSmooth[ ,sample ], main=paste(sample, ": CopyNumbers", sep="") )
|
|
|
317 mtext( paste( binSizeString, " bins", sep=""), 3 )
|
|
|
318 abline( h=c(-2,-1,1,2,3,4), lty=1, lwd=0.5, col="grey" )
|
|
|
319 dev.off()
|
|
|
320
|
|
|
321 plotted_images[[ paste(binSize, sample, type, sep="_" ) ]] <- img_file
|
|
|
322
|
|
|
323 if ( doSegment ){
|
|
|
324 type <- 'Segmented'
|
|
|
325 img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png', sep='')
|
|
|
326 img_file_path <- paste( outputPath, '/', img_file, sep='' )
|
|
|
327
|
|
|
328 ## COPYNUMBER PLOT
|
|
|
329 png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT, res=PLOT_RES );
|
|
|
330 par( PAR_SET )
|
|
|
331 plot( copyNumbersSegmented[ ,sample ], main=paste(sample, ": CopyNumbers and Segments", sep="") )
|
|
|
332 mtext( paste( "(", binSizeString, " bins)", sep=""), 3 )
|
|
|
333 abline( h=c(-2,-1,1,2,3,4), lty=1, lwd=0.5, col="grey" )
|
|
|
334 dev.off()
|
|
|
335
|
|
|
336 plotted_images[[ paste(binSize, sample, type, sep="_" ) ]] <- img_file
|
|
|
337
|
|
|
338 catMsg( c(" Fusing regions of sample: ", sample) )
|
|
|
339 regions[[ sample ]] <- fuseRegions( cgh[, sample] )
|
|
|
340
|
|
|
341 region_count <- nrow( data.frame( regions[[ sample ]] ) )
|
|
|
342 catMsg( c( ' sample "', sample, '" has ', region_count, " regions" ) )
|
|
|
343 plotted_images[[ paste(binSize, sample, 'region_count', sep="_" ) ]] <- region_count
|
|
|
344 }
|
|
|
345
|
|
|
346 ## add USED read counts
|
|
|
347 plotted_images[[ paste(binSize, sample, 'usedReads', sep="_" ) ]] <- usedReads
|
|
|
348 }
|
|
|
349
|
|
|
350 if ( doSegment ){
|
|
|
351 saveRDS( regions, robjRegionPath )
|
|
|
352 plotted_images <- outputRegionsFromList( regions, outputBasename=outputName, outputDir=outputPath, binSize=binSize, storeList=plotted_images )
|
|
|
353 }
|
|
|
354 }# end bin
|
|
|
355
|
|
|
356
|
|
|
357 ## ----- debug -----
|
|
|
358 #catMsg( "done" )
|
|
|
359 #q(status=0)
|
|
|
360 ## ---- /debug -----
|
|
|
361
|
|
|
362
|
|
|
363 ## ------------------------
|
|
|
364 ## prepare output
|
|
|
365 ## ------------------------
|
|
|
366 catMsg( "...zipping output")
|
|
|
367 zip_cmd <- paste( "zip -j", gzipOutputPath, paste(outputPath,'/*',sep='') ) ## -j is for removing dirs from the tree
|
|
|
368 system( zip_cmd )
|
|
|
369
|
|
|
370 ## ------------------------
|
|
|
371 ## get filesizes for report
|
|
|
372 ## ------------------------
|
|
|
373 zippedSize <- paste( round( file.info( gzipOutputPath )[["size"]] / 1e+6, digits=2 ), 'MB' )
|
|
|
374 endTime <- Sys.time()
|
|
|
375 timeDiff <- format( round( endTime - startTime, 3 ) )
|
|
|
376 analysisEnd <- as.character( endTime )
|
|
|
377
|
|
|
378 ## ------------------------
|
|
|
379 ## creating html output to be linked to from the middle galaxy pane
|
|
|
380 ## ------------------------
|
|
|
381 sink( file = htmlOutputPath, type = "output" )
|
|
|
382 cat( "<html>\n")
|
|
|
383 cat( "<head>\n")
|
|
|
384
|
|
|
385 cat( "\t", '<title>QDNAseq Report | ', outputName,'</title>', "\n", sep='' )
|
|
|
386 cat( "\t", '<link rel="stylesheet" href="', PURE_CSS, '">', "\n", sep='' )
|
|
|
387 cat( "\t<style>\n", sep='')
|
|
|
388 ## include CSS into html file, makes it more portable
|
|
|
389 cat( "\t\t", readLines( CSS_FILE ), sep="\n\t\t" )
|
|
|
390 #cat( "\t\th1 {color:red;}", "\n")
|
|
|
391 cat( "\n\t</style>\n" )
|
|
|
392
|
|
|
393 cat( "\n</head>\n")
|
|
|
394 cat( "\n<body>\n")
|
|
|
395
|
|
|
396 cat( "<h1>QDNAseq Report</h1>", "\n")
|
|
|
397
|
|
|
398 cat( '<h3 class="qdnaseq">About this analysis</h3>', "\n")
|
|
|
399 cat( '<p>This page provides access to all results. To have a local copy of this report just download the <a href="', gzipOutputName, '" class="button">zipfile</a> with all output (', zippedSize, ')</p>', "\n", sep='')
|
|
|
400
|
|
|
401 ## ------------------------
|
|
|
402 ## table with general info
|
|
|
403 ## ------------------------
|
|
|
404 cat( '<h3 class="qdnaseq">Settings</h3><p>', "\n")
|
|
|
405 cat( '<table class="pure-table pure-table-striped"><tbody>' )
|
|
|
406 cat( htmlTableRow( c( "AnalysisName", outputName ) ) )
|
|
|
407 cat( htmlTableRow( c( "AnalysisStart", analysisStart ) ) )
|
|
|
408 cat( htmlTableRow( c( "AnalysisEnd", analysisEnd ) ) )
|
|
|
409 cat( htmlTableRow( c( "AnalysisTime", timeDiff ) ) )
|
|
|
410 cat( htmlTableRow( c( "BinSizes (kbp)", paste(binSizes,collapse=", ") ) ) )
|
|
|
411 cat( htmlTableRow( c( "R info", rVersion ) ) )
|
|
|
412 cat( htmlTableRow( c( "QDNAseq info", qdnaseqVersion ) ) )
|
|
|
413
|
|
|
414 sampleStrings <- c()
|
|
|
415 for ( galaxyName in fileList[[ 'all_files' ]] ){
|
|
|
416 sampleName <- fileList[[ galaxyName ]]
|
|
|
417 sampleStrings <- c( sampleStrings, paste( galaxyName, ' (', sampleName, ')', sep='' ) )
|
|
|
418 }
|
|
|
419 cat( htmlTableRow( c( "InputBams", paste( sampleStrings, collapse=", ") ) ) )
|
|
|
420
|
|
|
421 cat( "</tbody></table></p>", "\n")
|
|
|
422
|
|
|
423 ## ------------------------
|
|
|
424 ## list with links to all output files
|
|
|
425 ## ------------------------
|
|
|
426 cat( '<h3 class="qdnaseq">Output files</h3><p>', "\n")
|
|
|
427 cat( '<p>This table contains output files that can be used for local downstream analysis with the bioconductor QDNAseq package. For each bin-size / data-level there is a R data structure file with data of all samples. See ', htmlLink( WEB_LINK, 'the bioconductor QDNAseq documentation' ), ' for more information on how to work with these files</p>', "\n", sep='')
|
|
|
428 cat( '<table class="pure-table pure-table-striped">', "\n" )
|
|
|
429 cat( '<thead><th>Type</th>', as.vector( mapply( paste, "<th>", binSizes, "kbp</th>", sep="" ) ),'</thead>', "\n" )
|
|
|
430 cat( "<tbody>", "\n")
|
|
|
431 files <- list()
|
|
|
432 fileTypes <- c( 'ReadCounts.rds', 'CopyNumbers.rds' )
|
|
|
433 if ( doSegment ){ fileTypes <- c( fileTypes, 'CopyNumbersSegmented.rds') }
|
|
|
434
|
|
|
435 for ( fileType in fileTypes ){
|
|
|
436 fileNames <- mapply( paste, binSizes, paste( 'kbp_QDNAseq', fileType, sep=''), sep='')
|
|
|
437 fileLinks <- mapply( htmlLink, fileNames, paste( binSizes, "kbp", sep="" ) )
|
|
|
438 cat( htmlTableRow( c( fileType, fileLinks ) ) )
|
|
|
439 }
|
|
|
440 cat( "\n</tbody></table></p>", "\n")
|
|
|
441
|
|
|
442 ## ------------------------
|
|
|
443 ## table with links to files
|
|
|
444 ## ------------------------
|
|
|
445 ratio <- PLOT_WIDTH / PLOT_HEIGHT
|
|
|
446 width <- 960; height <- width / ratio ## bigger img
|
|
|
447 width_t <- 100; height_t <- 40 ## thumb img
|
|
|
448
|
|
|
449 cat( '<h3 class="qdnaseq">Results: overview</h3><p>', "\n")
|
|
|
450 cat( '<p>This table contains the visual results of the copy number aberration analysis. You can click on an image to jump to the larger version. If segmentation was performed as well the number of segments is shown and a file with genomic regions can be downloaded (just remember to inspect the results carefully as this is a more exprimental step).</p>', "\n", sep='')
|
|
|
451 plots_html <- ''
|
|
|
452
|
|
|
453 colspan <- 1
|
|
|
454 binHeader <- "<th>Image</th>"
|
|
|
455 if ( doSegment ){ # extra column with segment info
|
|
|
456 colspan <- 2
|
|
|
457 binHeader <- "<th>Image</th><th>Segments</th>"
|
|
|
458 }
|
|
|
459 cat( '<table class="pure-table pure-table-striped">', "\n" )
|
|
|
460 cat( '<thead><tr><th></th><th></th>', as.vector( mapply( paste, "<th colspan=\"", colspan,"\">", binSizes, "kbp</th>", sep="" ) ), '</tr></thead>' )
|
|
|
461 cat( '<thead><tr><th>Sample / File</th><th>Reads</th>', rep( binHeader, length(binSizes) ), '</tr></thead>' )
|
|
|
462 cat( '<tbody>' )
|
|
|
463
|
|
|
464 for ( bam_file in bamsNames ){
|
|
|
465
|
|
|
466 usedReads <- plotted_images[[ paste(binSize, bam_file, 'usedReads', sep="_" ) ]]
|
|
|
467 usedReads <- format( as.integer(usedReads), digits=4, decimal.mark=".", big.mark="," )
|
|
|
468 htmlRow <- paste( '<tr><td>', bam_file, '</td><td>', usedReads, '</td>', sep='' )
|
|
|
469
|
|
|
470 for ( binSize in binSizes ){
|
|
|
471
|
|
|
472 ## add thumbnails to table with links to anchors on html page
|
|
|
473 copy_img <- plotted_images[[ paste(binSize, bam_file, 'CopyNumbers', sep="_" ) ]]
|
|
|
474 html_copy_thumb <- htmlLink( path=paste('#', copy_img, sep=''), paste('<img src="',copy_img,'" alt="', bam_file, '" width="', width_t, '" height="', height_t, '">', sep='') )
|
|
|
475 html_copy_img <- htmlLink( path=copy_img, paste('<img id="', copy_img,'" src="',copy_img,'" alt="',bam_file, '" width="', width, '" height="', height, '">', sep='') )
|
|
|
476 html_segm_img <- ''
|
|
|
477 html_bedGraph <- ''
|
|
|
478 region_count <- ''
|
|
|
479 htmlRow <- paste( htmlRow, '<td>', html_copy_thumb, '</td>' )
|
|
|
480
|
|
|
481 if ( doSegment ){
|
|
|
482 segm_img <- plotted_images[[ paste(binSize, bam_file, 'Segmented', sep="_" ) ]]
|
|
|
483 region_count <- plotted_images[[ paste(binSize, bam_file, 'region_count', sep="_" ) ]]
|
|
|
484
|
|
|
485 html_bedGraph <- htmlLink( path=plotted_images[[ paste(binSize, bam_file, 'bedgraph', sep="_" ) ]], 'bedGraph' )
|
|
|
486 html_segm_img <- htmlLink( path=segm_img, paste('<img id="', segm_img,'" src="', segm_img,'" alt="', bam_file, '" width="', width, '" height="', height,'">', sep='') )
|
|
|
487 htmlRow <- paste( htmlRow, '<td>', region_count, ' (', html_bedGraph, ')</td>', sep="" )
|
|
|
488 }
|
|
|
489 plots_html <- paste( plots_html, html_copy_img, "\n", html_segm_img, "\n<br \\>\n", sep='' )
|
|
|
490 }
|
|
|
491 plots_html <- paste( plots_html, "\n<hr \\>\n", sep='' )
|
|
|
492 ## add info to overview table, including small thumbnails
|
|
|
493 htmlRow <- paste( htmlRow, '</tr>', sep='' )
|
|
|
494 cat( htmlRow, "\n" )
|
|
|
495 }
|
|
|
496 cat( "</tbody></table></p>", "\n")
|
|
|
497
|
|
|
498 ## ------------------------
|
|
|
499 ## section with various output shown
|
|
|
500 ## ------------------------
|
|
|
501 cat( '<h3 class="qdnaseq">Results: Sample plots</h3><p>', "\n")
|
|
|
502 ## now include (large) images in html page
|
|
|
503 cat( plots_html, "\n")
|
|
|
504 cat( "\n</p></body>\n")
|
|
|
505 cat( "\n</html>\n")
|
|
|
506 sink()
|
|
|
507
|
|
|
508 ## ------------------------
|
|
|
509 ## creating main html output for galaxy history
|
|
|
510 ## ------------------------
|
|
|
511 if ( inGalaxy ){ # dont create when running outside Galaxy
|
|
|
512 sink( file = outputHtml, type = "output" )
|
|
|
513
|
|
|
514 cat( "<head>", "\n")
|
|
|
515 cat( "\t", '<link rel="stylesheet" href="', PURE_CSS, '">', "\n", sep='' )
|
|
|
516
|
|
|
517 cat( "<style>", "\n")
|
|
|
518 ## include CSS directly into html file
|
|
|
519 cat( paste( "\t", '/* the css here originates from ', CSS_FILE,' */', "\n") )
|
|
|
520 cat( paste( "\t", readLines( CSS_FILE, n = -1)), sep="\n" )
|
|
|
521 cat( "</style>", "\n")
|
|
|
522 cat( "</head>", "\n")
|
|
|
523
|
|
|
524 cat( '<h1>QDNAseq Results (', outputName,')</h1>', "\n", sep="")
|
|
|
525 cat( '<p>Explore <a href="', htmlOutputName, '" class="button">the results</a> directly within galaxy</p>', "\n", sep="")
|
|
|
526 cat( '<p>Or download a <a href="', gzipOutputName, '" class="button">zipfile</a> with all output (', zippedSize, ')</p>', "\n", sep="" )
|
|
|
527
|
|
|
528 sink()
|
|
|
529 }
|
|
|
530
|
|
|
531 ## ------------------------
|
|
|
532 ## create final zip and quit with status 0 to tell galaxy all was fine
|
|
|
533 ## ------------------------
|
|
|
534 catMsg( "zipping all output")
|
|
|
535 system( paste( "zip -j ", gzipOutputPath, paste(outputPath,'/', htmlOutputName, sep='') ) )
|
|
|
536 catMsg( "done" )
|
|
|
537 q(status=0)
|
