diff --git a/examples/README.html b/examples/README.html
index 792bf71effc68b1b099f2a631f8cab8c7e22fce7..1d017cdb9c5527cbedd80a3e17bdee1f029a2567 100644
--- a/examples/README.html
+++ b/examples/README.html
@@ -336,12 +336,14 @@
<p class="tab">
<b>NOTE:</b> the path for locating the chr_NC_gi/chr_hdr files is hard-coded in the <i>pwm_bowtie_wrapper</i> and <i>pwm_mscan_wrapper</i> scripts (variable chrNC_dir).
<p class="tab">
- The path to the <i>matrix_scan_parallel.py</i> script is defined by the bin_dir varaible.
+ The path to the <i>matrix_scan_parallel.py</i> script is defined by the bin_dir varaible.
<p class="tab">
These shell wrapper scripts have been written to make it easier to run the whole pwmscan pipeline, irrespecitve of what method one chooses.
<br>
They require to specify the matrix file containing the integer PWM, the p-value, the path to the assembly or index files, and the UCSC assembly name for the species (e.g. hg19, hg38, mm9 etc.).
- Optional parameters are the <non-overlapping matches> and <parallelize matrix_scan> flags that are by default set to 1.
+ Optional parameters are the overlapping matches flag (-o), the forward scanning (-f) option, the parallel execution (-p) and write to file (-w) flags.
+ <br>
+ If the '-o' option is specified, overlapping matches are retained. The '-f' option activates forward scanning. If the '-p' option is set, the matrix_scan program execution is done by parallel dispatch. If the -w option is set, the results are written to file, else output goes to STDOUT.
<br>
The output file is a match list in <b>BEDdetail</b> format.
<br>
@@ -367,51 +369,75 @@
<ul style='list-style-type: none;'>
<li><font color=firebrick><b>1) Bowtie (pwm_bowtie_wrapper)</b> script</font>
<p class="code">
- # Usage: pwm_bowtie_wrapper <matrix-file> <p-value> <bowtie-dir> <idx-file> <assembly[hg19|mm9|..]> [<non-overlapping matches [1/0]>]
+ # Usage: <br>pwm_bowtie_wrapper -m <matrix-file> -e <p-value> -d <bowtie-dir> -s <assembly[hg19|mm9|..]> [options]<br>
+ - version 1.1.4<br>
+ where options are:<br>
+ -f Scan sequences in forward direction [def=bidirectional]<br>
+ -o Allow for overlapping matches [def=non-overlapping matches]<br>
+ -w Write output to file [def=STDOUT]<br>
+ -p Distribute matrix_scan on multiple CPU-cores [def=non-parallel processing]<br>
<br>
<font color=firebrick># Example:</font>
<br>
- <i>pwm_bowtie_wrapper chen10_ctcf.mat 0.00001 /home/local/db/bowtie h_sapiens_hg19 hg19 0</i>
+ <i>pwm_bowtie_wrapper -m chen10_ctcf.mat -e 0.00001 -d /home/local/db/bowtie -s hg19 -o -w</i>
+ <br>
+ <font color=firebrick># List of matches (BED format) : chen10_ctcf_co1128_bowtie.bed</font>
<br>
- <font color=firebrick># Output file (BED) : chen10_ctcf_co1128_bowtie.bed</font>
+ <font color=firebrick># Total number of PWM matches : 143597</font>
<br>
</p>
</li><li><font color=firebrick><b>2) Matrix Scan (pwm_mscan_wrapper)</b> script</font>
<p class="code">
- # Usage: pwm_mscan_wrapper <matrix-file> <p-value> <genome-dir> <assembly[hg19|mm9|..]> [<non-overlapping matches [1/0]>] [<parallel [1/0]>]
+ # Usage:<br> pwm_mscan_wrapper -m <matrix-file> -e <p-value> -d <genome-dir> -s <assembly[hg19|mm9|..]> [options]<br>
+ - version 1.1.4<br>
+ where options are:<br>
+ -f Scan sequences in forward direction [def=bidirectional]<br>
+ -o Allow for overlapping matches [def=non-overlapping matches]<br>
+ -w Write output to file [def=STDOUT]<br>
+ -p Distribute matrix_scan on multiple CPU-cores [def=non-parallel processing]<br>
<br>
<font color=firebrick># Example:</font>
<br>
- <i>pwm_mscan_wrapper chen10_ctcf.mat 0.00001 /home/local/db/genome hg19 0 1</i>
+ <i>pwm_mscan_wrapper -m chen10_ctcf.mat -e 0.00001 -d /home/local/db/genome -s hg19 -o -p -w</i>
<br>
- <font color=firebrick># Output file (BED) : chen10_ctcf_co1128_matrix_scan.bed</font>
-
+ <font color=firebrick># List of matches (BED format) : chen10_ctcf_co1128_matrix_scan.bed</font>
+ <br>
+ <font color=firebrick># Total number of PWM matches : 143597</font>
<br>
</p>
</li><li><font color=firebrick><b>3) PWM Scan (pwm_scan)</b> script</font>
<p class="code">
- # Usage: pwm_scan <matrix-file> <p-value> <assembly[hg19|mm9|..]> <genome-root-dir [ex:/db]> [<non-overlapping matches [1/0]>] [<parallel [1/0]>]
+ # Usage:<br> pwm_scan -m <matrix-file> -e <p-value> -d <genome-root-dir> -s <assembly[hg19|mm9|..]> [options]<br>
+ - version 1.1.4<br>
+ where options are:<br>
+ -f Scan sequences in forward direction [def=bidirectional]<br>
+ -o Allow for overlapping matches [def=non-overlapping matches]<br>
+ -w Write output to file [def=STDOUT]<br>
+ -p Distribute matrix_scan on multiple CPU-cores [def=non-parallel processing]<br>
<br>
<font color=firebrick># Examples:</font>
<br>
- <i>pwm_scan chen10_ctcf.mat 0.00001 hg19 /home/local/db 0 0</i>
+ <i>pwm_scan -m examples/chen10_ctcf.mat -e 0.00001 -d /home/local/db -s hg19 -o -w</i>
+ <br>
+ <font color=firebrick># List of matches (BED format) : chen10_ctcf_co1128_bowtie.bed</font>
<br>
<font color=firebrick># Or if we set the <parallel> option:</font>
<br>
- <i>pwm_scan chen10_ctcf.mat 0.00001 hg19 /home/local/db 0 1</i>
+ <i>pwm_scan -m examples/chen10_ctcf.mat -e 0.00001 -d /home/local/db -s hg19 -o -w -p</i>
<br>
- <font color=firebrick># Output file (BED) : chen10_ctcf_co1128_matrix_scan.bed</font>
+ <font color=firebrick># List of matches (BED format) : chen10_ctcf_co1128_matrix_scan.bed</font>
<br>
<font color=firebrick># Also try the following commands (with p-value=10<sup>-6</sup>):</font>
<br>
- <i>pwm_scan chen10_ctcf.mat 0.000001 hg19 /home/local/db</i><br>
- <i>pwm_scan stat1.mat 0.000001 hg19 /home/local/db</i>
+ <i>pwm_scan -m chen10_ctcf.mat -e 0.000001 -d /home/local/db -s hg19 -w -p</i>
+ <br>
+ <font color=firebrick># List of matches (BED format) : chen10_ctcf_co1549_bowtie.bed</font>
<br>
- <font color=firebrick># Output file (BED) : chen10_ctcf_co1549_bowtie.bed</font>
+ <i>pwm_scan -m stat1.mat -e 0.000001 -d /home/local/db -s hg19 -o -w -p</i>
<br>
- <font color=firebrick># Output file (BED) : stat1_co1731_bowtie.bed</font>
+ <font color=firebrick># List of matches (BED format) : stat1_co1731_bowtie.bed</font>
<br>
<font color=firebrick># N.B.</font>
<br>
diff --git a/examples/README.txt b/examples/README.txt
index 8f8e6afd30c81bc52059f0a6498690e2ef00ffc0..5d214a0f4fd5f2055e946cd37b9e884752a3a13e 100644
--- a/examples/README.txt
+++ b/examples/README.txt
@@ -279,7 +279,8 @@ Shell (bash) Wrappers to execute the analysis pipeline
These shell wrapper scripts have been written to make it easier to run the whole pwmscan pipeline, irrespecitve of what method one chooses.
They require to specify the matrix file containing the integer PWM, the p-value, the path to the assembly or index files, and the UCSC assembly name for the species (e.g. hg19, hg38, mm9 etc.).
- Optional parameters are the <non-overlapping matches> and <parallelize matrix_scan> flags that are by default set to 1.
+ Optional parameters are the overlapping matches flag (-o), the forward scanning (-f) option, the parallel execution (-p) and write to file (-w) flags.
+ If the '-o' option is specified, overlapping matches are retained. The '-f' option activates forward scanning. If the '-p' option is set, the matrix_scan program execution is done by parallel dispatch. If the -w option is set, the results are written to file, else output goes to STDOUT.
The output file is a match list in BEDdetail format.
@@ -303,47 +304,74 @@ Shell (bash) Wrappers to execute the analysis pipeline
1) Bowtie wrapper script
- # Usage: pwm_bowtie_wrapper <matrix-file> <p-value> <bowtie-dir> <genome-idx-file> <assembly[hg19|mm9|..]> [<non-overlapping matches [1/0]>]
+ # Usage:
+ pwm_bowtie_wrapper -m <matrix-file> -e <p-value> -d <bowtie-dir> -s <assembly[hg19|mm9|..]> [options]
+ - version 1.1.4
+ where options are:
+ -f Scan sequences in forward direction [def=bidirectional]
+ -o Allow for overlapping matches [def=non-overlapping matches]
+ -w Write output to file [def=STDOUT]
+ -p Distribute matrix_scan on multiple CPU-cores [def=non-parallel processing]
# Example:
- pwm_bowtie_wrapper chen10_ctcf.mat 0.00001 /home/local/db/bowtie h_sapiens_hg19 hg19 0
+ pwm_bowtie_wrapper -m chen10_ctcf.mat -e 0.00001 -d /home/local/db/bowtie -s hg19 -o -w
- # Output file (BED) : chen10_ctcf_co1128_bowtie.bed
+ # List of matches (BED format) : chen10_ctcf_co1128_bowtie.bed
+ # Total number of PWM matches : 143597
------------------------------------------------------------------------------
2) Matrix Scan (matrix_scan) wrapper script
- # Usage: pwm_mscan_wrapper <matrix-file> <p-value> <genome-dir> <assembly[hg19|mm9|..]> [<non-overlapping matches [1/0]>] [<parallel [1/0]>]
+ # Usage:
+ pwm_mscan_wrapper -m <matrix-file> -e <p-value> -d <genome-dir> -s <assembly[hg19|mm9|..]> [options]
+ - version 1.1.4
+ where options are:
+ -f Scan sequences in forward direction [def=bidirectional]
+ -o Allow for overlapping matches [def=non-overlapping matches]
+ -w Write output to file [def=STDOUT]
+ -p Distribute matrix_scan on multiple cores [def=non-parallel processing]
# Example:
- pwm_mscan_wrapper chen10_ctcf.mat 0.00001 /home/local/db/genome hg19 0 1
+ pwm_mscan_wrapper -m chen10_ctcf.mat -e 0.00001 -d /home/local/db/genome -s hg19 -o -p -w
- # Output file (BED) : chen10_ctcf_co1128_matrix_scan.bed
+ # List of matches (BED format) : chen10_ctcf_co1128_matrix_scan.bed
+ # Total number of PWM matches : 143597
------------------------------------------------------------------------------
3) PWM Scan (pwm_scan) wrapper script
- # Usage: pwm_scan <matrix-file> <p-value> <assembly[hg19|mm9|..]> <genome-root-dir [ex:/db]> [<non-overlapping matches [1/0]>] [<parallel [1/0>]>]
+ # Usage:
+ pwm_scan -m <matrix-file> -e <p-value> -d <genome-root-dir> -s <assembly[hg19|mm9|..]> [options]
+ - version 1.1.4
+ where options are:
+ -f Scan sequences in forward direction [def=bidirectional]
+ -o Allow for overlapping matches [def=non-overlapping matches]
+ -w Write output to file [def=STDOUT]
+ -p Distribute matrix_scan on multiple CPU-cores [def=non-parallel processing]
# Examples:
- pwm_scan chen10_ctcf.mat 0.00001 hg19 /home/local/db 0 0
+ pwm_scan -m examples/chen10_ctcf.mat -e 0.00001 -d /home/local/db -s hg19 -o -w
+
+ # List of matches (BED format) : chen10_ctcf_co1128_bowtie.bed
# Or if we set the <parallel> option:
- pwm_scan chen10_ctcf.mat 0.00001 hg19 /home/local/db 0 1
+ pwm_scan -m examples/chen10_ctcf.mat -e 0.00001 -d /home/local/db -s hg19 -o -w -p
- # Output file (BED) : chen10_ctcf_co1128_matrix_scan.bed
+ # List of matches (BED format) : chen10_ctcf_co1128_matrix_scan.bed
# Also try the following commands (with p-value=10^6):
- pwm_scan chen10_ctcf.mat 0.000001 hg19 /home/local/db
- pwm_scan stat1.mat 0.000001 hg19 /home/local/db
+ pwm_scan -m chen10_ctcf.mat -e 0.000001 -d /home/local/db -s hg19 -w -p
+
+ # List of matches (BED format) : chen10_ctcf_co1549_bowtie.bed
+
+ pwm_scan -m stat1.mat -e 0.000001 -d /home/local/db -s hg19 -o -w -p
- # Output file (BED) : chen10_ctcf_co1549_bowtie.bed
- # Output file (BED) : stat1_co1731_bowtie.bed
+ # List of matches (BED format) : stat1_co1731_bowtie.bed
# N.B.
# If the <parallel> option is set (default), matrix_scan is run in parallel by splitting the processing of the entire genome on multiple processes for each chromosome file in parallel.
diff --git a/pwm_bowtie_wrapper b/pwm_bowtie_wrapper
index f2019ce8e846e23211e8908e2e843f07253b899e..a7058835507239e932957a45211fa5948a101e0b 100755
--- a/pwm_bowtie_wrapper
+++ b/pwm_bowtie_wrapper
@@ -1,30 +1,174 @@
#!/bin/bash
-# pwm_bowtie_wrapper: scan a genome with a PWM (using bowtie)
+# pwm_bowtie_wrapper: scan a genome with a PWM (using Bowtie)
# arguments:
# matrix-file
# p-value
# bowtie-dir
-# genome-idx-file
# assembly
-# non-overlapping matches [1/0]
+#
+# forward scanning [-f]
+# non-overlapping matches [DEFAULT]
+# parallelize matrix_scan [-p]
+# write output to file [-w]
#
# 25.11.2015 Giovanna Ambrosini
#
+# The <assembly> argument corresponds to the UCSC assembly name (e.g. hg19, mm9).
+# The <bowtie-dir> is the root directory of the genome files (i.e. the bowtie index
+# files).
+#
+# The program writes the output to STDOUT, unless the -w option is specified, in which
+# case results are written to file.
+#
+# The variable chrNC_dir (used for locating the the chr_NC_gi/chr_hdr files) is set to
+# chrNC_dir=/home/local/db/genome
+#
# 05.12.2017 Giovanna Ambrosini
# Optimize Bowtie pipeline
# Use new version of the mba program
-#
+# 12.12.2017 Giovanna Ambrosini
+# Use getopts to parse script options
+# Add forward scanning option
+# Add write to file option
+# Remove <genome-idx-file> argument
+# (the variable is defined via the assoctaive array BOWTIEIDX)
-ARGS=5 # Script requires 5 arguments.
E_BADARGS=85 # Wrong number of arguments passed to script.
-if [ $# -lt "$ARGS" ]
+fwd_flag=""
+fwd_str=""
+forward=0
+w_flag=0
+non_overlapping=1
+parallel=0
+
+matrix_file=""
+p_value=""
+bowtie_dir=""
+genome_idx_file=""
+assembly=""
+
+display_usage() {
+ echo "Usage:"
+ echo " pwm_bowtie_wrapper -m <matrix-file> -e <p-value> -d <bowtie-dir> -s <assembly[hg19|mm9|..]> [options]"
+ echo " - version 1.1.4"
+ echo " where options are:"
+ echo " -f Scan sequences in forward direction [def=bidirectional]"
+ echo " -o Allow for overlapping matches [def=non-overlapping matches]"
+ echo " -w Write output to file [def=STDOUT]"
+ echo -e " -p Distribute matrix_scan on multiple CPU-cores [def=non-parallel processing]\n"
+ echo "Please report bugs to Giovanna.Ambrosini@epfl.ch"
+}
+# Parse options
+while getopts ":hfopwm:e:d:s:" opt; do
+ case ${opt} in
+ h )
+ display_usage
+ exit 0
+ ;;
+ f )
+ forward=1
+ ;;
+ o )
+ non_overlapping=0
+ ;;
+ p )
+ parallel=1
+ ;;
+ m )
+ matrix_file=$OPTARG
+ ;;
+ e )
+ p_value=$OPTARG
+ ;;
+ d )
+ bowtie_dir=$OPTARG
+ ;;
+ s )
+ assembly=$OPTARG
+ ;;
+ w )
+ w_flag=1
+ ;;
+ :)
+ echo "Option -$OPTARG requires an argument." >&2
+ exit 1
+ ;;
+ \? )
+ echo "Invalid Option: -$OPTARG" 1>&2
+ exit 1
+ ;;
+ esac
+done
+shift $((OPTIND -1))
+
+parsed_args=$((OPTIND -1))
+
+if [ $parsed_args -le 1 ]
+then
+ display_usage
+ exit 1
+fi
+
+#echo "Matrix file : $matrix_file" >&2
+#echo "P-value : $p_value" >&2
+#echo "bowtie-dir : $bowtie_dir" >&2
+#echo "assembly : $assembly" >&2
+
+if [ -z "$matrix_file" ]
+then
+ echo "Please specify the PWM file (-m)" >&2
+ exit 1
+fi
+if [ -z "$p_value" ]
+then
+ echo "Please specify the p-value (-e)" >&2
+ exit 1
+fi
+if [ -z "$bowtie_dir" ]
+then
+ echo "Please specify the bowtie genome directory (-d)" >&2
+ exit 1
+fi
+if [ -z "$assembly" ]
then
- echo "Usage: `basename $0` <matrix-file> <p-value> <bowtie-dir> <genome-idx-file> <assembly[hg19|mm9|..]> [<non-overlapping matches [1/0]>]"
- exit $E_BADARGS
+ echo "Please specify the species assembly (-s)" >&2
+ exit 1
fi
+# Define path for files chr_NC_gi/chr_hdr (needed for BED-format conversion)
+chrNC_dir=/home/local/db/genome
+
+# Declare associative array for bowtie index files
+declare -A BOWTIEIDX
+BOWTIEIDX[hg19]=h_sapiens_hg19
+BOWTIEIDX[hg18]=h_sapiens_hg18
+BOWTIEIDX[hg38]=h_sapiens_hg38
+BOWTIEIDX[mm9]=m_musculus_mm9
+BOWTIEIDX[mm10]=m_musculus_mm10
+BOWTIEIDX[bosTau3]=b_taurus_bosTau3
+BOWTIEIDX[bosTau8]=b_taurus_bosTau8
+BOWTIEIDX[canFam2]=c_familiaris_canFam2
+BOWTIEIDX[canFam3]=c_familiaris_canFam3
+BOWTIEIDX[panTro2]=p_troglodytes_panTro2
+BOWTIEIDX[panTro5]=p_troglodytes_panTro5
+BOWTIEIDX[rn5]=r_norvegicus_rn5
+BOWTIEIDX[rn6]=r_norvegicus_rn6
+BOWTIEIDX[dm3]=d_melanogaster_dm3
+BOWTIEIDX[dm6]=d_melanogaster_dm6
+BOWTIEIDX[ce6]=c_elegans_ce6
+BOWTIEIDX[ce10]=c_elegans_ce10
+BOWTIEIDX[ce11]=c_elegans_ce11
+BOWTIEIDX[danRer5]=d_rerio_danRer5
+BOWTIEIDX[danRer7]=d_rerio_danRer7
+BOWTIEIDX[danRer10]=d_rerio_danRer10
+BOWTIEIDX[sacCer2]=s_cerevisiae_sacCer2
+BOWTIEIDX[sacCer3]=s_cerevisiae_sacCer3
+
+#echo ${BOWTIEIDX[$assembly]}
+genome_idx_file=${BOWTIEIDX[$assembly]}
+
# Declare associative array for background base composition
declare -A BGFREQ
BGFREQ[hg19]=0.29,0.21,0.21,0.29
@@ -51,92 +195,95 @@ BGFREQ[danRer10]=0.32,0.18,0.18,0.32
BGFREQ[sacCer2]=0.31,0.19,0.19,0.31
BGFREQ[sacCer3]=0.31,0.19,0.19,0.31
-echo "BG nucleotide composition: ${BGFREQ[$5]}"
-bg_freq=${BGFREQ[$5]}
-
-chrNC_dir=/home/local/db/genome
+echo "BG nucleotide composition: ${BGFREQ[$assembly]}" >&2
+bg_freq=${BGFREQ[$assembly]}
-non_overlapping=1
if [ $# == 6 ]
then
non_overlapping=$6
fi
-if [ -f "$1" ]
+if [ ! -f "$matrix_file" ]
then
- matrix_file=$1
-else
- echo "File \"$1\" does not exist."
+ echo "Matrix File \"$1\" does not exist." >&2
exit $E_BADARGS
fi
-p_value=$2
-bowtie_dir=$3
-
-if [ -f "$3/$4.1.ebwt" ]
+if [ ! -f "$bowtie_dir/$genome_idx_file.1.ebwt" ]
then
- genome_idx_file=$4
-else
- echo "File \"$3/$4.1.ebwt\" does not exist."
+ echo "Genome File \"$1\" does not exist." >&2
exit $E_BADARGS
fi
-assembly=$5
-
# Extract basename from matrix file (without path)
matrix_name=$(basename "$matrix_file")
extension="${matrix_name##*.}"
matrix_name="${matrix_name%.*}"
-echo "PWM name: $matrix_name"
+echo "PWM name: $matrix_name" >&2
-echo "======== Computing PWM length ========"
+echo "======== Computing PWM length ========" >&2
matrix_len=$(cat $matrix_file | perl -ane 'next if (/^#/ or /^>/); print;' | wc -l)
file_len=$(cat $matrix_file | wc -l)
-echo "PWM file: $matrix_file"
-echo "PWM length: $matrix_len"
-echo "PWM file length : $file_len lines"
+echo "PWM file: $matrix_file" >&2
+echo "PWM length: $matrix_len" >&2
+echo "PWM file length : $file_len lines" >&2
-echo "======== Calculating PWM score ========"
+echo "======== Calculating PWM score ========" >&2
matrix_score=$(matrix_prob -e $p_value --bg "$bg_freq" $matrix_file \
| grep SCORE | sed 's/:/\ /'\
| awk -F " " '{print $2}')
-echo "PWM score: $matrix_score"
+echo "PWM score: $matrix_score" >&2
+
+if [ $forward == 1 ]
+then
+ echo "Scanning in forward direction..." >&2
+ fwd_str="fwd_"
+ fwd_flag="--norc"
+fi
if [ $non_overlapping == 1 ]
then
- echo "Output non-overlapping matches..."
+ echo "Output non-overlapping matches..." >&2
fi
-echo "======== Generating PWM score distribution ========"
+echo "======== Generating PWM score distribution ========" >&2
# Generate the Matrix Score Cumulative Table
pwmScore_tab=${matrix_name}_co${matrix_score}_scoretab.txt
matrix_prob --bg "$bg_freq" $matrix_file 2>/dev/null > $pwmScore_tab
-echo "PWM distribution score: $pwmScore_tab"
+echo "PWM distribution score: $pwmScore_tab" >&2
# Run the Bowtie pipeline
-pwmout_bed=${matrix_name}_co${matrix_score}_bowtie.bed
+if [ $w_flag == 1 ]
+then
+ pwmout_bed=${matrix_name}_co${matrix_score}_${fwd_str}bowtie.bed
+else
+ pwmout_bed="/dev/stdout"
+fi
-echo "======== Bowtie-based pipeline ========"
+echo "======== Bowtie-based pipeline ========" >&2
if [ $non_overlapping == 0 ]
then
- echo "mba -c $matrix_score $matrix_file | awk '{print \">\"\$2\"\n\"\$1}' | bowtie --threads 4 -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' > $pwmout_bed"
- echo "..."
- mba -c $matrix_score $matrix_file | awk '{print ">"$2"\n"$1}' | bowtie --threads 4 -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' > $pwmout_bed
+ echo "mba -c $matrix_score $matrix_file | awk '{print \">\"\$2\"\n\"\$1}' | bowtie --threads 4 $fwd_flag -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' >$pwmout_bed" >&2
+ echo "..." >&2
+ mba -c $matrix_score $matrix_file | awk '{print ">"$2"\n"$1}' | bowtie --threads 4 $fwd_flag -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' >$pwmout_bed
else
- echo "mba -c $matrix_score $matrix_file | awk '{print \">\"\$2\"\n\"\$1}' | bowtie --threads 4 -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | filterOverlaps -l$matrix_len | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' > $pwmout_bed"
- echo "..."
- mba -c $matrix_score $matrix_file | awk '{print ">"$2"\n"$1}' | bowtie --threads 4 -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | filterOverlaps -l$matrix_len | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' > $pwmout_bed
+ echo "mba -c $matrix_score $matrix_file | awk '{print \">\"\$2\"\n\"\$1}' | bowtie --threads 4 $fwd_flag -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | filterOverlaps -l$matrix_len | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' >$pwmout_bed" >&2
+ echo "..." >&2
+ mba -c $matrix_score $matrix_file | awk '{print ">"$2"\n"$1}' | bowtie --threads 4 $fwd_flag -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | filterOverlaps -l$matrix_len | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' >$pwmout_bed
fi
-echo "List of matches (BED format) : $pwmout_bed"
-nb_hits=`cat $pwmout_bed|wc -l`
-echo "Total number of PWM matches : $nb_hits"
+if [ $w_flag == 1 ]
+then
+ echo "List of matches (BED format) : $pwmout_bed" >&2
+ nb_hits=`cat $pwmout_bed|wc -l`
+ echo "Total number of PWM matches : $nb_hits" >&2
+fi
if [ -e unmapped.dat ]
then
diff --git a/pwm_mscan_wrapper b/pwm_mscan_wrapper
index 832b66926cc7ff1cc77247092bd298d1e8c32a9e..6f8589485ccbf3266215477aa235e515b32b7416 100755
--- a/pwm_mscan_wrapper
+++ b/pwm_mscan_wrapper
@@ -6,26 +6,146 @@
# p-value
# genome-dir
# assembly
-# non-overlapping matches [1/0]
-# parallel processing [1/0]
+#
+# forward scanning [-f]
+# non-overlapping matches [DEFAULT]
+# parallel processing [-p]
+# write output to file [-w]
#
# 25.11.2015 Giovanna Ambrosini
#
-# 05.12.2017 Giovanna Ambrosini
-# Use new faster version of the matrix_scan program
+# The <assembly> argument corresponds to the UCSC assembly name (e.g. hg19, mm9).
+# The <genome-dir> is the directory of the genome files (i.e. the FASTA chromosome
+# files used by matrix_scan).
+# If specified (via the -p option), matrix_scan can be run in parallel.
+# This is only possible on a multi-core processor architecture.
+#
+# The variable chrNC_dir (used for locating the the chr_NC_gi/chr_hdr files) is set to
+# chrNC_dir=$genome_dir
+#
+# The variable bin_dir is set to /home/local/bin
+#
+# The program writes the output to STDOUT, unless the -w option is specified, in which
+# case results are written to file.
#
+# 05.12.2017 Giovanna Ambrosini
+# Use the new faster version of the matrix_scan program
+# 07.12.2017 Giovanna Ambrosini
+# Use getopts to parse script options
+# Add forward scanning option
+# 12.12.2017 Giovanna Ambrosini
+# Use getopts to parse script options
+# Add forward scanning option
+# Add write to file option
-ARGS=4 # Script takes 6 arguments, but at least 4 arguments are mandatory.
E_BADARGS=85 # Wrong number of arguments passed to script.
-if [ $# -lt "$ARGS" ]
+fwd_flag=""
+fwd_str=""
+forward=0
+w_flag=0
+non_overlapping=1
+parallel=0
+
+matrix_file=""
+p_value=""
+genome_dir=""
+assembly=""
+
+display_usage() {
+ echo "Usage:"
+ echo " pwm_mscan_wrapper -m <matrix-file> -e <p-value> -d <genome-dir> -s <assembly[hg19|mm9|..]> [options]"
+ echo " - version 1.1.4"
+ echo " where options are:"
+ echo " -f Scan sequences in forward direction [def=bidirectional]"
+ echo " -o Allow for overlapping matches [def=non-overlapping matches]"
+ echo " -w Write output to file [def=STDOUT]"
+ echo -e " -p Distribute matrix_scan on multiple cores [def=non-parallel processing]\n"
+ echo "Please report bugs to Giovanna.Ambrosini@epfl.ch"
+}
+# Parse options
+while getopts ":hfopwm:e:d:s:" opt; do
+ case ${opt} in
+ h )
+ display_usage
+ exit 0
+ ;;
+ f )
+ forward=1
+ ;;
+ o )
+ non_overlapping=0
+ ;;
+ p )
+ parallel=1
+ ;;
+ m )
+ matrix_file=$OPTARG
+ ;;
+ e )
+ p_value=$OPTARG
+ ;;
+ d )
+ genome_dir=$OPTARG
+ ;;
+ s )
+ assembly=$OPTARG
+ ;;
+ w )
+ w_flag=1
+ ;;
+ :)
+ echo "Option -$OPTARG requires an argument." >&2
+ exit 1
+ ;;
+ \? )
+ echo "Invalid Option: -$OPTARG" 1>&2
+ exit 1
+ ;;
+ esac
+done
+shift $((OPTIND -1))
+
+parsed_args=$((OPTIND -1))
+
+if [ $parsed_args -le 1 ]
then
- echo "Usage: `basename $0` <matrix-file> <p-value> <genome-dir> <assembly[hg19|mm9|..]> [<non-overlapping matches [1/0]>] [<parallel [1/0]>]"
- exit $E_BADARGS
+ display_usage
+ exit 1
+fi
+
+#echo "Matrix file : $matrix_file" >&2
+#echo "P-value : $p_value" >&2
+#echo "genome-dir : $genome_dir" >&2
+#echo "assembly : $assembly" >&2
+#echo "Write flag : $w_flag" >&2
+
+if [ -z "$matrix_file" ]
+then
+ echo "Please specify the PWM file (-m)" >&2
+ exit 1
+fi
+if [ -z "$p_value" ]
+then
+ echo "Please specify the p-value (-e)" >&2
+ exit 1
+fi
+if [ -z "$genome_dir" ]
+then
+ echo "Please specify the genome directory (-d)" >&2
+ exit 1
+fi
+if [ -z "$assembly" ]
+then
+ echo "Please specify the species assembly (-s)" >&2
+ exit 1
fi
bin_dir=$(echo /home/local/bin)
+# Define path for files chr_NC_gi/chr_hdr (needed for BED-format conversion)
+chrNC_dir=$genome_dir
+
# Declare associative array for background base composition
declare -A BGFREQ
BGFREQ[hg19]=0.29,0.21,0.21,0.29
@@ -52,124 +172,114 @@ BGFREQ[danRer10]=0.32,0.18,0.18,0.32
BGFREQ[sacCer2]=0.31,0.19,0.19,0.31
BGFREQ[sacCer3]=0.31,0.19,0.19,0.31
-echo "BG nucleotide composition: ${BGFREQ[$4]}"
-bg_freq=${BGFREQ[$4]}
+echo "BG nucleotide composition: ${BGFREQ[$assembly]}" >&2
+bg_freq=${BGFREQ[$assembly]}
widx_size=""
-non_overlapping=1
-
-if [ $# -ge 5 ]
+if [ ! -f "$matrix_file" ]
then
- non_overlapping=$5
-fi
-
-parallel=1
-
-if [ $# == 6 ]
-then
- parallel=$6
-fi
-
-if [ -f "$1" ]
-then
- matrix_file=$1
-else
- echo "File \"$1\" does not exist."
+ echo "Matrix File \"$1\" does not exist."
exit $E_BADARGS
fi
-p_value=$2
-genome_dir=$3
-
-# Define path for files chr_NC_gi/chr_hdr (needed for BED-format conversion)
-chrNC_dir=$genome_dir
-
-assembly=$4
-
# Extract basename from matrix file (without path)
matrix_name=$(basename "$matrix_file")
extension="${matrix_name##*.}"
matrix_name="${matrix_name%.*}"
-echo "PWM name: $matrix_name"
+echo "PWM name: $matrix_name" >&2
-echo "======== Computing PWM length ========"
+echo "======== Computing PWM length ========" >&2
matrix_len=$(cat $matrix_file | perl -ane 'next if (/^#/ or /^>/); print;' | wc -l)
file_len=$(cat $matrix_file | wc -l)
-echo "PWM file: $matrix_file"
-echo "PWM length: $matrix_len"
-echo "PWM file length : $file_len"
+echo "PWM file: $matrix_file" >&2
+echo "PWM length: $matrix_len" >&2
+echo "PWM file length : $file_len" >&2
if [ "$matrix_len" -gt "17" ]
then
- if [ $(echo "$p_value > 0.0001 " | bc) -eq 1 ]
+ if [ $(echo "$p_value >= 0.0001 " | bc) -eq 1 ]
then
widx_size="-i 10"
echo "Word index size : $widx_size"
fi
fi
-echo "======== Calculating PWM score ========"
+echo "======== Calculating PWM score ========" >&2
matrix_score=$(matrix_prob -e $p_value --bg "$bg_freq" $matrix_file \
| grep SCORE | sed 's/:/\ /'\
| awk -F " " '{print $2}')
-echo "PWM score: $matrix_score"
-#echo "parallelize: $parallel"
+echo "PWM score: $matrix_score" >&2
# Run matrix_scan pipeline
+if [ $forward == 1 ]
+then
+ echo "Scanning in forward direction..." >&2
+ fwd_flag="-f"
+ fwd_str="fwd_"
+fi
+
if [ $parallel == 1 ]
then
- echo "Parallelize matrix_scan..."
+ echo "Parallelize matrix_scan..." >&2
fi
if [ $non_overlapping == 1 ]
then
- echo "Output non-overlapping matches..."
+ echo "Output non-overlapping matches..." >&2
fi
-echo "======== Generating PWM score distribution ========"
+echo "======== Generating PWM score distribution ========" >&2
# Generate the Matrix Score Cumulative Table
pwmScore_tab=${matrix_name}_co${matrix_score}_scoretab.txt
matrix_prob --bg "$bg_freq" $matrix_file 2>/dev/null > $pwmScore_tab
-echo "PWM distribution score: $pwmScore_tab"
+echo "PWM distribution score: $pwmScore_tab" >&2
# Run the matrix_scan pipeline
-pwmout_bed=${matrix_name}_co${matrix_score}_matrix_scan.bed
+if [ $w_flag == 1 ]
+then
+ pwmout_bed=${matrix_name}_co${matrix_score}_${fwd_str}matrix_scan.bed
+else
+ pwmout_bed="/dev/stdout"
+fi
-echo "======== matrix_scan-based pipeline ========"
+echo "======== matrix_scan-based pipeline ========" >&2
if [ $parallel == 0 ]
then
if [ $non_overlapping == 0 ]
then
- echo "cat $genome_dir/$assembly/chrom[^M]*.seq | $bin_dir/matrix_scan -m $matrix_file -c $matrix_score $widx_size | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' > $pwmout_bed"
- echo "..."
- cat $genome_dir/$assembly/chrom[^M]*.seq | $bin_dir/matrix_scan -m $matrix_file -c $matrix_score $widx_size | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' > $pwmout_bed
+ echo "cat $genome_dir/$assembly/chrom[^M]*.seq | $bin_dir/matrix_scan -m $matrix_file -c $matrix_score $fwd_flag $widx_size | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' >$pwmout_bed" >&2
+ echo "..." >&2
+ cat $genome_dir/$assembly/chrom[^M]*.seq | $bin_dir/matrix_scan -m $matrix_file -c $matrix_score $fwd_flag $widx_size | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' >$pwmout_bed
else
- echo "cat $genome_dir/$assembly/chrom[^M]*.seq | $bin_dir/matrix_scan -m $matrix_file -c $matrix_score $widx_size | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' > $pwmout_bed"
- echo "..."
- cat $genome_dir/$assembly/chrom[^M]*.seq | $bin_dir/matrix_scan -m $matrix_file -c $matrix_score $widx_size | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' > $pwmout_bed
+ echo "cat $genome_dir/$assembly/chrom[^M]*.seq | $bin_dir/matrix_scan -m $matrix_file -c $matrix_score $fwd_flag $widx_size | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' >$pwmout_bed" >&2
+ echo "..." >&2
+ cat $genome_dir/$assembly/chrom[^M]*.seq | $bin_dir/matrix_scan -m $matrix_file -c $matrix_score $fwd_flag $widx_size | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' >$pwmout_bed
fi
else
if [ $non_overlapping == 0 ]
then
- echo "python $bin_dir/matrix_scan_parallel.py -m $matrix_file -f \"\$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)\" -c $matrix_score $widx_size -p 15 | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' > $pwmout_bed"
- echo "..."
- python $bin_dir/matrix_scan_parallel.py -m $matrix_file -f "$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)" -c $matrix_score $widx_size -p 15 | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' > $pwmout_bed
+ echo "python $bin_dir/matrix_scan_parallel.py -m $matrix_file -s \"\$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)\" -c $matrix_score $fwd_flag $widx_size -p 15 | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' >$pwmout_bed" >&2
+ echo "..." >&2
+ python $bin_dir/matrix_scan_parallel.py -m $matrix_file -s "$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)" -c $matrix_score $fwd_flag $widx_size -p 15 | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' >$pwmout_bed
else
- echo "python $bin_dir/matrix_scan_parallel.py -m $matrix_file -f \"\$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)\" -c $matrix_score $widx_size -p 15 | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' > $pwmout_bed"
- echo "..."
- python $bin_dir/matrix_scan_parallel.py -m $matrix_file -f "$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)" -c $matrix_score $widx_size -p 15 | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' > $pwmout_bed
+ echo "python $bin_dir/matrix_scan_parallel.py -m $matrix_file -s \"\$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)\" -c $matrix_score $fwd_flag $widx_size -p 15 | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' >$pwmout_bed" >&2
+ echo "..." >&2
+ python $bin_dir/matrix_scan_parallel.py -m $matrix_file -s "$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)" -c $matrix_score $fwd_flag $widx_size -p 15 | sort -s -k1,1 -k2,2n | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' >$pwmout_bed
fi
fi
-echo "List of matches (BED format) : $pwmout_bed"
-nb_hits=`cat $pwmout_bed|wc -l`
-echo "Total number of PWM matches : $nb_hits"
+if [ $w_flag == 1 ]
+then
+ echo "List of matches (BED format) : $pwmout_bed" >&2
+ nb_hits=`cat $pwmout_bed|wc -l`
+ echo "Total number of PWM matches : $nb_hits" >&2
+fi
rm $pwmScore_tab
diff --git a/pwm_scan b/pwm_scan
index b1aacb5cd2856e49098595f0c1c7159728228751..d24edab4225b05fe79eae0aecd283a07832e7484 100755
--- a/pwm_scan
+++ b/pwm_scan
@@ -6,8 +6,11 @@
# p-value
# assembly
# genome-root-dir
-# non-overlapping matches [1/0]
-# parallelize matrix_scan [1/0]
+#
+# forward scanning [-f]
+# non-overlapping matches [DEFAULT]
+# parallelize matrix_scan [-p]
+# write output to file [-w]
#
# 19.09.2017 Giovanna Ambrosini
#
@@ -18,32 +21,142 @@
#
# the Bowtie-based method or the more conventional matrix_scan algotrithm
#
-# If specified (via the <parallel> option), matrix_scan can be run in parallel.
+# If specified (via the -p option), matrix_scan can be run in parallel.
# This is only possible if the machine has enough CPU power (multi-core processor).
# The <assembly> argument corresponds to the UCSC assembly name (e.g. hg19, mm9).
# The <genome-root-dir> is the root directory of the genome files (the bowtie index
# files, and the FASTA chromosome files used by matrix_scan).
# The <genome-root-dir> is supposed to have two sub-directories, bowtie and genome,
# for storing separately the bowtie indices and the chromosome files respectively.
+# The bowtie index files are declared in the BOWTIEIDX associative array.
# The variable chrNC_dir (used for locating the the chr_NC_gi/chr_hdr files) is set to
# chrNC_dir=$genome_dir
#
+#
+# The variable bin_dir is set to /home/local/bin
+#
+# The program writes the output to STDOUT, unless the -w option is specified, in which
+# case results are written to file.
+#
# 05.12.2017 Giovanna Ambrosini
# Optimize PWMscan pipeline
# Use new faster version of the matrix_scan program
-#
+# 12.12.2017 Giovanna Ambrosini
+# Use getopts to parse script options
+# Add forward scanning option
+# Add write to file option
-ARGS=4 # The script requires 4 mandatory arguments.
E_BADARGS=85 # Wrong number of arguments passed to the script.
-if [ $# -lt "$ARGS" ]
+fwd_flag=""
+fwd_str=""
+forward=0
+w_flag=0
+non_overlapping=1
+parallel=0
+
+matrix_file=""
+p_value=""
+genome_root_dir=""
+assembly=""
+
+display_usage() {
+ echo "Usage:"
+ echo " pwm_scan -m <matrix-file> -e <p-value> -d <genome-root-dir> -s <assembly[hg19|mm9|..]> [options]"
+ echo " - version 1.1.4"
+ echo " where options are:"
+ echo " -f Scan sequences in forward direction [def=bidirectional]"
+ echo " -o Allow for overlapping matches [def=non-overlapping matches]"
+ echo " -w Write output to file [def=STDOUT]"
+ echo -e " -p Distribute matrix_scan on multiple CPU-cores [def=non-parallel processing]\n"
+ echo "Please report bugs to Giovanna.Ambrosini@epfl.ch"
+}
+# Parse options
+while getopts ":hfopwm:e:d:s:" opt; do
+ case ${opt} in
+ h )
+ display_usage
+ exit 0
+ ;;
+ f )
+ forward=1
+ ;;
+ o )
+ non_overlapping=0
+ ;;
+ p )
+ parallel=1
+ ;;
+ m )
+ matrix_file=$OPTARG
+ ;;
+ e )
+ p_value=$OPTARG
+ ;;
+ d )
+ genome_root_dir=$OPTARG
+ ;;
+ s )
+ assembly=$OPTARG
+ ;;
+ w )
+ w_flag=1
+ ;;
+ :)
+ echo "Option -$OPTARG requires an argument." >&2
+ exit 1
+ ;;
+ \? )
+ echo "Invalid Option: -$OPTARG" 1>&2
+ exit 1
+ ;;
+ esac
+done
+shift $((OPTIND -1))
+
+parsed_args=$((OPTIND -1))
+
+if [ $parsed_args -le 1 ]
+then
+ display_usage
+ exit 1
+fi
+
+#echo "Matrix file : $matrix_file" >&2
+#echo "P-value : $p_value" >&2
+#echo "genome-root-dir: $genome_root_dir" >&2
+#echo "assembly : $assembly" >&2
+#echo "Write flag : $w_flag" >&2
+
+if [ -z "$matrix_file" ]
+then
+ echo "Please specify the PWM file (-m)" >&2
+ exit 1
+fi
+if [ -z "$p_value" ]
+then
+ echo "Please specify the p-value (-e)" >&2
+ exit 1
+fi
+if [ -z "$genome_root_dir" ]
+then
+ echo "Please specify the genome root directory (-d)" >&2
+ exit 1
+fi
+if [ -z "$assembly" ]
then
- echo "Usage: `basename $0` <matrix-file> <p-value> <assembly[hg19|mm9|..]> <genome-root-dir [ex:/db]> [<non-overlapping matches [1/0]>] [<parallel [1/0>]>]"
- exit $E_BADARGS
+ echo "Please specify the species assembly (-s)" >&2
+ exit 1
fi
bin_dir=$(echo /home/local/bin)
+bowtie_dir=$genome_root_dir"/bowtie"
+genome_dir=$genome_root_dir"/genome"
+
+# Define path for files chr_NC_gi/chr_hdr (needed for BED-format conversion)
+chrNC_dir=$genome_dir
+
# Declare associative array for bowtie index files
declare -A BOWTIEIDX
BOWTIEIDX[hg19]=h_sapiens_hg19
@@ -70,8 +183,8 @@ BOWTIEIDX[danRer10]=d_rerio_danRer10
BOWTIEIDX[sacCer2]=s_cerevisiae_sacCer2
BOWTIEIDX[sacCer3]=s_cerevisiae_sacCer3
-#echo ${BOWTIEIDX[$3]}
-genome_idx_file=${BOWTIEIDX[$3]}
+#echo ${BOWTIEIDX[$assembly]}
+genome_idx_file=${BOWTIEIDX[$assembly]}
# Declare associative array for background base composition
declare -A BGFREQ
@@ -99,81 +212,67 @@ BGFREQ[danRer10]=0.32,0.18,0.18,0.32
BGFREQ[sacCer2]=0.31,0.19,0.19,0.31
BGFREQ[sacCer3]=0.31,0.19,0.19,0.31
-echo "BG nucleotide composition: ${BGFREQ[$3]}"
-bg_freq=${BGFREQ[$3]}
+echo "BG nucleotide composition: ${BGFREQ[$assembly]}" >&2
+bg_freq=${BGFREQ[$assembly]}
widx_size=""
-non_overlapping=1
-
-if [ $# -ge 5 ]
-then
- non_overlapping=$5
-fi
-
-parallel=1
-if [ $# == 6 ]
+if [ ! -f "$matrix_file" ]
then
- parallel=$6
+ echo "Matrix File \"$1\" does not exist." >&2
+ exit $E_BADARGS
fi
-if [ -f "$1" ]
+if [ ! -f "$bowtie_dir/$genome_idx_file.1.ebwt" ]
then
- matrix_file=$1
-else
- echo "File \"$1\" does not exist."
+ echo "Genome File \"$1\" does not exist." >&2
exit $E_BADARGS
fi
-p_value=$2
-assembly=$3
-genome_root_dir=$4
-bowtie_dir=$genome_root_dir"/bowtie"
-genome_dir=$genome_root_dir"/genome"
-
-chrNC_dir=$genome_dir
-
# Extract basename from matrix file (without path)
matrix_name=$(basename "$matrix_file")
extension="${matrix_name##*.}"
matrix_name="${matrix_name%.*}"
-echo "PWM name: $matrix_name"
+echo "PWM name: $matrix_name" >&2
-echo "======== Computing PWM length ========"
+echo "======== Computing PWM length ========" >&2
matrix_len=$(cat $matrix_file | perl -ane 'next if (/^#/ or /^>/); print;' | wc -l)
file_len=$(cat $matrix_file | wc -l)
-echo "PWM file: $matrix_file"
-echo "PWM length: $matrix_len"
-echo "PWM file length : $file_len lines"
+echo "PWM file: $matrix_file" >&2
+echo "PWM length: $matrix_len" >&2
+echo "PWM file length : $file_len lines" >&2
if [ "$matrix_len" -gt "17" ]
then
if [ $(echo "$p_value > 0.0001 " | bc) -eq 1 ]
then
widx_size="-i 10"
- echo "Word index size : $widx_size"
+ echo "Word index size : $widx_size" >&2
fi
fi
-echo "======== Calculating PWM score ========"
+echo "======== Calculating PWM score ========" >&2
matrix_score=$(matrix_prob -e $p_value --bg "$bg_freq" $matrix_file \
| grep SCORE | sed 's/:/\ /'\
| awk -F " " '{print $2}')
-echo "PWM score: $matrix_score"
+echo "PWM score: $matrix_score" >&2
# Decide on search engine strategy
#
-# Check whether the p-value is too high (i.e. the raw score too low) given the motif lenght.
+# Check whether the p-value is high (i.e. the raw score too low) given the motif lenght.
# For long motifs and high p-values, the bowtie-based strategy becomes inefficient because
# the list of tags (representing the matrix with the given cut-off) becomes too large
# (of the order of half a billion or more)
#
+# If matrix_scan execution is done by parallel dispatch, the p-value threshold at which
+# matrix_scan becomes more effcient than Bowtie gets much lower
+#
use_matrix_scan=0
-echo "Parallelize: $parallel"
+echo "Parallelize: $parallel" >&2
if [ $parallel -eq 1 ]
then
if [ $matrix_len -gt 17 ]
@@ -203,75 +302,99 @@ else
fi
fi
-echo "Use matrix_scan: $use_matrix_scan"
+echo "Use matrix_scan: $use_matrix_scan" >&2
+
+if [ $forward == 1 ]
+then
+ echo "Scanning in forward direction..." >&2
+ fwd_str="fwd_"
+ if [ $use_matrix_scan == 1 ]
+ then
+ fwd_flag="-f"
+ else
+ fwd_flag="--norc"
+ fi
+fi
if [ $non_overlapping == 1 ]
then
- echo "Output non-overlapping matches..."
+ echo "Output non-overlapping matches..." >&2
fi
-echo "======== Generating PWM score distribution ========"
+echo "======== Generating PWM score distribution ========" >&2
# Generate the Matrix Score Cumulative Table
pwmScore_tab=${matrix_name}_co${matrix_score}_scoretab.txt
matrix_prob --bg "$bg_freq" $matrix_file 2>/dev/null > $pwmScore_tab
-echo "PWM distribution score: $pwmScore_tab"
+echo "PWM distribution score: $pwmScore_tab" >&2
# Run the PWMScan pipeline
if [ $use_matrix_scan -eq 0 ]
-then
+then
# Run Bowtie pipeline
- pwmout_bed=${matrix_name}_co${matrix_score}_bowtie.bed
- echo "======== Bowtie-based pipeline ========"
+ if [ $w_flag == 1 ]
+ then
+ pwmout_bed=${matrix_name}_co${matrix_score}_${fwd_str}bowtie.bed
+ else
+ pwmout_bed="/dev/stdout"
+ fi
+ echo "======== Bowtie-based pipeline ========" >&2
if [ $non_overlapping == 0 ]
then
- echo "mba -c $matrix_score $matrix_file | awk '{print \">\"\$2\"\n\"\$1}' | bowtie --threads 4 -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len $chrNC_dir | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' > $pwmout_bed"
- echo "..."
- mba -c $matrix_score $matrix_file | awk '{print ">"$2"\n"$1}' | bowtie --threads 4 -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' > $pwmout_bed
+ echo "mba -c $matrix_score $matrix_file | awk '{print \">\"\$2\"\n\"\$1}' | bowtie --threads 4 $fwd_flag -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' >$pwmout_bed" >&2
+ echo "..." >&2
+ mba -c $matrix_score $matrix_file | awk '{print ">"$2"\n"$1}' | bowtie --threads 4 $fwd_flag -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' >$pwmout_bed
else
- echo "mba -c $matrix_score $matrix_file | awk '{print \">\"\$2\"\n\"\$1}' | bowtie --threads 4 -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len $chrNC_dir | filterOverlaps -l$matrix_len | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' > $pwmout_bed"
- echo "..."
- mba -c $matrix_score $matrix_file | awk '{print ">"$2"\n"$1}' | bowtie --threads 4 -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | filterOverlaps -l$matrix_len | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' > $pwmout_bed
+ echo "mba -c $matrix_score $matrix_file | awk '{print \">\"\$2\"\n\"\$1}' | bowtie --threads 4 $fwd_flag -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | filterOverlaps -l$matrix_len | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' >$pwmout_bed" >&2
+ echo "..." >&2
+ mba -c $matrix_score $matrix_file | awk '{print ">"$2"\n"$1}' | bowtie --threads 4 $fwd_flag -l $matrix_len -n0 -a $bowtie_dir/$genome_idx_file -f - --un unmapped.dat | sort -s -k3,3 -k4,4n | bowtie2bed -s $assembly -l $matrix_len -i $chrNC_dir | filterOverlaps -l$matrix_len | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' >$pwmout_bed
fi
- echo "List of matches (BED format) : $pwmout_bed"
if [ -e unmapped.dat ]
then
rm unmapped.dat
fi
else
- pwmout_bed=${matrix_name}_co${matrix_score}_matrix_scan.bed
- echo "======== matrix_scan-based pipeline ========"
+ if [ $w_flag == 1 ]
+ then
+ pwmout_bed=${matrix_name}_co${matrix_score}_${fwd_str}matrix_scan.bed
+ else
+ pwmout_bed="/dev/stdout"
+ fi
+ echo "======== matrix_scan-based pipeline ========" >&2
# Run matrix_scan
if [ $parallel == 0 ]
then
if [ $non_overlapping == 0 ]
then
- echo "cat $genome_dir/$assembly/chrom*.seq | matrix_scan -m $matrix_file -c $matrix_score $widx_size | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' > $pwmout_bed"
- echo "..."
- cat $genome_dir/$assembly/chrom*.seq | matrix_scan -m $matrix_file -c $matrix_score $widx_size | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' > $pwmout_bed
+ echo "cat $genome_dir/$assembly/chrom*.seq | matrix_scan -m $matrix_file -c $matrix_score $fwd_flag $widx_size | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' >$pwmout_bed" >&2
+ echo "..." >&2
+ cat $genome_dir/$assembly/chrom*.seq | matrix_scan -m $matrix_file -c $matrix_score $fwd_flag $widx_size | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' >$pwmout_bed
else
- echo "cat $genome_dir/$assembly/chrom*.seq | matrix_scan -m $matrix_file -c $matrix_score $widx_size | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' > $pwmout_bed"
- echo "..."
- cat $genome_dir/$assembly/chrom*.seq | matrix_scan -m $matrix_file -c $matrix_score $widx_size | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' > $pwmout_bed
+ echo "cat $genome_dir/$assembly/chrom*.seq | matrix_scan -m $matrix_file -c $matrix_score $fwd_flag $widx_size | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' >$pwmout_bed" >&2
+ echo "..." >&2
+ cat $genome_dir/$assembly/chrom*.seq | matrix_scan -m $matrix_file -c $matrix_score $fwd_flag $widx_size | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' >$pwmout_bed
fi
else
if [ $non_overlapping == 0 ]
then
- echo "python $bin_dir/matrix_scan_parallel.py -m $matrix_file -f \"\$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)\" -c $matrix_score $widx_size -p 15 | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' > $pwmout_bed"
- echo "..."
- python $bin_dir/matrix_scan_parallel.py -m $matrix_file -f "$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)" -c $matrix_score $widx_size -p 15 | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' > $pwmout_bed
+ echo "python $bin_dir/matrix_scan_parallel.py -m $matrix_file -s \"\$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)\" -c $matrix_score $fwd_flag $widx_size -p 15 | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' >$pwmout_bed" >&2
+ echo "..." >&2
+ python $bin_dir/matrix_scan_parallel.py -m $matrix_file -s "$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)" -c $matrix_score $fwd_flag $widx_size -p 15 | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' >$pwmout_bed
else
- echo "python $bin_dir/matrix_scan_parallel.py -m $matrix_file -f \"\$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)\" -c $matrix_score $widx_size -p 15 | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' > $pwmout_bed"
- echo "..."
- python $bin_dir/matrix_scan_parallel.py -m $matrix_file -f "$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)" -c $matrix_score $widx_size -p 15 | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' > $pwmout_bed
+ echo "python $bin_dir/matrix_scan_parallel.py -m $matrix_file -s \"\$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)\" -c $matrix_score $fwd_flag $widx_size -p 15 | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk 'BEGIN { while((getline line < \"$pwmScore_tab\") > 0 ) {split(line,f,\" \"); pvalue[f[1]]=f[2]} close(\"$pwmScore_tab\")} {print \$1\"\t\"\$2\"\t\"\$3\"\t\"\$4\"\t\"\$5\"\t\"\$6\"\t\"\"$matrix_name\"\"\t\"\"P-value=\"pvalue[\$5]}' >$pwmout_bed" >&2
+ echo "..." >&2
+ python $bin_dir/matrix_scan_parallel.py -m $matrix_file -s "$(ls $genome_dir/$assembly/chrom*.seq|grep -v chromMt)" -c $matrix_score $fwd_flag $widx_size -p 15 | sort -s -k1,1 -k2,2n -k6,6 | mscan2bed -s $assembly -i $chrNC_dir | filterOverlaps -l$matrix_len | awk -v scoretab="$pwmScore_tab" -v pwmname="$matrix_name" 'BEGIN { while((getline line < scoretab) > 0 ) {split(line,f," "); pvalue[f[1]]=f[2]} close(scoretab)} {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"pwmname"\t""P-value="pvalue[$5]}' >$pwmout_bed
fi
fi
fi
-echo "List of matches (BED format) : $pwmout_bed"
-nb_hits=`cat $pwmout_bed|wc -l`
-echo "Total number of PWM matches : $nb_hits"
+if [ $w_flag == 1 ]
+then
+ echo "List of matches (BED format) : $pwmout_bed" >&2
+ nb_hits=`cat $pwmout_bed|wc -l`
+ echo "Total number of PWM matches : $nb_hits" >&2
+fi
# Remove temporary PWM file
rm $pwmScore_tab
diff --git a/python_tools/matrix_scan_parallel.py b/python_tools/matrix_scan_parallel.py
index b226397ec6b5420a7751ab6378a8c562b0b46917..5815ced30eafa6295a2cb6ac8892d7db3ce9f9d9 100755
--- a/python_tools/matrix_scan_parallel.py
+++ b/python_tools/matrix_scan_parallel.py
@@ -7,14 +7,14 @@ import tempfile
def check_options(options_obj):
"""
- Checks the presence of some values in the option fetch from the command
- line. Terminates the program if an expected value cannot be found.
+ Checks the arguments passed via the command line.
+ Terminates the program if an argument value is not found or is wrong.
:param options_obj: the return value of parser.parse_args()
:return: nothing
"""
# no option given (all default) -> print to use -h
if (options_obj.file_matrix is None) and \
- (options_obj.files_seq is None) and \
+ (options_obj.seq_files is None) and \
(options_obj.cutoff is None) and \
(options_obj.n_proc is 1):
print "use -h option for help"
@@ -25,8 +25,8 @@ def check_options(options_obj):
if options_dict["file_matrix"] is None:
sys.stderr.write("Option error! No matrix file was given (-m)!\n")
sys.exit(1)
- elif options_dict["files_seq"] is None:
- sys.stderr.write("Option error! No sequence files was given (-f)!\n")
+ elif options_dict["seq_files"] is None:
+ sys.stderr.write("Option error! No sequence files was given (-s)!\n")
sys.exit(1)
elif options_dict["cutoff"] is None:
sys.stderr.write("Option error! No cutoff score was given (-c)!\n")
@@ -41,7 +41,7 @@ def check_options(options_obj):
if __name__ == "__main__":
- # the dependencies
+ # the dependencies (paths to programs)
python_path = "/usr/bin/python" # the python 2.7 interpreter
parallelize_path = "/local/python/parallelize/parallelize.py" # the parall. dispatcher
matrix_scan_path = "/local/bin/matrix_scan" # matrix_scan
@@ -55,7 +55,7 @@ if __name__ == "__main__":
sys.stderr.write("Error! Cannot find %s!\n" % matrix_scan_path)
sys.exit(1)
- # parses options
+ # parse options
parser = optparse.OptionParser(version="v1.0",
description="This program executes matrix_scan in a parallel way. It will "
"spawn one process per sequence file to scan.",
@@ -72,14 +72,18 @@ if __name__ == "__main__":
help="The score cutoff.",
default=None,
metavar="score")
+ parser.add_option("-f", "--forward",
+ dest="fwd",
+ action="store_true",
+ help="Scan sequences in forward direction.")
parser.add_option("-i", "--windex",
dest="windex",
type="int",
help="The word index length, by default 7.",
default=None,
metavar="widx")
- parser.add_option("-f", "--files",
- dest="files_seq",
+ parser.add_option("-s", "--seqs",
+ dest="seq_files",
type="string",
default=None,
help="A list of sequence files to scan. File should be '\\n' separated (e.g. file1\\nfile2"
@@ -98,32 +102,34 @@ if __name__ == "__main__":
file_matrix = options.file_matrix
cutoff = options.cutoff
windex = options.windex
- files_seq = options.files_seq
- files_seq_list = [f for f in files_seq.split('\n')]
-
+ fwd = options.fwd
+ seq_files = options.seq_files
+ seq_files_list = [f for f in seq_files.split('\n')]
# check the sequence files
- for file_seq in files_seq_list:
- if not os.path.isfile(file_seq):
- sys.stderr.write("Error! %s does not exist!\n" % file_seq)
+ for seq_file in seq_files_list:
+ if not os.path.isfile(seq_file):
+ sys.stderr.write("Error! Sequence file %s does not exist!\n" % seq_file)
sys.exit(1)
# check the matrix file
if not os.path.isfile(file_matrix):
- sys.stderr.write("Error! %s does not exist!\n" % file_matrix)
+ sys.stderr.write("Error! Matrix file %s does not exist!\n" % file_matrix)
sys.exit(1)
# write a temporary file containing the arguments for each call of the program
# will be deleted as soon as closed
file_arg = tempfile.NamedTemporaryFile(mode="w+t")
- for file_seq in files_seq_list:
- if windex is not None:
- line = "-i %s\t-m %s\t-c %s\t%s\n" % (str(windex), file_matrix, str(cutoff), file_seq)
- file_arg.write(line)
- else:
- line = "-m %s\t-c %s\t%s\n" % (file_matrix, str(cutoff), file_seq)
- file_arg.write(line)
+ windex_str = ""
+ if windex is not None:
+ windex_str = "-i%s \t" % str(windex)
+ fwd_str = ""
+ if fwd:
+ fwd_str = "-f\t"
+ for seq_file in seq_files_list:
+ line = "%s%s-m %s\t-c %s\t%s\n" % (fwd_str, windex_str, file_matrix, str(cutoff), seq_file)
+ file_arg.write(line)
file_arg.flush()
- # execute the program using a parallel dispatch
+ # execute the program by parallel dispatch
try:
# for a reason I don't understand, using --simulate option raises "IOError: [Errno 32] Broken pipe"
command = "%s %s -s %s -i %s -p %s" % \