diff QDNAseq.R @ 2:336697c6f7fa draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/QDNAseq.R	Fri Jun 13 09:42:10 2014 -0400
@@ -0,0 +1,548 @@
+#!/usr/bin/Rscript
+
+## --------------------
+## prints all arguments as msg
+## --------------------
+catMsg <- function( msg=c() ){	
+	cat( MAIN_NAME, paste( msg, collapse="" ), "\n", sep='')
+}
+## --------------------
+## return the directory this script exists
+## --------------------
+getScriptPath <- function(){
+    cmd.args <- commandArgs()
+    m <- regexpr("(?<=^--file=).+", cmd.args, perl=TRUE)
+    script.dir <- dirname(regmatches(cmd.args, m))
+    if(length(script.dir) == 0) stop("[ERR] Can't determine script dir: please call the script with Rscript\n")
+    if(length(script.dir) > 1) stop("[ERR] Can't determine script dir: more than one '--file' argument detected\n")
+    return(script.dir)
+}
+## --------------------
+## Some html helper functions
+## --------------------
+htmlTableRow <- function( string_array=c() ){
+	td_cells <- ''
+	for ( i in string_array ){ 
+		td_cells <- paste( td_cells, '<td>', i, '</td>', sep='' )	
+	}
+	return( paste( "<tr>", td_cells, "</tr>") )
+}
+htmlLink <- function( path, desc="LINK" ){
+	return( paste( '<a href="', path, '">', desc, "</a>", sep='') )
+}
+## --------------------
+## constructs a list with bam file info
+## --------------------
+makeBamFileList <- function( paths, names ){	
+	tmp <- list()
+	l1 <- length(paths)
+	l2 <- length(names)
+	if ( l1 != l2 ) stop( "Unequal amount of bam-paths (",l1,") and -names (",l2,") in makeBamFileList!!!\n" )
+	for ( i in 1:length(paths) ){
+		path <- paths[i]
+		name <- names[i]
+		file <- basename(path)
+
+		tmp[[ file ]] <- name
+		tmp[[ 'all_paths' ]] <- c( tmp[[ 'all_paths' ]], path )
+		tmp[[ 'all_files' ]] <- c( tmp[[ 'all_files' ]], file )
+		tmp[[ 'all_names' ]] <- c( tmp[[ 'all_names' ]], name )
+	}
+	return( tmp )
+}
+
+## --------------------
+## took a function from Matias/MarkVdWuel? for extracting the regions by call
+## --------------------
+fuse.regions_test <- function(cgh, onlyCalled=T) {
+	if (ncol(cgh) > 1) stop('Please specify which sample...')
+
+	x <- data.frame(cgh@featureData@data[,1:3], calls(cgh), copynumber(cgh), segmented(cgh), check.names=FALSE, stringsAsFactors=FALSE)
+	colnames( x ) <- c( "chr", "start", "end", "call", "log2", "segmentval" )
+	
+	fused.data <- data.frame()
+	curr.bin <- 1
+	for (chr in unique(x$chr)) {
+		
+		chr.data <- x[x$chr == chr,]
+		prev.bin <- curr.bin
+		prev.call <- chr.data[1, 'call']
+		prev.log2 <- chr.data[1, 'log2']
+		prev.segm <- chr.data[1, 'segmentval']
+		start <- chr.data[1, 'start']
+
+		if ( nrow(chr.data) > 1) {
+			for ( i in 2:nrow(chr.data) ) {
+
+				curr.bin <- curr.bin + 1
+				curr.call <- chr.data[ i, 'call']
+				curr.segm <- chr.data[ i, 'segmentval']
+
+				if ( curr.segm != prev.segm ) {
+					
+					fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], call=prev.call, segmentval=round(prev.segm, digits=DECIMALS) ) )
+					if ( prev.call != 0 ){
+						cat( MAIN_NAME, " ...found called/segmented region (", chr, ':', start, ' call=', prev.call, ' segment=', prev.segm, ")\n", sep="" )
+					}
+					prev.call <- curr.call
+					prev.segm <- curr.segm
+					prev.bin <- curr.bin
+					start <- chr.data[ i, 'start']
+				}
+			}
+			fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], call=prev.call, segmentval=round(prev.segm, digits=DECIMALS) ) )
+		}else {
+			fused.data <- rbind( fused.data, data.frame( chr=chr, start=start, end=chr.data[ i-1, 'end'], call=prev.call, segmentval=round(prev.segm, digits=DECIMALS) ) )
+		}
+	}
+	if ( onlyCalled == T ){
+		fused.data <- fused.data[ fused.data$call != 0, ]
+	}
+	fused.data
+}
+
+## DESC: takes the output of fuse.regions and outputs a txt file per sample
+outputRegionsFromList <- function ( regionsList, outputBasename, outputDir="./" ){
+	if ( missing(regionsList) ) stop( 'Please provide regionsList...' )
+	if ( missing(outputBasename) ) stop( 'Please provide outputBasename...' )
+	if ( !is.list(regionsList) ) stop( 'Input not a list...?' )
+	if ( length(regionsList) < 1 ) stop( 'List seems empty...?' )
+	if ( file.exists( outputDir ) ) cat( MAIN_NAME, " Using dir ", outputDir, " for output\n", sep="")
+	else dir.create( outputDir )
+	outFiles <- list()
+
+	## have to set R output options otherwise scientific method is used at some point
+	options( "scipen"=100 )
+
+	sampleCount <- length( regionsList )
+	sampleNames <- names( regionsList )
+	systemUser <- system("whoami",T)
+	bedgraphColumns <- c( 'chr', 'start', 'end', 'segmentval' )
+
+	cat( MAIN_NAME, " Hello ", systemUser, "!!\n", sep="")
+	cat( MAIN_NAME, " There are ", sampleCount, " samples found in input list...\n", sep="")
+
+	for ( sample in sampleNames ){		
+		cat( MAIN_NAME, " Working on sample ", sample, "\n", sep="")
+		regionCount <- nrow( regionsList[[sample]] )
+		
+		outSampleBase   <- paste( outputBasename, '_', sample, '_QDNAseqRegions', sep='')
+		outBedFile      <- paste( outSampleBase, '.bed', sep="" )
+		outBedPath      <- paste( outputDir, '/', outBedFile, sep="" )
+		outBedgraphFile <- paste( outSampleBase, '.bedGraph', sep="" )
+		outBedgraphPath <- paste( outputDir, '/', outBedgraphFile, sep="" )
+
+		## ---------- BED ----------
+		txt <- "#"
+		sink( outBedPath )
+			cat( txt )
+		sink()
+		write.table( regionsList[[sample]], outBedPath, quote=F, sep="\t", row.names=F, append=T)
+
+		## ---------- BEDGRAPH ----------
+		txt <- paste( "track type=bedGraph color=0,100,0 altColor=255,0,0 name=", sample,"_segmReg description=segmented_regions_from_QDNAseq\n", sep="")
+		sink( outBedgraphPath )
+			cat( txt )
+		sink()
+		write.table( regionsList[[sample]][,bedgraphColumns], outBedgraphPath, quote=F, sep="\t", row.names=F, append=T, col.names=F)
+		outFiles[[sample]] <- c( outBedFile, outBedgraphFile )
+	}
+	outFiles
+}
+
+printIgvFile <- function( dat, filename ){
+	
+	if ( 'calls' %in% assayDataElementNames(dat) ) {
+	 	#output <- paste(output, '-called.igv', sep="")
+	  	cat('#type=COPY_NUMBER\n#track coords=1\n', file=filename)
+	  	df <- data.frame(chromosome=as.character(chromosomes(dat)), start=bpstart(dat), end=bpend(dat), feature=featureNames(dat), calls(dat), check.names=FALSE, stringsAsFactors=FALSE)
+	}else {
+	  #output <- paste(output, '-normalized.igv', sep="")
+	  cat('#type=COPY_NUMBER\n#track coords=1\n', file=filename)
+	  df <- data.frame(chromosome=as.character(chromosomes(dat)), start=bpstart(dat), end=bpend(dat), feature=featureNames(dat), round(copynumber(dat), digits=2), check.names=FALSE, stringsAsFactors=FALSE)
+	}
+
+	df$chromosome[df$chromosome == '23'] <- 'X'
+	df$chromosome[df$chromosome == '24'] <- 'Y'
+	df$chromosome[df$chromosome == '25'] <- 'MT'
+	#return( df )
+	write.table( df, file=filename, append=TRUE, quote=FALSE, sep='\t', row.names=FALSE )	
+}
+
+
+## ==================================================
+## Start of analysis
+## ==================================================
+TOOL_PATH <- getScriptPath()
+MAIN_NAME <- '[INFO] '
+CSS_FILE  <- paste( TOOL_PATH, '/QDNAseq.css', sep="" )
+DECIMALS  <- 3
+WEB_LINK  <- 'http://www.bioconductor.org/packages/release/bioc/html/QDNAseq.html'
+
+catMsg( "Starting QDNAseq wrapper" )	
+catMsg( "Loading R libraries" )
+suppressWarnings( suppressMessages( library( QDNAseq, quietly = TRUE ) ) )
+suppressWarnings( suppressMessages( library( CGHcall, quietly = TRUE ) ) )
+suppressWarnings( suppressMessages( library( matrixStats, quietly = TRUE ) ) )
+
+## only one param: the tmp config file
+cmdLineArgs <- commandArgs(TRUE)
+config      <- cmdLineArgs[1]
+
+## sourcing the config file will load all input params
+source( config )
+
+## if call output requested, set doCall such that we will segment and call
+if ( doOutputCallsRds == TRUE ){ doCall <- TRUE }
+
+## get the comma separated list of chromosomes to exclude
+excludeChrs <- unlist( strsplit( excludeChrsString, ",") )
+
+## ------------------------
+## DEBUG
+#catMsg( "PARAM" )
+#catMsg( galaxy_path )
+#catMsg( repos_path )
+#catMsg( instal_path )
+## /DEBUG
+## ------------------------
+
+## setup bam filelist for easy retrieval later
+catMsg( "Setting up input bam list" )
+cat( bamsPaths, "\n" )
+catMsg( "Namews" )
+cat( bamsNames, "\n" )
+fileList    <- makeBamFileList( bamsPaths, bamsNames )
+bamCount    <- length( fileList[[ 'all_paths' ]] )
+
+## help msg still needs work!
+if ( length(cmdLineArgs) == 0 || cmdLineArgs[1] == "-h" || cmdLineArgs[1] == "--help"){ 
+	cat( paste( MAIN_NAME, "Usage: ", params_help, sep=''), "\n" )
+	quit( save = 'no', status=0, runLast=F )
+}
+
+if ( !file.exists( outputPath) ){
+	dir.create( outputPath )
+}
+
+## because we circumvent params that galaxy can save, we want to
+## copy source config file to output dir to include it in output zip
+file.copy( config, paste(outputPath, 'qdnaseq_config_file.R', sep='/') )
+
+## ------------------------
+## construct output file-names and -paths
+## ------------------------
+htmlOutputName <- 'index.html'
+gzipOutputName <- paste( 'QDNAseqResults_', outputName, '.zip', sep='' )
+robjReadCoName <- paste( binSize, 'kbp_QDNAseqReadCounts.rds', sep='')
+robjCopyNrName <- paste( binSize, 'kbp_QDNAseqCopyNumbers.rds', sep='')
+robjCalledName <- paste( binSize, 'kbp_QDNAseqCalledSegments.rds', sep='')
+regiOutputName <- paste( binSize, 'kbp_QDNAseqRegions.rds', sep='')
+igvCalledName  <- paste( binSize, 'kbp_QDNAseq-calls.igv', sep='')
+igvCopyNrName  <- paste( binSize, 'kbp_QDNAseq-normalized.igv', sep='')
+
+gzipOutputPath <- paste( outputPath, '/', gzipOutputName, sep="")
+htmlOutputPath <- paste( outputPath, '/', htmlOutputName, sep="")
+robjReadCoPath <- paste( outputPath, '/', robjReadCoName, sep="")
+robjCopyNrPath <- paste( outputPath, '/', robjCopyNrName, sep="")
+robjCalledPath <- paste( outputPath, '/', robjCalledName, sep="")
+robjRegionPath <- paste( outputPath, '/', regiOutputName, sep="")
+igv_calledPath <- paste( outputPath, '/', igvCalledName, sep="")
+igv_copyNrPath <- paste( outputPath, '/', igvCopyNrName, sep="")
+
+
+## ------------------------
+## performing QDNAseq analysis steps
+## ------------------------
+if ( debug ){
+	## in case of debug just use inbuilt LGG data for speedup
+	data(LGG150)
+	readCounts <- LGG150
+}else{
+	if ( nchar(binAnnotations) == 0 ){
+	binAnnotations <- getBinAnnotations( binSize=binSize, type=experimentType )
+	}else{
+		## if user provided file, check if correct class
+		if ( class(binAnnotations)[1] != 'AnnotatedDataFrame' ){
+			stop( "Provided binAnnotations file is not of class 'AnnotatedDataFrame'\n" )
+		}
+		binAnnotations <- readRDS( binAnnotations )
+	}
+	## provide bamnames because in galaxy everyting is called "dataset_###"
+	readCounts <- binReadCounts( binAnnotations, bamfiles=fileList[[ 'all_paths' ]], bamnames=bamsNames )
+}
+
+readCountsFiltered    <- applyFilters( readCounts, residual=TRUE, blacklist=filterBlacklistedBins, mappability=mappabilityCutoff, chromosomes=excludeChrs )
+readCountsFiltered    <- estimateCorrection( readCountsFiltered )
+copyNumbers           <- correctBins( readCountsFiltered )
+copyNumbersNormalized <- normalizeBins( copyNumbers )
+copyNumbersSmooth     <- smoothOutlierBins( copyNumbersNormalized )
+
+## save objects to output dir
+saveRDS( readCounts, robjReadCoPath );
+saveRDS( copyNumbersSmooth, robjCopyNrPath );
+#printIgvFile( copyNumbersSmooth, igvCopyNrName )
+#exportBins(copyNumbersSmooth, file=igvCopyNrName, format="igv")
+
+## also save objects for galaxy history output if requested
+if ( doOutputReadcountsRds ){
+	saveRDS( readCountsFiltered, readCountsDatasetFile );	
+}
+if ( doOutputCopynumbersRds ){
+	saveRDS( copyNumbersSmooth, copyNumbersDatasetFile );
+}
+
+## proceed with calling if requested
+if ( doCall ){
+	copyNumbersSegmented  <- segmentBins( copyNumbersSmooth, undo.splits=undoSplits, undo.SD=undoSD )
+	copyNumbersSegmented  <- normalizeSegmentedBins( copyNumbersSegmented )
+	copyNumbersCalled     <- callBins( copyNumbersSegmented )
+	cgh <- makeCgh( copyNumbersCalled )
+	saveRDS( copyNumbersCalled, robjCalledPath );
+	if ( doOutputCallsRds ){
+		saveRDS( copyNumbersCalled, calledSegmentsDatasetFile );
+	}
+	#df <- getIGVdf( copyNumbersCalled )
+	#printIgvFile( copyNumbersCalled, igvCalledName )
+	#write.table( df, file=igvCalledName, append=TRUE, quote=FALSE, sep='\t', row.names=FALSE )
+}
+
+sampleNames <- readCountsFiltered@phenoData@data$name
+
+## ------------------------
+## create output files
+## ------------------------
+plotted_images <- list() # to keep track of images for later linking
+regions <- list() # will contain the (called) segments
+
+noise_img_file <- paste( binSize, 'kbp_QDNAseqNoisePlot.png',  sep='')
+noise_img_file_path <- paste( outputPath, '/', noise_img_file, sep='' )
+png( noise_img_file_path, width=480, height=480 ); 
+	noisePlot( readCountsFiltered )
+dev.off()
+
+for (i in 1:length(sampleNames) ){
+#for (sample in sampleNames(copyNumbersSmooth) ){
+	sample <- sampleNames[i]
+	usedReads  <- readCountsFiltered@phenoData@data$used.reads[i]
+	catMsg( c("Creating plots for sample: ", sample ) )	
+
+	type <- 'CopyNumbers'
+	img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png',  sep='')
+	img_file_path <- paste( outputPath, '/', img_file, sep='' )
+	png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT ); plot( copyNumbersSmooth[ ,sample ] ); dev.off()
+	plotted_images[[ sample ]][[ type ]] <- img_file
+
+	if ( doCall ){
+		type <- 'Called'
+		img_file <- paste( sample, '_', binSize, 'kbp_QDNAseq', type, '.png',  sep='')
+		img_file_path <- paste( outputPath, '/', img_file, sep='' )
+		png( img_file_path, width=PLOT_WIDTH, height=PLOT_HEIGHT ); 
+			plot( copyNumbersCalled[ ,sample ] ); 
+		dev.off()
+		plotted_images[[ sample ]][[ type ]] <- img_file
+
+		cat( MAIN_NAME, " Fusing regions of sample: ", sample, "\n", sep="")
+		regions[[ sample ]] <- fuse.regions_test( cgh[, sample] )
+		region_count <- nrow( data.frame( regions[[ sample ]] ) )
+		cat( MAIN_NAME, ' sample "', sample, '" has ', region_count, " regions\n", sep="" )
+		plotted_images[[ sample ]][[ 'region_count' ]] <- region_count
+	}
+
+	## add read counts
+	catMsg( "Used")
+	catMsg( usedReads )
+
+	plotted_images[[ sample ]][[ 'usedReads'  ]] <- usedReads
+}
+#cat( MAIN_NAME, "PLOTTED_IMAGES: ", names(plotted_images), "\n", sep="" )
+
+if ( doCall ){
+	saveRDS( regions, robjRegionPath )
+	printedFiles <- outputRegionsFromList( regions, outputBasename=outputName, outputDir=outputPath )	
+}
+
+## ------------------------
+## prepare output
+## ------------------------
+cat( MAIN_NAME, "...zipping output\n")
+zip_cmd <- paste( "zip -j", gzipOutputPath, paste(outputPath,'/*',sep='') ) ## -j is for removing dirs from the tree
+system( zip_cmd )
+
+## ------------------------
+## get filesizes for report
+## ------------------------
+zippedSize <- paste( round( file.info( gzipOutputPath )[["size"]] / 1000000, digits=2 ), 'MB' )
+readCoSize <- paste( round( file.info( robjReadCoPath )[["size"]] / 1000000, digits=2 ), 'MB' )
+copyNrSize <- paste( round( file.info( robjCopyNrPath )[["size"]] / 1000000, digits=2 ), 'MB' )
+calledSize <- paste( round( file.info( robjCalledPath )[["size"]] / 1000000, digits=2 ), 'MB' )
+regionSize <- paste( round( file.info( robjRegionPath )[["size"]] / 1000000, digits=2 ), 'MB' )
+
+## ------------------------
+## creating html output to be linked to from the middle galaxy pane
+## ------------------------
+
+#sink( file = outputHtml, type = "output" )
+sink( file = htmlOutputPath, type = "output" )
+		cat( "<html>\n")
+		cat( "<head>\n")
+			#cat( '<link rel="stylesheet" type="text/css" href="test.css" media="screen" />', "\n" )
+			#cat( '<link rel="stylesheet" type="text/css" href="../../../../static/style/test.css" media="screen" />', 
+			cat( "\t", '<link rel="stylesheet" href="http://yui.yahooapis.com/pure/0.4.2/pure-min.css">', "\n" )
+			cat( "\t<style>\n")
+				## have to include CSS into html file, because css referencing outside own dir doesn't seem to work...
+				cat( paste( "\t\t", '/* the css here originates from ', CSS_FILE,' */', "\n") )
+				cat( paste( "\t\t", readLines( CSS_FILE, n = -1)), sep="\n" )
+				#cat( "\t\th1 {color:red;}", "\n")
+
+			cat( "\t</style>\n" )
+			
+		cat( "\n</head>\n")
+		cat( "\n<body>\n")
+
+		cat( "<h1>QDNAseq Report</h1>", "\n")
+		
+		cat( '<h3 class="qdnaseq">About this analysis</h3>', "\n")
+		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-success button-small pure-button">zipfile</a> with all output (', zippedSize, ')</p>', "\n", sep='')
+		#cat( '<a href="#" class="button-success button-xsmall pure-button">test</a>' )
+		
+		
+		## ------------------------
+		## table with general info
+		## ------------------------
+		cat( '<h3 class="qdnaseq">Settings</h3><p>', "\n")
+		cat( '<table class="pure-table pure-table-striped"><thead><tr><th>setting</th><th>value</th></tr></thead><tbody>' )
+			cat( htmlTableRow( c( "AnalysisName", outputName ) ) )
+			cat( htmlTableRow( c( "AnalysisDate", 'todo' ) ) )
+			cat( htmlTableRow( c( "BinSize (kb)", binSize ) ) )
+			#cat( htmlTableRow( c( "undoSD", undoSD ) ) )
+			#cat( htmlTableRow( c( "useBlacklist", filterBlacklistedBins ) ) )
+			
+			for ( galaxyName in fileList[[ 'all_files' ]] ){
+				sampleName <- fileList[[ galaxyName ]]
+				cat( htmlTableRow( c( "InputBam", paste( galaxyName, ' (', sampleName, ')', sep='' ) ) ) )
+			}
+		cat( "</tbody></table></p>", "\n")
+		
+		## ------------------------
+		## list with links to all output files
+		## ------------------------
+		r_code <- '<p>'
+		r_code <- paste( r_code, '<code class="code">## R code to load files</code><br />', sep="\n" )
+		r_code <- paste( r_code, '<code class="code">library( QDNAseq )</code><br />', sep="\n")
+
+		cat( '<h3 class="qdnaseq">Output files</h3><p>', "\n")
+		cat( '<dl>', "\n" )
+            #cat( '<dt>Definition term</dt>', "\n", '<dd>Definition term</dd>', "\n" )
+            cat( '<dt>', htmlLink( path=robjReadCoName, robjReadCoName ), '</dt>', "\n" ) 
+            cat( '<dd>QDNAseq object with read counts per bin, ', readCoSize,'</dd>', "\n" )
+            r_code <- paste( r_code, '<code class="code">readCounts <- readRDS(', robjReadCoName, ")</code><br />", sep="")
+
+            cat( '<dt>', htmlLink( path=robjCopyNrName, robjCopyNrName ), '</dt>', "\n" ) 
+            cat( '<dd>QDNAseq object with copy numbers per bin, ', copyNrSize,'</dd>', "\n" )
+            r_code <- paste( r_code, '<code class="code">copyNumbersSmooth <- readRDS(', robjCopyNrName, ")</code><br />", sep="")
+            
+            if ( doCall ){
+            	cat( '<dt>', htmlLink( path=robjCalledName, robjCalledName ), '</dt>', "\n" ) 
+            	cat( '<dd>QDNAseq object with segment and call status per bin, ', calledSize,'</dd>', "\n" )
+            	r_code <- paste( r_code, '<code class="code">copyNumbersCalled <- readRDS(', robjCalledName, ")</code><br />", sep="")
+
+            	cat( '<dt>', htmlLink( path=regiOutputName, regiOutputName ), '</dt>', "\n" ) 
+            	cat( '<dd>list with segmented/called regions for each sample, ', regionSize, '</dd>', "\n" )
+            	r_code <- paste( r_code, '<code class="code">calledRegions <- readRDS(', regiOutputName, ")</code><br />", sep="")
+
+            }
+        cat( '</dl></p>', "\n" )
+
+        cat( r_code, "</p>\n", sep="\n")
+        cat( '<p>See ', htmlLink( WEB_LINK, 'the bioconductor QDNAseq documentation' ), ' for more information on how to work with these files</p>', "\n", sep='' )
+
+		## ------------------------
+		## table with links to files	
+		## ------------------------
+		cat( '<h3 class="qdnaseq">Results: overview</h3><p>', "\n")
+		plots_html <- ''
+		
+		cat( '<table class="pure-table pure-table-striped"><thead><tr><th>Sample</th><th>CopyNumber</th><th>Called</th><th>ReadCount</th><th>RegionCount</th><th>Files</th></tr></thead><tbody>' )
+			
+			for ( bam_file in sampleNames ){
+				
+				#width <- 600; height <- 240
+				width <- PLOT_WIDTH; height <- PLOT_HEIGHT
+				width_t <- 100; height_t <- 40
+
+				## add thumbnails to table with links to anchors on html page
+				copy_img <- plotted_images[[ bam_file ]][[ 'CopyNumbers' ]]
+				usedReads <- plotted_images[[ bam_file ]][[ 'usedReads' ]]
+				usedReads <- format( as.integer(usedReads), digits=4, decimal.mark=".", big.mark="," )
+
+				html_copy_thumb <- htmlLink( path=paste('#', copy_img, sep=''), paste('<img src="',copy_img,'" alt="', bam_file, '" width="', width_t, '" height="', height_t, '">', sep='') )
+				html_copy_img <- htmlLink( path=copy_img, paste('<img id="', copy_img,'" src="',copy_img,'" alt="',bam_file, '" width="', width, '" height="', height, '">', sep='') )
+				html_call_thumb <- 'NA'
+				html_call_img <- ''
+				html_bed <- 'NA'
+				html_bedGraph <- 'NA'
+				region_count <- 'NA'
+
+				if ( doCall ){
+					call_img <- plotted_images[[ bam_file ]][[ 'Called' ]]
+					region_count <- plotted_images[[ bam_file ]][[ 'region_count' ]]
+					html_call_thumb <- htmlLink( path=paste('#', call_img, sep=''), paste('<img src="', call_img, '" alt="', bam_file, '" width="', width_t,'" height="', height_t,'">', sep='') )
+
+					files <- printedFiles[[ bam_file ]]
+					html_bed <- htmlLink( path=files[1], 'bed' )	
+					html_bedGraph <- htmlLink( path=files[2], 'bedGraph' )	
+					html_call_img <- htmlLink( path=call_img, paste('<img id="', call_img,'" src="', call_img,'" alt="', bam_file, '" width="', width, '" height="', height,'">', sep='') )
+				}
+
+				## add info to overview table, including small thumbnails
+				cat( htmlTableRow( c(bam_file, html_copy_thumb, html_call_thumb, usedReads, region_count, paste( html_bed, html_bedGraph, sep=", ")) ), "\n" )
+				## now include (large) images in html page
+				plots_html <- paste( plots_html, html_copy_img, "\n", html_call_img, "\n<hr \\>\n", sep='' )
+			}
+		cat( "</tbody></table></p>", "\n")
+
+		## add noise plot
+		html_noise_img <- htmlLink( path=noise_img_file, paste('<img id="', noise_img_file,'" src="',noise_img_file,'" alt="NoisePlot">', sep='') )
+		plots_html <- paste( plots_html, html_noise_img, "\n<hr \\>\n", sep='' )
+		
+		## ------------------------
+		## section with various output shown
+		## ------------------------
+		cat( '<h3 class="qdnaseq">Results: plots</h3><p>', "\n")
+		cat( plots_html, "\n")
+		cat( "\n</p></body>\n")
+		cat( "\n</html>\n")
+sink()
+
+## ------------------------
+## creating html output to be viewed in middle galaxy pane
+## ------------------------
+#sink( file = htmlOutputPath, type = "output" )
+sink( file = outputHtml, type = "output" )
+		
+		cat( "<head>", "\n")
+			cat( "\t", '<link rel="stylesheet" href="http://yui.yahooapis.com/pure/0.4.2/pure-min.css">', "\n" )
+
+			cat( "<style>", "\n")
+				## have to include CSS into html file, because css referencing outside own dir doesn't seem to work...makes it more portable anyway :P
+				cat( paste( "\t", '/* the css here originates from ', CSS_FILE,' */', "\n") )
+				cat( paste( "\t", readLines( CSS_FILE, n = -1)), sep="\n" )
+			cat( "</style>", "\n")
+		cat( "</head>", "\n")
+
+
+		cat( '<h1>QDNAseq Results (', outputName,')</h1>', "\n", sep="")
+		cat( '<p>Explore <a href="', htmlOutputName, '" class="button-success button-small pure-button">the results</a> directly within galaxy</p>', "\n", sep="")
+		cat( '<p>Or download a <a href="', gzipOutputName, '" class="button-success button-small pure-button">zipfile</a> with all output (', zippedSize, ')</p>', "\n", sep="" )
+
+		cat( '<p>The zip file contains all output files, including *.rds files allowing you to load the R copyNumber object(s) and perform further detailed analysis or create your own output for further processing. You can load the rds file with <code class="code">loadRDS(file.rds)</code></p>', "\n", sep="")
+		
+sink()
+
+
+cat( MAIN_NAME, "...zipping output\n")
+zip_cmd <- paste( "zip -j ", gzipOutputPath, paste(outputPath,'/', htmlOutputName, sep='') ) ## -j is for removing dirs from the tree
+system( zip_cmd )
+cat( MAIN_NAME, "done...\n", sep="")
+q(status=0)
\ No newline at end of file