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