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