Mercurial > repos > sblanck > mpagenomics
view extractCN.R @ 6:7076911e5c64 draft default tip
"planemo upload for repository https://github.com/sblanck/MPAgenomics4Galaxy/tree/master/mpagenomics_wrappers commit a644ed69951bcc1ac46426c5e6c9a0af1003a9a8-dirty"
author | sblanck |
---|---|
date | Tue, 20 Apr 2021 15:00:42 +0000 |
parents | 3fcbb8030fcc |
children |
line wrap: on
line source
#!/usr/bin/env Rscript # setup R error handling to go to stderr options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) # we need that to not crash galaxy with an UTF8 error on German LC settings. # loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") library("optparse") library("zip") ##### Read options option_list=list( make_option("--chrom",type="character",default=NULL, dest="chrom"), make_option("--input",type="character",default=NULL, dest="input"), make_option("--zip",type="character",default=NULL, dest="zip"), make_option("--output",type="character",default=NULL, dest="output"), make_option("--new_file_path",type="character",default=NULL, dest="new_file_path"), make_option("--settings_type",type="character",default=NULL, dest="settings_type"), make_option("--settings_tumor",type="character",default=NULL, dest="settings_tumor"), make_option("--symmetrize",type="character",default=NULL, dest="symmetrize"), make_option("--settings_signal",type="character",default=NULL, dest="settings_signal"), make_option("--settings_snp",type="character",default=NULL, dest="settings_snp"), make_option("--outputlog",type="character",default=NULL, dest="outputlog"), make_option("--log",type="character",default=NULL, dest="log"), make_option("--userid",type="character",default=NULL, dest="userid") ); opt_parser = OptionParser(option_list=option_list); opt = parse_args(opt_parser); if(is.null(opt$input)){ print_help(opt_parser) stop("input required.", call.=FALSE) } #loading libraries chrom=opt$chrom input=opt$input zip=opt$zip tmp_dir=opt$new_file_path output=opt$output settingsType=opt$settings_type tumorcsv=opt$settings_tumor symmetrize=opt$symmetrize signal=opt$settings_signal snp=type.convert(opt$settings_snp) outputlog=opt$outputlog log=opt$log user=opt$userid library(MPAgenomics) library(aroma.affymetrix) library(R.utils) #workdir=file.path(tmp_dir, "mpagenomics",user) tmp_dir tmp_dir=file.path(tmp_dir) if (!dir.exists(tmp_dir)) dir.create(tmp_dir, showWarnings = TRUE, recursive = TRUE) setwd(tmp_dir) # tmpzip=file.copy(from = zip,to=paste0(workdir,"/tmp.zip")) # tmpzip unzip(zipfile = zip,exdir = ".") # if (file.exists(tmpzip)) { # #Delete file if it exists # file.remove(fn) # } workdir=file.path(tmp_dir,user) setwd(workdir) inputDataset=read.table(file=input,stringsAsFactors=FALSE) dataset=inputDataset[1,2] if (outputlog){ sinklog <- file(log, open = "wt") sink(sinklog ,type = "output") sink(sinklog, type = "message") } if (grepl("all",tolower(chrom)) | chrom=="None") { chrom_vec=c(1:25) } else { chrom_tmp <- strsplit(chrom,",") chrom_vecstring <-unlist(chrom_tmp) chrom_vec <- as.numeric(chrom_vecstring) } if (signal == "CN") { if (settingsType == "dataset") { if (tumorcsv== "None") { CN=getCopyNumberSignal(dataset,chromosome=chrom_vec, onlySNP=snp) } else { CN=getCopyNumberSignal(dataset,chromosome=chrom_vec, normalTumorArray=tumorcsv, onlySNP=snp) } } else { input_tmp <- strsplit(settingsType,",") input_tmp_vecstring <-unlist(input_tmp) input_vecstring = sub("^([^.]*).*", "\\1", input_tmp_vecstring) if (tumorcsv== "None") { CN=getCopyNumberSignal(dataset,chromosome=chrom_vec, listOfFiles=input_vecstring, onlySNP=snp) } else { CN=getCopyNumberSignal(dataset,chromosome=chrom_vec, normalTumorArray=tumorcsv, listOfFiles=input_vecstring, onlySNP=snp ) } } list_chr=names(CN) CN_global=data.frame(check.names = FALSE) for (i in list_chr) { chr_data=data.frame(CN[[i]],check.names = FALSE) CN_global=rbind(CN_global,chr_data) } names(CN_global)[names(CN_global)=="featureNames"] <- "probeName" write.table(format(CN_global), output, row.names = FALSE, quote = FALSE, sep = "\t") } else { if (symmetrize=="TRUE") { if (settingsType == "dataset") { input_vecstring = getListOfFiles(dataset) } else { input_tmp <- strsplit(settingsType,",") input_tmp_vecstring <-unlist(input_tmp) input_vecstring = sub("^([^.]*).*", "\\1", input_tmp_vecstring) } symFracB_global=data.frame(check.names = FALSE) tumorFile=read.csv(tumorcsv,header=TRUE) tumor=tumorFile$tumor input_vecstring=input_vecstring[which(input_vecstring %in% tumor)] for (currentFile in input_vecstring) { cat(paste0("extracting signal from ",currentFile,".\n")) currentSymFracB=data.frame() symFracB=getSymFracBSignal(dataset,chromosome=chrom_vec,file=currentFile,normalTumorArray=tumorcsv) list_chr=names(symFracB) for (i in list_chr) { cat(paste0(" extracting ",i,".\n")) chr_data=data.frame(symFracB[[i]]$tumor,check.names = FALSE) currentSymFracB=rbind(currentSymFracB,chr_data) } if (is.null(symFracB_global) || nrow(symFracB_global)==0) { symFracB_global=currentSymFracB } else { #symFracB_global=cbind(symFracB_global,currentFile=currentSymFracB[[3]]) symFracB_global=merge(symFracB_global,currentSymFracB[,c(3,4)],by="featureNames") symFracB_global=symFracB_global[c(2:ncol(symFracB_global),1)] symFracB_global=symFracB_global[order(symFracB_global$chromosome, symFracB_global$position),] } } names(symFracB_global)[names(symFracB_global)=="featureNames"] <- "probeName" write.table(format(symFracB_global), output, row.names = FALSE, quote = FALSE, sep = "\t") } else { if (settingsType == "dataset") { if (tumorcsv== "None") { fracB=getFracBSignal(dataset,chromosome=chrom_vec) } else { fracB=getFracBSignal(dataset,chromosome=chrom_vec, normalTumorArray=tumorcsv) } } else { input_tmp <- strsplit(settingsType,",") input_tmp_vecstring <-unlist(input_tmp) input_vecstring = sub("^([^.]*).*", "\\1", input_tmp_vecstring) if (tumorcsv== "None") { fracB=getFracBSignal(dataset,chromosome=chrom_vec, listOfFiles=input_vecstring) } else { fracB=getFracBSignal(dataset,chromosome=chrom_vec, normalTumorArray=tumorcsv, listOfFiles=input_vecstring) } } #formatage des données list_chr=names(fracB) fracB_global=data.frame(check.names = FALSE) for (i in list_chr) { chr_data=data.frame(fracB[[i]]$tumor,check.names = FALSE) fracB_global=rbind(fracB_global,chr_data) } names(fracB_global)[names(fracB_global)=="featureNames"] <- "probeName" write.table(format(fracB_global), output, row.names = FALSE, quote = FALSE, sep = "\t") } } if (dir.exists(workdir)) system(paste0("rm -r ", workdir)) if (outputlog){ sink(type="output") sink(type="message") close(sinklog) }