configure_dnaseq.R

cwd = getwd()

species = 'B37' # genome
wesImpactTarget = 'AgilentExon_51MB_b37_v3' # baits
assay = 'wes'

delFiles = F # delete files after each step 

if(!file.exists(groupfile)){
	stop('groupfile not exist.')
}

if(!file.exists(mapfile)){
	stop('mapfile not exist.')
}

if(!exists('species')){
	cat('Species must provided\n\n')
	cat('Environment parameters are not imported!!\n\n')
	return
}

if(!exists('wesImpactTarget')){
	cat('wesImpactTarget must provided\n\n')
	cat('Environment parameters are not imported!!\n\n')
	return
}

statsDir = paste0(resDir, '/stats')
initDir = paste0(resDir, '/initFiles')
progressDir = paste0(resDir, '/progress')
matricsDir = paste0(resDir, '/matrics')
varDir = paste0(resDir, '/variation')
haploDir = paste0(varDir, '/haplotypecaller')
mutectDir = paste0(varDir, '/mutect')
sniperDir = paste0(varDir, '/somaticsniper')
facetsDir = paste0(varDir, '/facets')
strvarDir = paste0(varDir, '/strvar')
system(paste0('mkdir -p ', resDir)) ## r_fang
system(paste0('mkdir -p ', statsDir)) ## err std files
system(paste0('mkdir -p ', progressDir)) ## jobname.done
system(paste0('mkdir -p ', matricsDir))  ## collected metrics
system(paste0('mkdir -p ', varDir))  
system(paste0('mkdir -p ', haploDir)) 
system(paste0('mkdir -p ', mutectDir)) 
system(paste0('mkdir -p ', sniperDir)) 
system(paste0('mkdir -p ', facetsDir)) 
system(paste0('mkdir -p ', facetsDir, '/sample')) 
system(paste0('mkdir -p ', strvarDir)) 
system(paste0('mkdir -p ', initDir)) 

##logfile = paste0('_', resDir, '_log')
##if(file.exists(logfile)){
##	stop(paste0(logfile, 'exist, already running?'))
##}else{
##	create.logger() -> logger
##	logfile(logger) = logfile
##	level(logger) = 'INFO' 
##	#info(logger, 'message') others like warn(), error(), fatal()
##}


## this configure for the locations for common programs
## the task.conf under the working directory provide the tasks. 
## save the input and base as special words here
## _genome will sub with $genome
#>Locations
source(paste0(resDir, '/configure_dir.R'))

## this for set up @ENV
path			= '/home/huw/program/homer/bin:/home/huw/local/bin:/home/huw/perl5/CPAN/bin:/home/huw/program/cufflink221/:/home/huw/program/tophat213/:/home/huw/program/bowtie2/:/home/huw/program/bin:/home/huw/program/blat:/home/huw/program/weblogo:/home/huw/program/gs/bin:/home/huw/program/ngsplot:$PATH'
JAVA_HOME 		= '/home/huw/program/jdk7/bin/java'
JAVA_HOME 		= '/opt/common/CentOS_6-dev/java/jdk1.8.0_31/'
NGSPLOT 		= '/home/huw/program/ngsplot'
TMPDIR			= '/scratch/huw'

if(species == 'mm9') {
	genomeFasta = MM9_FASTA
	genomeBWA =  MM9_BWA_INDEX
	genomeFAI = MM9_FAI
	DB_SNP = "" 
	FP_INT = "" 
	FP_TG  = ""
}

if(species == 'mm10_custom') {
	genomeFasta = MM10_CUSTOM_FASTA
	genomeBWA =  MM10_BWA_INDEX
	DB_SNP = paste0(dataDir, "/mm10/mm10_snp142.vcf") 
	FP_INT = "" 
	FP_TG  = ""
}

if(species == 'mm10') {
	genomeFasta = MM10_FASTA
	genomeBWA =  MM10_BWA_INDEX
	DB_SNP = paste0(dataDir, "/mm10/mm10_snp142.vcf") 
	FP_INT = "" 
	FP_TG  = ""
}

if(species == 'hybrid') {
	genomeFasta = B37_MM10_HYBRID_FASTA
	#genomeBWA =  B37_BWA_INDEX
	DB_SNP = paste0(dataDir, "/b37/dbsnp_138.b37.vcf")
	FP_INT = paste0(dataDir, "/b37/Agilent51MBExome__b37__FP_intervals.list")
	FP_TG  = paste0(dataDir, "/b37/Agilent51MBExome__b37__FP_tiling_genotypes.txt")
}

baits_ilist = paste0(targetsDir, '/', wesImpactTarget, '/', wesImpactTarget, '_baits.ilist')
targets_ilist = paste0(targetsDir, '/', wesImpactTarget, '/', wesImpactTarget, '_targets.ilist')
targets_bed = paste0(targetsDir, '/', wesImpactTarget, '/', wesImpactTarget, '_targets.bed')
targets5bp_ilist = paste0(targetsDir, '/', wesImpactTarget, '/', wesImpactTarget, '_targets_plus5bp.ilist')
targets5bp_bed = paste0(targetsDir, '/', wesImpactTarget, '/', wesImpactTarget, '_targets_plus5bp.bed')

if(grepl('hg19', species, ignore.case=T)){
	   genomeFasta = HG19_FASTA
	   genomeBWA = HG19_BWA_INDEX
	   genomeFAI = HG19_FAI
	   DB_SNP = paste0(dataDir, "/hg19/dbsnp_138.hg19.vcf")
	   FP_INT = paste0(dataDir, "/hg19/Agilent51MBExome__hg19__FP_intervals.list")
	   FP_TG  = paste0(dataDir, "/hg19/Agilent51MBExome__hg19__FP_tiling_genotypes.txt")
	   ExAC_VCF = paste0(VEP, "/ExAC_nonTCGA.r0.3.1.sites.vep.vcf.gz");
	   FACETS_DB_SNP = paste0(dataDir, "/hg19/dbsnp_137.hg19__RmDupsClean__plusPseudo50__DROP_SORT.vcf");
	   MILLS_1000G = paste0(dataDir, "/hg19/Mills_and_1000G_gold_standard.indels.hg19.vcf");
	   HAPMAP = paste0(dataDir, "/hg19/hapmap_3.3.hg19.vcf");
	   OMNI_1000G = paste0(dataDir, "/hg19/1000G_omni2.5.hg19.vcf");
	   PHASE1_SNPS_1000G = paste0(dataDir, "/hg19/1000G_phase1.snps.high_confidence.hg19.vcf");
	   COSMIC = paste0(dataDir, "/hg19/CosmicCodingMuts_v67_20131024.vcf");
}

if(grep('b37', species, ignore.case=T)){
	genomeFasta = B37_FASTA
	genomeBWA =  B37_BWA_INDEX
	genomeFAI = B37_FAI
	DB_SNP = paste0(dataDir, "/b37/dbsnp_138.b37.vcf")
	FP_INT = paste0(dataDir, "/b37/Agilent51MBExome__b37__FP_intervals.list")
	FP_TG  = paste0(dataDir, "/b37/Agilent51MBExome__b37__FP_tiling_genotypes.txt")
	ExAC_VCF = paste0(VEP, "/ExAC_nonTCGA.r0.3.1.sites.vep.vcf.gz");
	FACETS_DB_SNP = paste0(dataDir, "/b37/dbsnp_137.b37__RmDupsClean__plusPseudo50__DROP_SORT.vcf");
	MILLS_1000G = paste0(dataDir, "/b37/Mills_and_1000G_gold_standard.indels.b37.vcf");
	HAPMAP = paste0(dataDir, "/b37/hapmap_3.3.b37.vcf");
	OMNI_1000G = paste0(dataDir, "/b37/1000G_omni2.5.b37.vcf");
	PHASE1_SNPS_1000G = paste0(dataDir, "/b37/1000G_phase1.snps.high_confidence.b37.vcf");
	COSMIC = paste0(dataDir, "/b37/CosmicCodingMuts_v67_b37_20131024__NDS.vcf");
	COSMIC_HOTSPOTS = paste0(dataDir, "/b37/dmp_cosmic_for_hotspots.vcf");
	vep.species = 'homo_sapiens'
}

exe.jobs = function(jobnames=NULL, logger=logger, outputs=NULL, size=2000){
	if(is.null(jobnames) & !is.vector(jobnames)){
		msg = "jobnames parameter must vector of character"
		fatal(logger, msg)
		stop(msg)
	}
	for(i in 1:length(jobnames)){
		if(is.null(outputs)){ # run job without check
			jobname = jobnames[i]
			info(logger, jobname)
			ret = system(jobname)
			info(logger, paste0('return value for the above cmd: ', ret))
		}else if(!is.null(outputs[i]) & file.exists(outputs[i]) & file.size(output[i]) > size){ # no run if output exist
			msg = paste0('Output file ', outputs[i], ' exist, job not submitted for this')
			info(logger, msg)
		}else{ # either output file not exists or is too small
			jobname = jobnames[i]
			info(logger, jobname)
			ret = system(jobname)
			info(logger, paste0('return value for the above cmd: ', ret))
		}
	}
}

deldir = function(ddir=NULL){
	if(!is.null(ddir) & dir.exists(ddir)){
		msg = paste0('deleting folder ', ddir)
		warn(logger, msg)
		system(paste0('rm -rf ', ddir))
		return(0)
	}
	return(1)
}

delfile.2 = function(dfile=NULL, logger){
	if(!is.null(dfile) & file.exists(dfile)){
		msg = paste0('deleting file ', dfile)
		warn(logger, msg)
		system(paste0('rm -f ', dfile))
		return(0)
	}
	return(1)
}

delfile = function(dfile=NULL, logger, cfile=NULL, s){
	if(!is.null(dfile) & !file.exists(dfile)){
		msg = paste0(dfile, ' not exists')
		info(logger, msg);stop(msg)
	}
	if(!is.null(cfile) & !file.exists(cfile)){
		msg = paste0(cfile, ' not exists')
		info(logger, msg);stop(msg)
	}else if(file.size(cfile) < s){
		msg = paste0(cfile, ' is too small')
		info(logger, msg);stop(msg)
	}else{
		msg = paste0('deleting file ', dfile)
		warn(logger, msg)
		system(paste0('rm -f ', dfile))
		return(0)
	}
	return(1)
}

bjoblist = function(){
	tmpfile = tempfile()
	system(paste0('bjobs -w > ', tmpfile), wait=T)
	if(file.exists(tmpfile) & file.size(tmpfile) > 5 ) {
		joblist = suppressWarnings(fread(tmpfile))
		system(paste0('rm -f ', tmpfile))
		return(joblist$V7)
	}else{
		system(paste0('rm -f ', tmpfile))
		return('')
	}
}

wait4jobs = function(jobnames, logger){
	tmpfile = tempfile()
	system(paste0('bjobs -w > ', tmpfile), wait=T)
	joblist = 1
	if(file.exists(tmpfile) & file.size(tmpfile) > 5 ) joblist = fread(tmpfile)
	while(joblist != 1 & length(intersect(joblist$V7, jobnames)) > 0){
		info(logger, 'waiting for the below jobs to finish')
		info(logger, paste(intersect(joblist$V7, jobnames), collapse = ' '))
		Sys.sleep(100);
	}
	system(paste0('rm -f ', tmpfile))
}

#>Parameters
isPE			= '1'
clusters		= 'LSF'
#>Versions

#>Convention
##name1fastq		= 'base_xx_R1.fastq.gz'
##name2fastq		= 'base_xx_R2.fastq.gz'
##name1mrgedfq		= 'base_R1.fq.gz'
##name2mrgedfq		= 'base_R2.fq.gz'
##name1trimmedfq		= 'base_ptrimmed_R1.fq.gz'
##name2trimmedfq		= 'base_ptrimmed_R2.fq.gz'
##sai1			= 'base_R1.sai'
##sai2			= 'base_R2.sai'
##sam			= 'base_genome.sam'
##bam			= 'base_genome'
##bamext			= 'base_genome.bam'
##bamsorted		= 'base_genome_sorted'
##bamsortedext		= 'base_genome_sorted.bam'
##bamsortedrmdup		= 'base_genome_sorted_rmdup.bam'
##htseqcount		= 'base.htseq.count'
##bigwig			= 'base_genome_sorted_rmdup_10m.bw'
##MACS14_peaks_file 	= 'base_genome_macs14_peaks.bed'
##MACS2_peaks_file	= 'base_genome_macs2_peaks.xls'
##MACS2_broad_peaks_file  = 'base_genome_macs2broad_peaks.broadPeak'
##peakannoout		= 'base_genome_macs2_peaks_anno.txt'
##mcspeaks		= 'base_genome_macs2_peaks.narrowPeak'
##mcsbroadpeaks		= 'base_genome_macs2broad_peaks.broadPeak'
##peakannout		= 'base_genome_macs2_peaks_anno.txt'
##peakannout2		= 'base_genome_macs2broad_peaks_anno.txt'
##peakannoFiltered	= 'base_genome_macs2_peaks_anno_filtered.txt'
##peakannoFiltered2	= 'base_genome_macs2broad_peaks_anno_filtered.txt'
##ngsplot_file		= 'ngsplot_base.tss.avgprof.pdf'
##ngsplot_file_heatmap	= 'ngsplot_base.tss.heatmap.pdf'
##TSSCounts		= 'base_genome_tss2kb.bed'
##tssout			= 'ngsplot.base.tss'
##genebodyout		= 'ngsplot.base.genebody'
##enhancerout		= 'ngsplot.base.enhancer'
##motifout		= 'motif'
##starbam        		= 'base.Aligned.sortedByCoord.out.bam'
##
##>Checkpoints
##mergefastq		= 'name1mrgedfq'
##trimfq			= 'name1trimmedfq'
##alignStarAlign		= 'bamsortedext'
##alignBowtie2		= 'bamsrotedext'
##alignBwaAln		= 'sai1'
##alignBwaSampe		= 'sam'
##sortSam			= 'bamsortedext'
##sam2bam			= 'bamsortedext'
##removeDup		= 'bamsortedrmdup'
##bigwig			= 'bigwig'
##indexBam		= 'bamInput_bai'
##MACS2			= 'mcspeaks'
##MACS2BroadPeak		= 'mcsbroadpeaks'
##annotatePeaksHomer	= 'peakannoout_peakannoout2'
##annotateFilter		= 'peakannoFiltered_peakannoFiltered2'
##macs14			= 'mcspeaks'
##TSSCounts		= 'TSSCounts'
##Ngsplot			= 'tssout.avgprof.pdf'
##motifHomerFinder	= 'motifout'
##htseqcount		= 'htseqcount'

bsub.head2 = function(jobname, mem, cpu, We, cwd, postdone, submit, resDir = 'r_fang'){
	cmd = paste0(BSUB, " -J ", jobname, " -e ", statsDir, '/', jobname, ".err -o ", statsDir, '/', jobname, ".std ")
	cmd = paste0(cmd, " -cwd ", cwd, ' -We ', We, ' -R "rusage[mem=', mem, ']" -R "rusage[iounits=0]" -n ', cpu)
	if(!missing(postdone)){
		postdone = unlist(strsplit(postdone, split=" "))
		for( i in 1:length(postdone)){
			cmd = paste0(cmd, ' -w "post_done(', postdone[i], ')" ')
		}
	}
	if(!missing(submit)){
		submit2 = unlist(strsplit(submit, " "))
		submit3 = paste0('ls ', cwd, '/progress/', resDir, '/', submit2)
		submit4 = paste0(submit3, collapse = ' && ')
		cmd = paste0(cmd, ' -E "', submit4, '"')
	}
	cmd
}

bsub.head = function(jobname, mem, cpu, We, cwd, postdone, submit, resDir = 'r_fang'){
	cmd = paste0(BSUB, " -J ", jobname, " -e ", statsDir, '/', jobname, ".err -o ", statsDir, '/', jobname, ".std ")
	cmd = paste0(cmd, " -cwd ", cwd, ' -We ', We, ' -R "rusage[mem=', mem, ']" -R "rusage[iounits=0]" -n ', cpu)
	if(!missing(postdone)){
		postdone = unlist(strsplit(postdone, split=" "))
		runningjobs = bjoblist()
		for( i in 1:length(postdone)){
			if(length(intersect(runningjobs, postdone[i])) > 0){
				cmd = paste0(cmd, ' -w "post_done(', postdone[i], ')" ')
			}
		}
	}
	if(!missing(submit)){
		submit2 = unlist(strsplit(submit, " "))
		submit3 = paste0('ls ', cwd, '/progress/', resDir, '/', submit2)
		submit4 = paste0(submit3, collapse = ' && ')
		cmd = paste0(cmd, ' -E "', submit4, '"')
	}
	cmd
}

checkfiles = function(files, sizes){
	if(!(all(file.exists(files))) & !(all(file.size(files) > sizes))){
		msg = paste0('files not passed check: ', paste(files, collapse=' '))
		fatal(logger, msg)
		stop(msg)
	}
}