code refactoring, bug fix for unary operator error

Author: leepc12

chipseq.bds

@@ -24,6 +24,8 @@ use_macs2_gapped_peak := false	help Use gapped peak instead of narrow peak (for
 
 help() // print help and exit if no parameters are given
 
+include "modules/pipeline_template.bds"
+
 include "modules/input.bds"
 
 include "modules/align_bwa.bds"
@@ -37,8 +39,6 @@ include "modules/callpeak_idr.bds"
 
 include "modules/signal.bds"
 
-include "modules/report.bds"
-
 
 
 // Important file names are stored in global variables (usually a string map string{} with a key with replicate id and peakcaller name)
@@ -96,8 +96,6 @@ void init_chipseq() {
 
 	read_conf_chipseq()
 
-	print_chipseq()
-
 	init_filetable()
 }
 
@@ -121,9 +119,7 @@ void read_conf_chipseq() {
 	use_macs2_gapped_peak = get_conf_val_bool( use_macs2_gapped_peak, 		["use_macs2_gapped_peak"] )
 
 	if ( histone ) type = "histone"
-}
-
-void print_chipseq() {
+	if ( idr_rank == "" ) idr_rank = is_TF_chipseq() ? "signal.value" : "p.value"
 
 	print( "\n\n== chipseq pipeline settings\n")
 	print( "Type of ChIP-Seq pipeline\t\t: $type\n")
@@ -282,7 +278,7 @@ void align() {
 		for ( int rep=1; rep <= get_num_rep( ctl ); rep++) {		
 			if ( !input_file_exists( ctl, rep ) ) continue
 
-			// check file size to distribute_nth nth to each nth_app
+			// check file size to distribute_nonzero nth to each nth_app
 			// determine # threads for each app related to alignment
 
 			group := get_group_name( ctl, rep )
@@ -300,7 +296,7 @@ void align() {
 		}
 	}
 
-	nth_rep := distribute_nth( nth, filesize )
+	nth_rep := distribute_nonzero( nth, filesize )
 
 	// align
 
@@ -717,13 +713,13 @@ void call_peaks() {
 
 		for ( int rep=1; rep<=get_num_rep(); rep++ ) {
 
-			nlines{"rep$rep"} = get_no_lines( tag{"rep$rep"} )
+			nlines{"rep$rep"} = get_num_lines( tag{"rep$rep"} )
 
 			//print("DEBUG: # lines rep $rep tag.: " + nlines{"rep$rep"} + ", " + tag{"rep$rep"} + "\n")
 
 			if ( input_file_exists( 1, rep ) ) { // if control exists
 
-				nlines{"ctl$rep"} = get_no_lines( tag{"ctl$rep"} )
+				nlines{"ctl$rep"} = get_num_lines( tag{"ctl$rep"} )
 				nlines_ctl.add( nlines{"ctl$rep"} )
 
 				//print("DEBUG: # lines ctl $rep tag.: " + nlines{"ctl$rep"} + ", " + tag{"ctl$rep"} + "\n")
@@ -799,7 +795,7 @@ void call_peaks() {
 
 	// distribute # of threads for each peak calling
 	int nth_pooled, nth_other
-	( nth_pooled, nth_other ) = distribute_nth( nth, [ get_num_rep(), 1, 1, 1 ] )
+	( nth_pooled, nth_other ) = distribute_nonzero( nth, [ get_num_rep(), 1, 1, 1 ] )
 
 	// create directories for peaks (spp, macs2), signal tracks (macs2)
 	spp_o_dir 	:= mkdir( "$out_dir/peak/spp") 
@@ -1096,11 +1092,8 @@ void do_idr() {
 
 	if ( !is_final_stage_idr() ) return
 	if ( is_histone_chipseq() ) return
-
 	if ( !peak_exists() ) return // if no peaks from peak caller exit
 
-	if ( idr_rank == "" ) idr_rank = is_TF_chipseq() ? "signal.value" : "p.value"
-
 	// peak caller
 	pc := get_peak_caller()
 
@@ -1308,7 +1301,7 @@ string html_chipseq_tracks() {
 		if ( idr_pr.hasKey(rep) ) {trk_types += "hammock"; trk_names += "$title peak idr (rep$rep-pr)"; trk_files += peak_to_hammock( _get_idr_peak_trk( idr_pr{rep} ) ) }
 	}
 
-	html := html_epg_browser_viz( trk_files, trk_types, trk_names )
+	html := html_epg_browser_viz( trk_files, trk_types, trk_names, species )
 
 	return html
 }

modules/align_bwa.bds

@@ -2,7 +2,7 @@
 #vim: syntax=java
 
 include "species.bds"
-include "report.bds"
+include "module_template.bds"
 
 
 help == align bwa settings (requirements: -bwa_idx)
@@ -72,7 +72,7 @@ string[] bwa_PE( string fastq1, string fastq2, string o_dir, string log_o_dir, s
 
 	if ( out <- in ) { // compare file timestamps of in and out (to check if job is already done or not)
 
-		nth_bwa_aln := distribute_nth( nth_bwa, [1,1] )
+		nth_bwa_aln := distribute_nonzero( nth_bwa, [1,1] )
 
 		// parallel jobs
 		sai1 := bwa_aln( fastq1, o_dir, group+"_1", nth_bwa_aln[0] )

modules/align_multimapping.bds

@@ -8,12 +8,10 @@ help == align multimapping settings
 multimapping 	:= 0 		help # alignments reported for multimapping (default: 0).
 
 
-
 init_align_multimapping()
 
 
 void init_align_multimapping() {
-
 	multimapping 	= get_conf_val_int( multimapping,		["multimapping"] )
 
 	print("\n\n== align multimapping settings\n")

modules/callpeak_idr.bds

@@ -2,7 +2,7 @@
 #vim: syntax=java
 
 include "species.bds"
-include "report.bds"
+include "module_template.bds"
 
 
 help == IDR settings
@@ -19,26 +19,20 @@ void init_idr() {
 
 	idr_suffix 	= get_conf_val_bool( idr_suffix, 	["idr_suffix"] )
 
-	print_idr()
-}
-
-void print_idr() {
-
 	print( "\n\n== IDR settings\n")
 	print( "Append IDR threshold to IDR out_dir\t: $idr_suffix\n" )
 }
 
 void chk_idr() {
-	if ( !path_exists( blacklist ) ) print("\nWarning: Blacklist bed (-blacklist) for final IDR QC is missing! (file: $blacklist)\n\n")
+	if ( !path_exists( blacklist ) ) print("\nWarning: Blacklist bed (-blacklist) for final IDR QC"+\
+						" is missing! (file: $blacklist)\n\n")
 }
 
-
 // 12-col.bed.gz (10 col for narrowPeak + local IDR, global IDR) will be generated for browser tracks
 // A temporary 13th column bed should be generated due to bedtools intersect compatibility issue
 // use function _get_idr_peak_trk( idr_peak ) to find 12-col.bed.gz for browser tracks
 // rank: 'signal.value' for SPP peaks, 'p.value' for MACS2 peaks
 string[] idr2( string peak1, string peak2, string pooled, string idr_thresh, string rank, string o_dir, string group ) {
-
 	if ( idr_suffix ) o_dir = mkdir( o_dir + "_$idr_thresh" )
 
 	prefix 			:= "$o_dir/" + (title=="" ? "" : (title+"_") ) + group	
@@ -111,12 +105,11 @@ string[] idr2( string peak1, string peak2, string pooled, string idr_thresh, str
 // idr_pr: IDR peaks for pseudo replicates (key:"rep") e.g. "1" for pseudo replicates of replicate 1
 string[] idr_final_qc( string{} idr_tr, string{} idr_pr, string idr_ppr, \
 			string idr_o_dir, string qc_o_dir, string group ) {
-
 	//// get # of lines in each IDR peak files
 
 	// pseudo replicates
 	int[] N
-	for ( string peak : idr_pr ) N.add( get_no_lines( peak ) )
+	for ( string peak : idr_pr ) N.add( get_num_lines( peak ) )
 
 	// true replicates
 	Nt := 0
@@ -125,7 +118,7 @@ string[] idr_final_qc( string{} idr_tr, string{} idr_pr, string idr_ppr, \
 	for ( string key : idr_tr.keys() ) { // key = rep id 1, rep id 2
 
 		peak := idr_tr{ key }
-		nlines := ( peak != "" ) ? get_no_lines( peak ) : 0
+		nlines := ( peak != "" ) ? get_num_lines( peak ) : 0
 
 		if ( nlines >= Nt ) {
 			Nt = nlines
@@ -135,7 +128,7 @@ string[] idr_final_qc( string{} idr_tr, string{} idr_pr, string idr_ppr, \
 	}
 
 	// pooled pseudo replicates
-	Np := ( idr_ppr != "" ) ? get_no_lines( idr_ppr ) : 0
+	Np := ( idr_ppr != "" ) ? get_num_lines( idr_ppr ) : 0
 
 	//// find optimal set and conservative set
 	string optimal_set, conservative_set
@@ -180,11 +173,11 @@ string[] idr_final_qc( string{} idr_tr, string{} idr_pr, string idr_ppr, \
 
 	//// compute IDR scores
 
-	real max_Np_Nt = max( Np, Nt )
-	real min_Np_Nt = min( Np, Nt )
+	real max_Np_Nt = _max( Np, Nt )
+	real min_Np_Nt = _min( Np, Nt )
 
-	real max_N = N.size() > 0 ? max( N ) : 0
-	real min_N = N.size() > 0 ? min( N ) : 0
+	real max_N = N.size() > 0 ? _max( N ) : 0
+	real min_N = N.size() > 0 ? _min( N ) : 0
 
 	rescue_ratio 		:= max_Np_Nt / min_Np_Nt
 	self_consistency_ratio 	:= max_N / min_N
@@ -227,14 +220,32 @@ string[] idr_final_qc( string{} idr_tr, string{} idr_pr, string idr_ppr, \
 }
 
 string _get_group_from_key( string key ) { //parse "i,j" to "repi-repj"
-
 	tmp := key.split(",")
 	i := tmp[0]
 	j := tmp[1]
 	return "rep$i-rep$j"
 }
 
 string _get_idr_peak_trk( string idr_peak ) {
-
 	return rm_ext( idr_peak, ["narrowPeak","regionPeak"] ) + ".12-col.bed.gz"	
 }
+
+int _max( int a, int b ) {
+	return max( a, b )
+}
+
+int _min( int a, int b ) {
+	return min( a, b )
+}
+
+int _max( int[] arr ) {
+	tmp := arr[0]
+	for( int i=1; i<arr.size();i++ ) tmp = _max( tmp, arr[i] )
+	return tmp
+}
+
+int _min( int[] arr ) {
+	tmp := arr[0]
+	for( int i=1; i<arr.size();i++ ) tmp = _min( tmp, arr[i] )
+	return tmp
+}

modules/callpeak_macs2.bds

@@ -2,7 +2,7 @@
 #vim: syntax=java
 
 include "species.bds"
-include "report.bds"
+include "module_template.bds"
 
 
 help == callpeak macs2 settings (requirements: -chrsz -gensz)

modules/callpeak_naive_overlap.bds

@@ -2,7 +2,7 @@
 #vim: syntax=java
 
 include "species.bds"
-include "report.bds"
+include "module_template.bds"
 
 
 grp_color_naive_overlap := "skyblue"

modules/callpeak_spp.bds

@@ -2,7 +2,7 @@
 #vim: syntax=java
 
 include "species.bds"
-include "report.bds"
+include "module_template.bds"
 
 
 help == callpeak spp settings

modules/cluster.bds

@@ -0,0 +1,83 @@
+#!/usr/bin/env bds
+#vim: syntax=java
+
+include "conf.bds"
+
+
+help == cluster/system/resource settings
+wt 		:= "5h50m"	help Walltime for all single-threaded tasks (example: 8:10:00, 3h, 3600, default: 5h50m, 5:50:00).
+memory 		:= "7G"		help Maximum memory for all single-threaded tasks (equivalent to '-mem', example: 4.5G, 1024M, default: 7G).
+use_system 	:= "" 		help Force to use a system (equivalent to 'bds -s [SYSTEM_NAME] ...', any system defined in bds.config can be used).
+nice 		:= 0 		help Set process priority for all tasks (default: 0; -20 (highest) ~ 19 (lowest) ).
+retrial		:= 0 		help # of Retrial for failed tasks (default: 0).
+q 		:= ""		help Submit tasks to a specified cluster queue.
+unlimited_mem_wt:= false 	help Use unlimited max. memory and walltime.
+
+
+init_cluster()
+
+
+void init_cluster() { 
+	wt 		= get_conf_val( wt, 			["wt"] )
+	memory 		= get_conf_val( memory, 		["memory","mem"] )
+	use_system 	= get_conf_val( use_system, 		["use_system","system"] )
+	nice 		= get_conf_val_int( nice, 		["nice"] )
+	retrial 	= get_conf_val_int( retrial, 		["retrial","retry"] )
+	q 		= get_conf_val( q, 			["q"] )
+	unlimited_mem_wt= get_conf_val_bool( unlimited_mem_wt, 	["unlimited_mem_wt"] )
+
+	if ( cmd_line_arg_has_key("mem") ) memory = get_cmd_line_arg_val( "mem" )
+	if ( cmd_line_arg_has_key("system") ) use_system = get_cmd_line_arg_val( "system" )
+	if ( nice <= -20 ) nice = -20
+	if ( nice > 19 ) nice = 19	
+	if ( use_system != "" ) system = use_system.toLower()
+	if ( system == "slurm" ) system = "generic"
+	if ( q != "" ) queue = q
+
+	// cpus, mem and timeout are pre-declared BDS variables for default resource settings
+	mem 		= get_res_mem(memory,1)
+	timeout 	= get_res_wt(wt)
+	retry 		= retrial
+
+	// do not modify this (BDS timeout; how long BDS will wait for tasks to be queued on the cluster)
+	walltimeout 	= 3600*24*100 	// timeout var. in BigDataScript (100 days, jobs will never be stopped by BDS due to BDS timeout)
+
+	print("\n\n== cluster/system info\n")
+	print( "Walltime (general)\t\t: $wt\n" )
+	print( "Max. memory (general)\t\t: $memory\n" )
+	print( "Force to use a system\t\t: $use_system\n" )
+	print( "Process priority (niceness)\t: $nice\n" )
+	print( "Retiral for failed tasks\t: $retrial\n" )
+	print( "Submit tasks to a cluster queue\t: $q\n" )
+	print( "Unlimited cluster mem./walltime\t: $unlimited_mem_wt\n")	
+}
+
+int get_res_wt( string str ) {
+	return (unlimited_mem_wt || is_system_local() ) ? -1 : parse_time( str )
+}
+
+int get_res_mem( string str, int n ) {
+	if ( n < 1 ) n = 1
+	return (unlimited_mem_wt || is_system_local() ) ? -1 : parse_mem( str )/n
+}
+
+int get_res_mem( string str ) {
+	return get_res_mem( str , 1 )
+}
+
+bool is_system_sge() {
+	return system == "sge"
+}
+
+bool is_system_local() {
+	return system == "local"
+}
+
+bool is_system_generic() {
+	return system == "generic"
+}
+
+bool is_system_slurm() { 
+	// slurm uses generic cluster, it's configured in bds.config and ./utils/clusterGeneral
+	return system == "generic"
+}