qdnaseq: QDNAseq.R comparison

comparison QDNAseq.R @ 28:40ae95ad9d8c draft

Uploaded

author	stef
date	Wed, 18 Jun 2014 05:29:59 -0400
parents	2efda0ed7ccf
children	4c4e4d88779c

comparison

equal deleted inserted replaced

-:47c3ecce5544
+:40ae95ad9d8c
 ## --------------------
 catMsg <- function( msg=c() ){
 	cat( MAIN_NAME, paste( msg, collapse="" ), "\n", sep='')
 }
 ## --------------------
-## return the directory this script exists
+## return the location of this script
 ## --------------------
 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
+## Some html creation functions
 ## --------------------
 htmlTableRow <- function( string_array=c() ){
 	td_cells <- ''
 	for ( i in string_array ){
 		td_cells <- paste( td_cells, '<td>', i, '</td>', sep='' )
 }
 htmlLink <- function( path, desc="LINK" ){
 	return( paste( '<a href="', path, '">', desc, "</a>", sep='') )
 }
 ## --------------------
-## constructs a list with bam file info
+## constructs a list with input bam file info
 ## --------------------
 makeBamFileList <- function( paths, names ){
 	tmp <- list()
 	l1 <- length(paths)
 	l2 <- length(names)
 	}
 	return( tmp )
 }
 ## --------------------
-## took a function from Matias/MarkVdWuel? for extracting the regions by call
+## copied code for extracting the regions by segment call status
 ## --------------------
 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)
 		outFiles[[sample]] <- c( outBedFile, outBedgraphFile )
 	}
 	outFiles
 }
-printIgvFile <- function( dat, filename ){
+#printIgvFile <- function( dat, filename ){
+#
-	if ( 'calls' %in% assayDataElementNames(dat) ) {
+#	if ( 'calls' %in% assayDataElementNames(dat) ) {
-	 	#output <- paste(output, '-called.igv', sep="")
+#	 	#output <- paste(output, '-called.igv', sep="")
-	  	cat('#type=COPY_NUMBER\n#track coords=1\n', file=filename)
+#	  	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)
+#	  	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 {
+#	}else {
-	  #output <- paste(output, '-normalized.igv', sep="")
+#	  #output <- paste(output, '-normalized.igv', sep="")
-	  cat('#type=COPY_NUMBER\n#track coords=1\n', file=filename)
+#	  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 <- 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 == '23'] <- 'X'
-	df$chromosome[df$chromosome == '24'] <- 'Y'
+#	df$chromosome[df$chromosome == '24'] <- 'Y'
-	df$chromosome[df$chromosome == '25'] <- 'MT'
+#	df$chromosome[df$chromosome == '25'] <- 'MT'
-	#return( df )
+#	#return( df )
-	write.table( df, file=filename, append=TRUE, quote=FALSE, sep='\t', row.names=FALSE )
+#	write.table( df, file=filename, append=TRUE, quote=FALSE, sep='\t', row.names=FALSE )
-}
+#}
 ## ==================================================
 ## Start of analysis
 ## ==================================================
 #catMsg( instal_path )
 ## /DEBUG
 ## ------------------------
 ## setup bam filelist for easy retrieval later
-catMsg( "Setting up input bam list" )
+#catMsg( "Setting up input bam list" )
-cat( bamsPaths, "\n" )
+#cat( bamsPaths, "\n" )
-catMsg( "Namews" )
+#catMsg( "Namews" )
-cat( bamsNames, "\n" )
+#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"){
+#if ( length(cmdLineArgs) == 0 || cmdLineArgs[1] == "-h" || cmdLineArgs[1] == "--help"){
-	cat( paste( MAIN_NAME, "Usage: ", params_help, sep=''), "\n" )
+#	cat( paste( MAIN_NAME, "Usage: ", params_help, sep=''), "\n" )
-	quit( save = 'no', status=0, runLast=F )
+#	quit( save = 'no', status=0, runLast=F )
-}
+#}
 if ( !file.exists( outputPath) ){
 	dir.create( outputPath )
 }
 readCountsFiltered    <- applyFilters( readCounts, residual=TRUE, blacklist=filterBlacklistedBins, mappability=mappabilityCutoff, chromosomes=excludeChrs )
 readCountsFiltered    <- estimateCorrection( readCountsFiltered )
 copyNumbers           <- correctBins( readCountsFiltered )
 copyNumbersNormalized <- normalizeBins( copyNumbers )
 copyNumbersSmooth     <- smoothOutlierBins( copyNumbersNormalized )
+sampleNames           <- readCountsFiltered@phenoData@data$name
 ## save objects to output dir
 saveRDS( readCounts, robjReadCoPath );
-saveRDS( copyNumbersSmooth, robjCopyNrPath );
+saveRDS( copyNumbersSmooth, robjCopyNrPath );s
-#printIgvFile( copyNumbersSmooth, igvCopyNrName )
+exportBins(copyNumbersSmooth, file=igvCopyNrName, format="igv")
-#exportBins(copyNumbersSmooth, file=igvCopyNrName, format="igv")
 ## also save objects for galaxy history output if requested
 if ( doOutputReadcountsRds ){
 	saveRDS( readCountsFiltered, readCountsDatasetFile );
 }
 	cgh <- makeCgh( copyNumbersCalled )
 	saveRDS( copyNumbersCalled, robjCalledPath );
 	if ( doOutputCallsRds ){
 		saveRDS( copyNumbersCalled, calledSegmentsDatasetFile );
 	}
-	#df <- getIGVdf( copyNumbersCalled )
+	exportBins( copyNumbersCalled, file=igvCalledName, format="igv")
-	#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

Mercurial > repos > stef > qdnaseq

comparison QDNAseq.R @ 28:40ae95ad9d8c draft