# HG changeset patch
# User messersc
# Date 1424342385 18000
# Node ID 243f75d0ed6edc3b62d9788c99d363cf0056bd0e
# Parent d42f4d78c85eed6aa4d3b0ca3871cc43afe7a7fb
Uploaded.
Includes new release 1.0.7 with fixed optional controls.
diff -r d42f4d78c85e -r 243f75d0ed6e JAMM.sh
--- a/JAMM.sh Wed Dec 17 10:40:23 2014 -0500
+++ b/JAMM.sh Thu Feb 19 05:39:45 2015 -0500
@@ -1,10 +1,23 @@
-#!/bin/bash
-########################################
-#### Peak finding Pipeline for NGS Data
-#### Bash script
-########################################
+########################################################################
+# JAMMv1.0.7rev1 is a peak finder for joint analysis of NGS replicates.
+# Copyright (C) 2014-2015 Mahmoud Ibrahim
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# Contact: mahmoud.ibrahim@mdc-berlin.de
+########################################################################
-export LC_ALL=en_US.UTF-8
##Finding out the path
sPath="`dirname \"$0\"`"
@@ -15,7 +28,7 @@
usage()
{
cat << EOF
-Welcome to JAMM v1.0.6rev2 (GNU GPLv3). Copyright (C) 2014 Mahmoud Ibrahim.
+Welcome to JAMM v1.0.7rev1 (GNU GPLv3). Copyright (C) 2014-2015 Mahmoud Ibrahim.
This program comes with ABSOLUTELY NO WARRANTY; for details visit http://www.gnu.org/licenses/gpl.html. This is free software, and you are welcome to redistribute it under certain conditions; visit http://www.gnu.org/licenses/gpl.html for details.
@@ -24,14 +37,17 @@
-g Genome size file (required)
-o Output directory (required)
-c directory containing input or Control files
+ -f Fragment length(s) (default: estimated)
+ -r Resolution, peak or region or window (default: peak)
-m Mode, normal or narrow (default: normal)
- -r Resolution, peak or region (default: peak)
- -f Fragment length(s) (default: estimated)
+ -i clustering Initialization window selection, deterministic or stochastic (default: deterministic)
+ -b Bin Size (default: estimated)
-w minimum Window size (default: 2 --- Note: this means minimum_window_size = bin_size x the_value_of_-w)
+ -e window Enrichment cutoff, auto or any numeric value (default: 1 --- Set this to "auto" to estimate the window enrichment cutoff)
+ -d keep PCR Dupicates in single-end mode, y or n (default: n --- if -t is "paired", this option has no effect)
+ -t Type, single or paired (default: single, requires BED files. paired requires BEDPE files)
-p Number of processors used by R scripts (default: 1)
- -b Bin Size (default: estimated)
- -t Type, single or paired (default: single, requires BED files. paired requires BEDPE files)
- -e window Enrichment cutoff, auto or any numeric value (default: 1 --- Set this to auto to estimate the appropriate window enrichment cutoff)
+
EOF
}
@@ -48,7 +64,8 @@
window="2"
type="single"
windowe="1"
-
+initModel="deterministic"
+uniq="n"
#Defaults -- Do not change
sdir=""
@@ -58,8 +75,10 @@
fraglen="ns"
ran=$RANDOM
wdir=$(mktemp -d)
+export LANG=C #locale defaults
+export LC_ALL=C #locale defaults
-while getopts "s:g:o:c:m:r:f:p:w:b:t:e:" OPTION
+while getopts "s:g:o:c:m:r:f:p:w:b:t:e:i:d:" OPTION
do
case $OPTION in
s) sdir=$OPTARG
@@ -86,6 +105,10 @@
;;
e) windowe=$OPTARG
;;
+ i) initModel=$OPTARG
+ ;;
+ d) uniq=$OPTARG
+ ;;
?)
usage
exit
@@ -104,14 +127,25 @@
usage
exit 1
fi
+if [[ -d "$out/peaks" ]]; then
+ printf "\n\nOutput directory $out/peaks already exists. I can't override existing results!\n\n"
+ exit 0
+fi
+if [ $fraglen == "ns" ]; then
+ if [[ -d "$out/xcorr" ]]; then
+ printf "\n\nOutput directory $out/xcorr already exists. I can't override existing results!\n\n"
+ exit 0
+ fi
+fi
#=======================> DONE!
+
# =============================
# Step One: Initial Processing
# =============================
-printf "\n\n========================================\nStarted JAMM Pipeline v1.0.6rev2...Hang on!\n========================================\n\n"
+printf "\n\n============================================\nStarted JAMM Pipeline v1.0.7rev1...Hang on!\n============================================\n\n"
if [ ! -d "$wdir" ]; then
mkdir $wdir #make working directory
@@ -161,17 +195,21 @@
fi
#determine average read lengths
-printf "Getting average read lengths..."
- if [ ! -z $bdir ]; then
- readC=$(awk '{a=$3-$2;print a;}' "$wdir/bkgd.$ran/ctrl.bed" | perl -lane '$a+=$_;END{print $a/$.}' | awk '{a=$1+0.5;print a;}' | cut -d"." -f1)
- fi
- readL=""
- for s in $sdir/*.bed; do #and for each sample file
- read=$(awk '{a=$3-$2;print a;}' "$s" | perl -lane '$a+=$_;END{print $a/$.}' | awk '{a=$1+0.5;print a;}' | cut -d"." -f1)
- readL="$readL,$read"
- done
- readL=${readL#","}
-printf "Done!\n"
+printf "Getting average read lengths...\n"
+readL=""
+if [ ! -z $bdir ]; then
+ readC=$(awk '{a=$3-$2;print a;}' "$wdir/bkgd.$ran/ctrl.bed" | perl -lane '$a+=$_;END{print $a/$.}' | awk '{a=$1+0.5;print a;}' | cut -d"." -f1)
+ printf "Control: $readC\n"
+fi
+readL=""
+for s in $sdir/*.bed; do #and for each sample file
+ file=$(basename $s)
+ samplefile=$(echo $file | awk -F"." '{print $1}');
+ read=$(awk '{a=$3-$2;print a;}' "$s" | perl -lane '$a+=$_;END{print $a/$.}' | awk '{a=$1+0.5;print a;}' | cut -d"." -f1)
+ printf "$samplefile: $read\n"
+ readL="$readL,$read"
+done
+readL=${readL#","}
#=======================> DONE!
@@ -299,7 +337,6 @@
#remove leading comma
samplelist=${samplelist#","}
frag=${frag#","}
-
Rscript "$sPath/bincalculator.r" -ibed="$samplelist" -s="$gsize" -rl="$readL" -d="$wdir" -p="$cores" -f="$frag" -type="$type"
fi
if [ $binsize != "ns" ]; then
@@ -360,7 +397,12 @@
read=$(echo $readL | cut -f "$k" -d ",")
k=$(($k+1))
fi
- perl "$sPath/readshifter.pl" "$wdir/samples.$ran/$samplefile" $shift $read > "$wdir/samples.$ran/ext.$samplefile"
+ if [ $uniq == "y" ]; then
+ perl "$sPath/readshifter.pl" "$wdir/samples.$ran/$samplefile" $shift $read > "$wdir/samples.$ran/ext.$samplefile"
+ fi
+ if [ $uniq == "n" ]; then
+ perl "$sPath/readshifter.pl" "$wdir/samples.$ran/$samplefile" $shift $read | uniq > "$wdir/samples.$ran/ext.$samplefile"
+ fi
fi
done
@@ -375,7 +417,13 @@
l=$(($dupnum+1))
bshift=$(echo $frag | cut -f "$l" -d ",")
fi
- perl "$sPath/readshifter.pl" "$wdir/bkgd.$ran/bkgd.$chr.ctrl.bed" $bshift $readC > "$wdir/bkgd.$ran/ext.bkgd.$chr.ctrl.bed"
+
+ if [ $uniq == "y" ]; then
+ perl "$sPath/readshifter.pl" "$wdir/bkgd.$ran/bkgd.$chr.ctrl.bed" $bshift $readC > "$wdir/bkgd.$ran/ext.bkgd.$chr.ctrl.bed"
+ fi
+ if [ $uniq == "n" ]; then
+ perl "$sPath/readshifter.pl" "$wdir/bkgd.$ran/bkgd.$chr.ctrl.bed" $bshift $readC | uniq > "$wdir/bkgd.$ran/ext.bkgd.$chr.ctrl.bed"
+ fi
bkgdfile="$wdir/bkgd.$ran/ext.bkgd.$chr.ctrl.bed"
fi
@@ -384,7 +432,7 @@
frag=${frag#","}
#call the peak calling R script
- Rscript "$sPath/peakfinder.r" -sfile="$f" -bednames="$samplelist" -frag="$frag" -bkgd=$bkgdfile -out="$wdir/peaks.$ran/" -clustnummer="$clustno" -resolution="$resol" -window="$window" -p="$cores" -bin="$binsize" -type="$type" -chrcount="$counting" -windowe="$windowe"
+ Rscript "$sPath/peakfinder.r" -sfile="$f" -bednames="$samplelist" -frag="$frag" -bkgd=$bkgdfile -out="$wdir/peaks.$ran/" -clustnummer="$clustno" -resolution="$resol" -window="$window" -p="$cores" -bin="$binsize" -type="$type" -chrcount="$counting" -initModel="$initModel" -windowe="$windowe"
counting=$(($counting+1));
cp "$wdir/peaks.$ran/$chr.peaks.bed" "$out/peaks/$chr.peaks.bed"
rm "$wdir/peaks.$ran/$chr.peaks.bed"
@@ -434,7 +482,7 @@
frag=${frag#","}
#call the peak calling R script
- Rscript "$sPath/peakfinder.r" -sfile=$f -bednames=$samplelist -frag="NA" -bkgd=$bkgdfile -out="$wdir/peaks.$ran/" -clustnummer="$clustno" -resolution="$resol" -window="$window" -p="$cores" -bin="$binsize" -type="$type" -chrcount="$counting" -windowe="$windowe"
+ Rscript "$sPath/peakfinder.r" -sfile=$f -bednames=$samplelist -frag="NA" -bkgd=$bkgdfile -out="$wdir/peaks.$ran/" -clustnummer="$clustno" -resolution="$resol" -window="$window" -p="$cores" -bin="$binsize" -type="$type" -chrcount="$counting" -initModel="$initModel" -windowe="$windowe"
counting=$(($counting+1));
cp "$wdir/peaks.$ran/$chr.peaks.bed" "$out/peaks/$chr.peaks.bed"
rm "$wdir/peaks.$ran/$chr.peaks.bed"
@@ -447,10 +495,12 @@
#concatenate, sort and filter
-cat $out/peaks/*.bed > $out/peaks/all.narrowPeak
-Rscript "$sPath/peakhelper.r" -filelist="$out/peaks/all.narrowPeak"
-perl "$sPath/peakfilter.pl" $out/peaks/all.narrowPeak | sort -nr -k7 > $out/peaks/filtered.peaks.narrowPeak
-cut -f1-10 $out/peaks/all.narrowPeak | awk -F"\t" -v j=0 '$7 > j' | sort -nr -k7 > $out/peaks/all.peaks.narrowPeak
+cat $out/peaks/*.bed > $out/peaks/all.narrowPeak
+if [[ -s $out/peaks/all.narrowPeak ]]; then
+ Rscript "$sPath/peakhelper.r" -filelist="$out/peaks/all.narrowPeak"
+ perl "$sPath/peakfilter.pl" $out/peaks/all.narrowPeak | sort -nr -k7 > $out/peaks/filtered.peaks.narrowPeak
+ cut -f1-10 $out/peaks/all.narrowPeak | awk -F"\t" -v j=0 '$7 > j' | sort -nr -k7 > $out/peaks/all.peaks.narrowPeak
+fi
rm $out/peaks/all.narrowPeak
rm $out/peaks/*.bed
rm $out/peaks/min.peaksize
diff -r d42f4d78c85e -r 243f75d0ed6e README.md
--- a/README.md Wed Dec 17 10:40:23 2014 -0500
+++ b/README.md Thu Feb 19 05:39:45 2015 -0500
@@ -7,6 +7,15 @@
determine enrichment site widths accurately and resolve neighboring narrow peaks.
JAMM is a universal peak finder that is applicable to different types of datasets
+Stand-alone version
+===================
+
+The latest release of the stand-alone version can be found here:
+http://code.google.com/p/jamm-peak-finder/
+
+(You can also use the stand-alone version from this repository, ignoring the
+Galaxy-specific files. The upstream version might be never, though)
+
Galaxy
====
diff -r d42f4d78c85e -r 243f75d0ed6e SignalGenerator.sh
--- a/SignalGenerator.sh Wed Dec 17 10:40:23 2014 -0500
+++ b/SignalGenerator.sh Thu Feb 19 05:39:45 2015 -0500
@@ -1,7 +1,22 @@
-#!/usr/bin/bash
-
-
-
+########################################################################
+# JAMMv1.0.7rev1 is a peak finder for joint analysis of NGS replicates.
+# Copyright (C) 2014-2015 Mahmoud Ibrahim
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# Contact: mahmoud.ibrahim@mdc-berlin.de
+########################################################################
##Finding out the path
sPath="`dirname \"$0\"`"
@@ -12,7 +27,7 @@
usage()
{
cat << EOF
-Welcome to JAMM v1.0.6rev2 Signal Generator Script (GNU GPLv3). Copyright (C) 2014 Mahmoud Ibrahim.
+Welcome to JAMM v1.0.7rev1 Signal Generator Script (GNU GPLv3). Copyright (C) 2014-2015 Mahmoud Ibrahim.
This program comes with ABSOLUTELY NO WARRANTY; for details visit http://www.gnu.org/licenses/gpl.html. This is free software, and you are welcome to redistribute it under certain conditions; visit http://www.gnu.org/licenses/gpl.html for details.
@@ -23,8 +38,9 @@
-c directory containing input or Control files
-r file with Regions to get signal for (required)
-b Bin size for signal generation (default: 10)
- -f Fragment lengths (required)
+ -f Fragment lengths (required if -t is "single")
-p Number of processors used by R scripts (default: 1)
+ -t Alignment type, paired or single (default: single)
EOF
}
@@ -44,9 +60,9 @@
wdir=$(mktemp -d)
ran=$RANDOM
cores="1"
-
+type="single"
-while getopts "s:o:c:r:b:f:g:p:" OPTION
+while getopts "s:o:c:r:b:f:g:p:t:" OPTION
do
case $OPTION in
s) sdir=$OPTARG
@@ -65,6 +81,8 @@
;;
p) cores=$OPTARG
;;
+ t) type=$OPTARG
+ ;;
?)
usage
exit
@@ -73,11 +91,24 @@
done
-if [[ -z $fraglen ]] || [[ -z $sdir ]] || [[ -z $regs ]] || [[ -z $out ]]
+if [[ -z $sdir ]] || [[ -z $regs ]] || [[ -z $out ]]
then
usage
exit 1
fi
+
+if [[ -z $fraglen ]]
+then
+ if [ $type == "single" ]
+ then
+ usage
+ exit 1
+ fi
+fi
+if [[ -d "$out/signal" ]]; then
+ printf "\n\nOutput directory $out/signal already exists. I can't override existing results!\n\n"
+ exit 0
+fi
#=======================> DONE!
@@ -85,7 +116,7 @@
# =============================
# Step One: Initial Processing
# =============================
-printf "\n\n============================================\nStarted JAMM Signal Generator Pipeline...Hang on!\n============================================\n\n"
+printf "\n\n=======================================================\nStarted JAMM1.0.7rev1 Signal Generator Pipeline...Hang on!\n=======================================================\n\n"
if [ ! -d "$wdir" ]; then
mkdir $wdir #make working directory
@@ -131,13 +162,13 @@
awk -v chr="$chr" -v ext="$wdir/bkgd.$ran/" -F"\t" '$1 == chr { print $2"\t"$6 >> ext"bkgd."chr".ctrl.bed" }' "$wdir/bkgd.$ran/ctrl.bed"
done
-rm $wdir/bkgd.$ran/ctrl.bed
printf "Done!\n"
fi
#determine average read lengths
+if [ $type == "single" ]; then
printf "Getting average read lengths..."
if [ ! -z $bdir ]; then
readC=$(awk '{a=$3-$2;print a;}' "$wdir/bkgd.$ran/ctrl.bed" | perl -lane '$a+=$_;END{print $a/$.}' | awk '{print int($1)}')
@@ -149,6 +180,7 @@
done
readL=${readL#","}
printf "Done!\n"
+fi
#=======================> DONE!
@@ -165,6 +197,7 @@
printf "Generating Signal...(bin size: $signal)\n"
+if [ $type == "single" ]; then
counting=1;
for f in $wdir/sizes.$ran/*; do #for each chromosome
samplelist=""
@@ -207,12 +240,57 @@
frag=${frag#","}
#call the peak calling R script
- Rscript "$sPath/signalmaker.r" -chromoS="$chrSize" -chromo="$chr" -bednames=$samplelist -frag=$frag -bkgd=$bkgdfile -out="$wdir/signal.$ran/" -sig="$signal" -regions="$regs" -p="$cores" -chrcount="$counting"
+ Rscript "$sPath/signalmaker.r" -chromoS="$chrSize" -chromo="$chr" -bednames=$samplelist -frag=$frag -bkgd=$bkgdfile -out="$wdir/signal.$ran/" -sig="$signal" -regions="$regs" -p="$cores" -chrcount="$counting" -t="$type"
cat $wdir/signal.$ran/*.bedSignal | awk -F"\t" -v j=0 '$4 > j' > "$out/signal/$chr.bedGraph"
rm $wdir/signal.$ran/*.bedSignal
counting=$(($counting+1));
done
counting=1;
+fi
+
+
+
+if [ $type == "paired" ]; then
+counting=1;
+for f in $wdir/sizes.$ran/*; do #for each chromosome
+ samplelist=""
+
+ sizefile=$(basename $f)
+ chr=$(echo $sizefile | awk -F"." '{print $2}' | awk -F"." '{print $1}');
+ chrSize=$(cat $f | cut -f2);
+ printf "Chromosome $chr: "
+
+ #list of sample bed files and fragment lengths
+ for s in $wdir/samples.$ran/*.bed; do #and for each sample file
+ samplefile=$(basename $s)
+ chr2=$(echo $samplefile | awk -F"." '{print $2}');
+ if [ $chr == $chr2 ] #belonging to this chromosome
+ then
+ samplelist="$samplelist,$wdir/samples.$ran/$samplefile"
+ samplename=$(echo $samplefile | awk -F"." '{ print $3 }')
+ samplefilename=$(echo $samplefile | cut -d'.' -f 3-)
+ x="$sdir/$samplefilename"
+ fi
+ done
+
+ #control file
+ bkgdfile="None"
+ if [ ! -z $bdir ]; then
+ bkgdfile="$wdir/bkgd.$ran/bkgd.$chr.ctrl.bed"
+ fi
+
+ #remove leading comma
+ samplelist=${samplelist#","}
+ frag=${frag#","}
+
+ #call the peak calling R script
+ Rscript "$sPath/signalmaker.r" -chromoS="$chrSize" -chromo="$chr" -bednames=$samplelist -frag=$frag -bkgd=$bkgdfile -out="$wdir/signal.$ran/" -sig="$signal" -regions="$regs" -p="$cores" -chrcount="$counting" -t="$type"
+ cat $wdir/signal.$ran/*.bedSignal | awk -F"\t" -v j=0 '$4 > j' > "$out/signal/$chr.bedGraph"
+ rm $wdir/signal.$ran/*.bedSignal
+ counting=$(($counting+1));
+done
+counting=1;
+fi
#=======================> DONE!
rm -rf $wdir
diff -r d42f4d78c85e -r 243f75d0ed6e bincalculator.r
--- a/bincalculator.r Wed Dec 17 10:40:23 2014 -0500
+++ b/bincalculator.r Thu Feb 19 05:39:45 2015 -0500
@@ -1,7 +1,22 @@
-#############################################################
-#### Bin Size Calculator after Shimazaki and Shinomoto 2007
-#### R script
-#############################################################
+########################################################################
+# JAMMv1.0.7rev1 is a peak finder for joint analysis of NGS replicates.
+# Copyright (C) 2014-2015 Mahmoud Ibrahim
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# Contact: mahmoud.ibrahim@mdc-berlin.de
+########################################################################
# =======================
@@ -171,9 +186,9 @@
-# =====================
-# Shimazaki Procedure
-# =====================
+# ===================================================
+# Shimazaki Procedure (Shimazaki and Shinomoto 2007)
+# ===================================================
for (i in 1:length(ibed)) {
if (frags[i] > rl[i]) {
minbin = floor(frags[i] / 2)
diff -r d42f4d78c85e -r 243f75d0ed6e changelog
--- a/changelog Wed Dec 17 10:40:23 2014 -0500
+++ b/changelog Thu Feb 19 05:39:45 2015 -0500
@@ -1,6 +1,6 @@
===================
-Version 1.0.6rev2
+Version 1.0.7rev1
===================
-- fixed various bugs in v1.0.6rev1 that caused paired-end runs to malfunction.
-
+- default seed is 1011414
+-
diff -r d42f4d78c85e -r 243f75d0ed6e jamm.xml
--- a/jamm.xml Wed Dec 17 10:40:23 2014 -0500
+++ b/jamm.xml Thu Feb 19 05:39:45 2015 -0500
@@ -1,20 +1,28 @@
-
+Calls peaks on NGS data
-
+ R
+ perl
-jammwrapper.py -i
- #for $x in $replicate
- ${x.input.file_name}
- #end for
- -g ${genomeSizeFile}
- -z $output
+
@@ -24,12 +32,15 @@
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
@@ -37,19 +48,27 @@
+
-
-
+
+
+
+
+
+
+
+
+
-
-
+
+
+
+
+
+
+
+
+
+
@@ -64,13 +95,23 @@
10.1093/bioinformatics/btu568
-.. class:: infomark
+**What it does**
+JAMM calls peaks on NGS data, i.e. for finding binding sites.
+
+.. class:: warningmark
This tool requires files in *bed* format.
-----
+.. class:: infomark
+
+If you are just starting out, you can find some sensible parameter choices here:
+
+http://code.google.com/p/jamm-peak-finder/wiki/Usage#Examples
-**What it does**
-JAMM calls peaks on NGS data, i.e. for finding binding sites.
+* ChIP-Seq TF: everything default
+* ChIP-Seq HM: Binsize=200
+* ChIP-Seq Broad HM: Binsize=200 Resolution=region
+* ChIP-Seq Very Broad HM: Binsize=10000 Resolution=window
+* DNase-seq: Fragment length=1
diff -r d42f4d78c85e -r 243f75d0ed6e jammwrapper.py
--- a/jammwrapper.py Wed Dec 17 10:40:23 2014 -0500
+++ b/jammwrapper.py Thu Feb 19 05:39:45 2015 -0500
@@ -6,7 +6,7 @@
# map the files (from Galaxy history) to a directory
# pass the parameters from the GUI to JAMM.sh
# call JAMM.sh
-# maybe map the resulting tabular files back to history items, the XML could take care of that if output is fixed
+# map the resulting tabular files back to history items
#import optparse
import argparse, os, shutil, subprocess, sys, tempfile
@@ -22,48 +22,77 @@
#Read command line arguments
parser = argparse.ArgumentParser()
- parser.add_argument('-i', dest = 'input', nargs='+')
+ parser.add_argument('-i', dest = 'input', nargs='+') # + allows for 1 or more arguments
+ # parser.add_argument('-c', dest = 'control', nargs='*', default=None) # allow for not supplying -c
+ parser.add_argument('-c', dest = 'control', action='append') # allow for not supplying -c
parser.add_argument('-g', dest = 'gsize')
- parser.add_argument('-z', dest = 'peakfile')
+ parser.add_argument('-o', dest = 'peakfile')
+ parser.add_argument('-of', dest = 'filteredpeakfile')
parser.add_argument('-m', dest = 'mode')
parser.add_argument('-r', dest = 'resolution')
parser.add_argument('-p', dest = 'processes')
- #parser.add_argument('-b', dest = 'binSize')
+ parser.add_argument('-t', dest = 'type')
+ parser.add_argument('-f', dest = 'fraglen')
+ parser.add_argument('-b', dest = 'binsize')
args = parser.parse_args()
- print "\n"
+ print "################################################"
print "Wrapper debugging"
- print "Files to be used:"
+ print "################################################"
+ print "Sample files:"
for j in args.input:
print j
+ if args.control is not None:
+ print "Control files:"
+ for j in args.control:
+ print j
- print "output file:"
+ print "output files:"
print args.peakfile
+ print args.filteredpeakfile
+ print "current working dir:"
+ print os.getcwd()
+ print "dir with jammwrapper in it:"
+ path = os.path.abspath(os.path.dirname(sys.argv[0]))
+ print path
- # depracted, can still be found in may of the example wrappers
+ # optparse was depracted, can still be found in may of the example wrappers
#parser = optparse.OptionParser()
#parser.add_option( '-i', dest='input', help='input bed files' )
#parser.add_option( '-c', dest='csize', help='chr sizes' )
#(options, args) = parser.parse_args()
-
+
# create temp dir
tmp_dir = tempfile.mkdtemp()
- os.chdir(tmp_dir)
+ os.mkdir(tmp_dir + "/sample")
# symlink creation
for file in args.input:
- filen = tmp_dir + "/" + os.path.basename(os.path.splitext(file)[0])+".bed"
+ filen = tmp_dir + "/sample/" + os.path.basename(os.path.splitext(file)[0])+".bed"
+ print "input files mapped: %s" % filen
os.symlink(file, filen)
-
- # in case temp files should have random names
- # ref_file = tempfile.NamedTemporaryFile( dir=tmp_dir )
- # ref_file_name = ref_file.name
- # ref_file.close()
- # os.symlink( file, ref_file_name )
-
- command = ( "/home/cmesser/galaxy/tools/jamm/JAMM.sh -s %s -g %s -o results -m %s -r %s -p %s"
- % ( tmp_dir, args.gsize, args.mode, args.resolution, args.processes ) )
+
+ # Here comes some unnecessary repetition
+ # if control files are supplied, we make another tmp subdir and put those files there.
+ # JAMM.sh is then called with one additional switch -c
+ if args.control is not None:
+ # symlink creation
+ os.mkdir(tmp_dir + "/control")
+ for file in args.control:
+ filen = tmp_dir + "/control/" + os.path.basename(os.path.splitext(file)[0])+".bed"
+ print "input files mapped: %s" % filen
+ os.symlink(file, filen)
+ command = ( "bash %s/JAMM.sh -s %s/sample -c %s/control -g %s -o results -m %s -r %s -p %s -t %s -f %s -b %s"
+ % ( path, tmp_dir, tmp_dir, args.gsize, args.mode, args.resolution, args.processes, \
+ args.type, args.fraglen, args.binsize ) )
+ else:
+ command = ( "bash %s/JAMM.sh -s %s/sample -g %s -o results -m %s -r %s -p %s -t %s -f %s -b %s"
+ % ( path, tmp_dir, args.gsize, args.mode, args.resolution, args.processes, \
+ args.type, args.fraglen, args.binsize ) )
+
+
+ print "Command called by bash:"
print command
# depending on how your programm is called, it may be necessary to use shlex.split on the command string before
# in this case, this was actually harmful. idk why
@@ -75,13 +104,15 @@
returncode = proc.wait()
#mv files to a place where galaxy wrapper can find them
- mvcommand = "mv %s/results/peaks/all.peaks.narrowPeak %s" % ( tmp_dir, args.peakfile )
+ # mvcommand = "mv %s/results/peaks/all.peaks.narrowPeak %s" % ( tmp_dir, args.peakfile )
+ mvcommand = "mv results/peaks/all.peaks.narrowPeak %s" % args.peakfile
+ os.system(mvcommand)
+ mvcommand = "mv results/peaks/filtered.peaks.narrowPeak %s" % args.filteredpeakfile
os.system(mvcommand)
# clean up temp dir
if os.path.exists( tmp_dir ):
shutil.rmtree( tmp_dir )
-
if __name__ == "__main__":
main()
diff -r d42f4d78c85e -r 243f75d0ed6e peakfilter.pl
--- a/peakfilter.pl Wed Dec 17 10:40:23 2014 -0500
+++ b/peakfilter.pl Thu Feb 19 05:39:45 2015 -0500
@@ -1,8 +1,22 @@
-#!/usr/bin/perl
-########################################
-#### Filters peaks by p-value and width
-#### Perl script
-#######################################
+########################################################################
+# JAMMv1.0.7rev1 is a peak finder for joint analysis of NGS replicates.
+# Copyright (C) 2014-2015 Mahmoud Ibrahim
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# Contact: mahmoud.ibrahim@mdc-berlin.de
+########################################################################
use feature qw(say);
diff -r d42f4d78c85e -r 243f75d0ed6e peakfinder.r
--- a/peakfinder.r Wed Dec 17 10:40:23 2014 -0500
+++ b/peakfinder.r Thu Feb 19 05:39:45 2015 -0500
@@ -1,16 +1,35 @@
-###############################
-#### Peak finder for NGS Data
-#### R script
-###############################
+########################################################################
+# JAMMv1.0.7rev1 is a peak finder for joint analysis of NGS replicates.
+# Copyright (C) 2014-2015 Mahmoud Ibrahim
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# Contact: mahmoud.ibrahim@mdc-berlin.de
+########################################################################
+
+
# =======================
# User-defined variables
# =======================
-cutoff = NA #To enforce an SNR cutoff ratio, delete NA and enter the number you want
-strict = 1 #to make peak calling more / less strict, increase / decrease this number
-meanAdjust = "n" #adjust the initialization mean vector for each window before clustering? If you want to do so, change this to "y".
-options(warn = -1, scipen = 1000) #R will not report any warnings (warn = -1), R will not use scientific notation (scipen = 1000)
+samplingSeed = 1011414 #the seed makes the random sampling that JAMM does "deterministic". You can change this to your integer seed of choice to reproduce exact same results or randomize it this way: ceiling(rnorm(1,50,300)).
+reportNoClust = "n" #report windows for which clustering failed? Reported windows will be marked by "NoClust" in the peak name (4th column).
+cutoff = NA #To enforce an SNR cutoff ratio for bin enrichment calling, delete NA and enter the number you want.
+strict = 1 #To make bin enrichment calling more / less strict, increase / decrease this number.
+meanAdjust = "n" #Adjust the initialization mean vector for each window before clustering? If you want to do so, change this to "y".
+options(warn = -1, scipen = 1000) #R will not report any warnings (warn = -1), R will not use scientific notation (scipen = 1000).
#=======================> DONE!
@@ -115,7 +134,7 @@
}
counts = counts/mCount
}
-
+
return(counts)
}
@@ -146,7 +165,7 @@
}
counts = counts/mCount
}
-
+
return(counts)
}
@@ -168,11 +187,8 @@
Cs = counts[[numdup+1]][winStart:tempend]
mCs = mean(Cs)
dCs = sd(Cs)
- } else {
- Cs = sample(C, binSizeTemp, replace = TRUE)
}
-
go = rep(0, numdup)
for (g in 1:numdup) {
mS = (mean(counts[[g]][winStart:tempend]))
@@ -236,6 +252,7 @@
Cs = Cs[winStart:winEnd]
Cs = filtfilt(rep(1,80)/80,1,Cs) + Cmin #gets rid of Inf in the fold change
} else {
+ set.seed(samplingSeed)
Cs = sample(C, rWinSizeTemp, replace = TRUE)
Cs = filtfilt(rep(1,80)/80,1,Cs)
}
@@ -265,6 +282,7 @@
return(Rs)
}
cluster = function(model, sig, init, clustnummer, noise) { #helper function2
+ set.seed(samplingSeed)
noisy = sample(noise, length(sig[,1]), replace = TRUE)
clust = me(model, sig+noisy, init)
bicc = bic(model, clust$loglik, length(sig[,1]), length(sig[1,]), clustnummer)
@@ -283,7 +301,8 @@
sig = do.call(rbind, sig)
#kmeans initialization
- init = kmeans(sig, clustnummer)
+ set.seed(samplingSeed)
+ init = kmeans(sig, clustnummer, nstart = 20)
init = unmap(init$cluster)
@@ -307,15 +326,16 @@
#find peaks
-findpeak = function(winStart, coffeeshopSud, numdup, C, param, bkgd, resol, counts, noise, startlist, meanAdjust) {
+findpeak = function(winStart, coffeeshopSud, numdup, C, param, bkgd, resol, counts, noise, startlist, meanAdjust, clustnummer) {
-
+
plz = which(startlist == winStart)
winEnd = coffeeshopSud[plz]
#will store peak information
writethis = list()
+ ccx = 1 #default is clustering didNOT work
rWinSizeTemp = winEnd - winStart + 1
@@ -338,20 +358,32 @@
}
if (resol != "window") {
-
- ####CLUSTERING START HERE#####
+
+ #clustering (take 1)
+ take = 1
+ set.seed(samplingSeed)
noisy = sample(noise, rWinSizeTemp, replace = TRUE)
clust = em(param$modelname, Rs+noisy, param$initparam)
clust$classification = map(clust$z)
- ####CLUSTERING DONE####
-
- #check whether clutering succeeded
- if((any(diff(clust$classification)) != 0) && (!(any(is.na(clust$classification))))) {
-
-
- #check whether all components replicates agreed on clustering assignments
- ccx = 1
- if (numdup > 1) {
+ if (!((any(diff(clust$classification)) != 0) && (!(any(is.na(clust$classification)))))) { #clustering didn't work, take1
+
+ #repeat clustering from scratch, take 2!
+ set.seed(samplingSeed)
+ init = kmeans(Rs, clustnummer, nstart = 20)
+ init = unmap(init$cluster)
+ set.seed(samplingSeed)
+ noisy = sample(noise, rWinSizeTemp, replace = TRUE)
+ clust = me(param$modelname, Rs+noisy, init)
+ clust$classification = map(clust$z)
+ if ((any(diff(clust$classification)) != 0) && (!(any(is.na(clust$classification))))) {
+ ccx = 0 #clustering worked, take2
+ take = 2
+ }
+ } else {ccx = 0} #clustering worked, take1
+
+ if (ccx != 1) { #clustering worked either in take1 or take2
+
+ if (numdup > 1) { #check whether all components replicates agreed on clustering assignments
cc = vector(mode = "numeric", length = numdup)
for (g in 1:numdup) {
cc[g] = which.max(clust$parameters$mean[g,]) #which cluster has the largest mean (this is the peak cluster, hopefully!)
@@ -359,98 +391,123 @@
ccx = sum(diff(cc))
cluster = cc[1]
rm(cc)
- } else {
- cluster = which.max(clust$parameters$mean) #which cluster has the largest mean (this is the peak cluster, hopefully!)
- ccx = 0
- }
- if (ccx == 0) {
- #extract subset of the background
- if (bkgd != "None") {
- Cs = counts[[numdup+1]][winStart:winEnd]
- Cs = filtfilt(rep(1,80)/80,1,Cs)
- } else {
- Cs = sample(C, rWinSizeTemp, replace = TRUE)
- Cs = filtfilt(rep(1,80)/80,1,Cs)
+ if ((ccx != 0) && (take == 1)) { #not all replicates agreed? Repeat the clustering with from scratch if not already done!
+ set.seed(samplingSeed)
+ init = kmeans(Rs, clustnummer, nstart = 20)
+ init = unmap(init$cluster)
+ set.seed(samplingSeed)
+ noisy = sample(noise, rWinSizeTemp, replace = TRUE)
+ clust = me(param$modelname, Rs+noisy, init)
+ clust$classification = map(clust$z)
+ if ((any(diff(clust$classification)) != 0) && (!(any(is.na(clust$classification))))) { #clustering worked? check whether replicates agreed take 3
+ cc = vector(mode = "numeric", length = numdup)
+ for (g in 1:numdup) {
+ cc[g] = which.max(clust$parameters$mean[g,]) #which cluster has the largest mean (this is the peak cluster, hopefully!)
+ }
+ ccx = sum(diff(cc))
+ cluster = cc[1]
+ rm(cc)
+ take = 3
+ }
}
+ } else { #no replicates!
+ cluster = which.max(clust$parameters$mean) #which cluster has the largest mean (this is the peak cluster, hopefully!)
+ }
+ }
+
+ if ((ccx != 0) && (reportNoClust=="y")) { resol = "window" } #clustering did not work and windows should be reported
+
+
+ if (ccx == 0) { #clustering worked and all replicates agree on the cluster assignments
+
+ #extract subset of the background
+ if (bkgd != "None") {
+ Cs = counts[[numdup+1]][winStart:winEnd]
+ Cs = filtfilt(rep(1,80)/80,1,Cs)
+ } else {
+ set.seed(samplingSeed)
+ Cs = sample(C, rWinSizeTemp, replace = TRUE)
+ Cs = filtfilt(rep(1,80)/80,1,Cs)
+ }
- #find region boundaries
- loc = 1:length(clust$classification)
- gmclass = cbind(loc, clust$classification)
- locPeak = gmclass[gmclass[,2] == cluster,,drop=FALSE]
- rStart = locPeak[1] #start position of the region
- rEnd = locPeak[length(locPeak[,1]),1] #end position of the region
+ #find region boundaries
+ loc = 1:length(clust$classification)
+ gmclass = cbind(loc, clust$classification)
+ locPeak = gmclass[gmclass[,2] == cluster,,drop=FALSE]
+ rStart = locPeak[1] #start position of the region
+ rEnd = locPeak[length(locPeak[,1]),1] #end position of the region
- #peak resolution check
- if (resol == "region") {
- pSize = rEnd - rStart
- signal = (geomean(Rs[rStart:rEnd,]))
- signal2 = (signal) - (Cs[rStart:rEnd])
+ #peak resolution check
+ if (resol == "region") {
+ pSize = rEnd - rStart
+ signal = (geomean(Rs[rStart:rEnd,]))
+ signal2 = (signal) - (Cs[rStart:rEnd])
+ gm = mean(signal2)
+ summit = which.max(geomean(Rsr[rStart:rEnd,])) - 1
+ will2k = wilcox.test(signal, Cs[rStart:rEnd])
+
+ #Is there signal in the region above background
+ if (gm > 0) {
+ writethis[[1]] = rStart + winStart - 1
+ writethis[[2]] = rEnd + winStart
+ writethis[[3]] = paste0(chromName, ".", rStart+winStart -1)
+ writethis[[4]] = "1000"
+ writethis[[5]] = "."
+ writethis[[6]] = gm
+ writethis[[7]] = will2k$p.value
+ writethis[[8]] = "-1"
+ writethis[[9]] = summit
+ }
+ } else if (resol == "peak") {
+ #find out where separate peaks are
+ d = diff(locPeak[,1])
+ d[length(d)+1] = 0
+ locPeak = cbind(locPeak, d)
+ bound1 = which(locPeak[,3] > 1, arr.in=TRUE)
+ bound2 = bound1 + 1
+ bound = locPeak[sort(c(bound1,bound2))]
+ bound = c(rStart, bound, rEnd)
+ w = 1
+ warum = 0
+ while (w < length(bound)) {
+ pStart = bound[w] + winStart - 1
+ pEnd = bound[w+1] + winStart
+ pSize = pEnd - pStart
+ signal = (geomean(Rs[(bound[w]):(bound[w+1]),]))
+ signal2 = (signal) - (Cs[bound[w]:bound[w+1]])
gm = mean(signal2)
- summit = which.max(geomean(Rsr[rStart:rEnd,])) - 1
- will2k = wilcox.test(signal, Cs[rStart:rEnd])
-
+ summit = which.max(geomean(Rsr[(bound[w]):(bound[w+1]),])) - 1
+ will2k = wilcox.test(signal, Cs[(bound[w]):(bound[w+1])])
+
+ weil = warum * 9
#Is there signal in the region above background
if (gm > 0) {
- writethis[[1]] = rStart + winStart - 1
- writethis[[2]] = rEnd + winStart
- writethis[[3]] = paste0(chromName, ".", rStart+winStart -1)
- writethis[[4]] = "1000"
- writethis[[5]] = "."
- writethis[[6]] = gm
- writethis[[7]] = will2k$p.value
- writethis[[8]] = "-1"
- writethis[[9]] = summit
+ writethis[[1+weil]] = pStart
+ writethis[[2+weil]] = pEnd
+ writethis[[3+weil]] = paste0(chromName, ".", pStart)
+ writethis[[4+weil]] = "1000"
+ writethis[[5+weil]] = "."
+ writethis[[6+weil]] = gm
+ writethis[[7+weil]] = will2k$p.value
+ writethis[[8+weil]] = "-1"
+ writethis[[9+weil]] = summit
}
- } else if (resol == "peak") {
- #find out where separate peaks are
- d = diff(locPeak[,1])
- d[length(d)+1] = 0
- locPeak = cbind(locPeak, d)
- bound1 = which(locPeak[,3] > 1, arr.in=TRUE)
- bound2 = bound1 + 1
- bound = locPeak[sort(c(bound1,bound2))]
- bound = c(rStart, bound, rEnd)
- w = 1
- warum = 0
- while (w < length(bound)) {
- pStart = bound[w] + winStart - 1
- pEnd = bound[w+1] + winStart
- pSize = pEnd - pStart
- signal = (geomean(Rs[(bound[w]):(bound[w+1]),]))
- signal2 = (signal) - (Cs[bound[w]:bound[w+1]])
- gm = mean(signal2)
- summit = which.max(geomean(Rsr[(bound[w]):(bound[w+1]),])) - 1
- will2k = wilcox.test(signal, Cs[(bound[w]):(bound[w+1])])
-
- weil = warum * 9
- #Is there signal in the region above background
- if (gm > 0) {
- writethis[[1+weil]] = pStart
- writethis[[2+weil]] = pEnd
- writethis[[3+weil]] = paste0(chromName, ".", pStart)
- writethis[[4+weil]] = "1000"
- writethis[[5+weil]] = "."
- writethis[[6+weil]] = gm
- writethis[[7+weil]] = will2k$p.value
- writethis[[8+weil]] = "-1"
- writethis[[9+weil]] = summit
- }
- w = w + 2
- warum = warum + 1
- }
-
- } #peak resolution check
- } #all replicates agree on clustering assignments?
- } #clustering worked?
-
- } else if (resol == "window") { #window resolution check
+ w = w + 2
+ warum = warum + 1
+ }
+ } #peak resolution check
+ } #clustering worked and all replicates agree on clustering assignments?
+ } #window resolution check
+
+ if (resol == "window") {
#extract subset of the background
if (bkgd != "None") {
Cs = counts[[numdup+1]][winStart:winEnd]
Cs = filtfilt(rep(1,80)/80,1,Cs)
} else {
+ set.seed(samplingSeed)
Cs = sample(C, rWinSizeTemp, replace = TRUE)
Cs = filtfilt(rep(1,80)/80,1,Cs)
}
@@ -467,7 +524,7 @@
if (gm > 0) {
writethis[[1]] = winStart - 1
writethis[[2]] = winEnd
- writethis[[3]] = paste0(chromName, ".", winStart -1)
+ writethis[[3]] = paste0(chromName, ".", winStart -1, ".NoClust")
writethis[[4]] = "1000"
writethis[[5]] = "."
writethis[[6]] = gm
@@ -475,7 +532,7 @@
writethis[[8]] = "-1"
writethis[[9]] = summit
}
- } #window clustering
+ } #window reporting
return(writethis)
}
@@ -506,7 +563,7 @@
#Parsing arguments and storing values
for (each.arg in args) {
#chormosome size file
- if (grepl('^-sfile=',each.arg)) {
+ if (grepl('-sfile=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
size.file <- arg.split[2]
@@ -515,7 +572,7 @@
}
}
#bed file names
- if (grepl('^-bednames=',each.arg)) {
+ if (grepl('-bednames=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
bednames <- arg.split[2]
@@ -524,7 +581,7 @@
}
}
#bed files directory
- if (grepl('^-frag=',each.arg)) {
+ if (grepl('-frag=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
frag <- arg.split[2]
@@ -533,7 +590,7 @@
}
}
#bakcground files directory
- if (grepl('^-bkgd=',each.arg)) {
+ if (grepl('-bkgd=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
bkgd <- arg.split[2]
@@ -542,7 +599,7 @@
}
}
#bakcground files directory
- if (grepl('^-out=',each.arg)) {
+ if (grepl('-out=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
out <- arg.split[2]
@@ -551,61 +608,69 @@
}
}
#Cluster number
- if (grepl('^-clustnummer=',each.arg)) {
+ if (grepl('-clustnummer=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
clustnummer <- as.numeric(arg.split[2])
}
}
#resolution
- if (grepl('^-resolution=',each.arg)) {
+ if (grepl('-resolution=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
resol <- arg.split[2]
}
}
#processor cores
- if (grepl('^-p=',each.arg)) {
+ if (grepl('-p=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
cornum <- as.numeric(arg.split[2])
}
}
#minimum window size
- if (grepl('^-window=',each.arg)) {
+ if (grepl('-window=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
winSize <- arg.split[2]
}
}
#window size
- if (grepl('^-bin=',each.arg)) {
+ if (grepl('-bin=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
binsize <- arg.split[2]
}
}
#type (paired / single)
- if (grepl('^-type=',each.arg)) {
+ if (grepl('-type=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
type <- arg.split[2]
}
}
#chromosome number
- if (grepl('^-chrcount=',each.arg)) {
+ if (grepl('-chrcount=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
chrcount <- as.numeric(arg.split[2])
}
}
#window enrichment cutoff
- if (grepl('^-windowe=',each.arg)) {
+ if (grepl('-windowe=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
windowe <- arg.split[2]
}
}
+ #initialize
+ if (grepl('-initModel=',each.arg)) {
+ arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
+ if (! is.na(arg.split[2]) ) {
+ initialize <- arg.split[2]
+ }
+ }
+
}
##Parse in variables
@@ -614,6 +679,11 @@
chromSize = chromosomes$V2; #chromosomes size
rm(chromosomes)
+if (chrcount == 1) {
+ write(paste(samplingSeed), file = paste0(out, "/seed.info"))
+} else {
+ samplingSeed = as.numeric(read.table(paste0(out, "/seed.info"), header = FALSE)[[1]])
+}
readsFiles = as.list(strsplit(bednames, ",", fixed = TRUE)[[1]])
numdup = length(readsFiles) #number of replicates
if (bkgd != "None") {
@@ -703,7 +773,6 @@
# ============================
# Estimating Background Model
# ============================
-
if (chrcount == 1){ #first chromosome, estimate bkgd (includes SNR cutoff)
if (is.na(cutoff)) {
if (bkgd != "None") {
@@ -722,13 +791,16 @@
mmM = mean(geomeanL(counts))
sigma = log(1+((mmV) / ((mmM)^2)))
mu = (log(mmM)) - (0.5 * (sigma))
+ set.seed(samplingSeed)
C = rlnorm(100000, mu, sqrt(sigma))
for (x in 1:numdup) {
cutoff[x] = (mean(counts[[x]]))/(sd(C))
}
cutoff = max(cutoff)
- mCs = (mean(sample(C, binSize*5, replace = TRUE)))
- dCs = (sd(sample(C, binSize*5, replace = TRUE)))
+ set.seed(samplingSeed)
+ snow = sample(C, binSize*5, replace = TRUE)
+ mCs = mean(snow)
+ dCs = sd(snow)
write(paste(c(cutoff,sigma,mu)), file = paste0(out, "/bkgd.info"), append = TRUE)
}
}
@@ -744,9 +816,12 @@
} else {
sigma = as.numeric(bkgdInfo[[1]][2])
mu = as.numeric(bkgdInfo[[1]][3])
+ set.seed(samplingSeed)
C = rlnorm(100000, mu, sqrt(sigma))
- mCs = (mean(sample(C, binSize*5, replace = TRUE)))
- dCs = (sd(sample(C, binSize*5, replace = TRUE)))
+ set.seed(samplingSeed)
+ snow = sample(C, binSize*5, replace = TRUE)
+ mCs = mean(snow)
+ dCs = sd(snow)
}
}
#=======================> DONE!
@@ -764,6 +839,9 @@
}
coffeeshop = as.numeric(unlist(coffeeshop))
coffeeshop[coffeeshop != numdup] = 0
+
+if (sum(coffeeshop) != 0) { #Any enriched bins?
+
coffeeshop = c(0, diff(coffeeshop))
coffeeshop = cbind(coffeeshop, bins)
coffeeshopNord = coffeeshop[coffeeshop[,1] == numdup,,drop=FALSE]
@@ -800,16 +878,26 @@
# ===================================
if (length(coffeeshop[,1]) > 0) { #any enriched windows detected?
-yummy = ceiling(length(coffeeshop[,1]) / 4)
-if (yummy > 20) {
- yummy = sample(1:yummy, 20)
-} else if (yummy > 0) {
- yummy = 1:yummy
-} else {
- yummy = 1
+if (initialize == "deterministic") {
+ yummy = ceiling(length(coffeeshop[,1]) / 1000)
+ if (yummy == 0) {
+ yummy = 1
+ }
+}
+if (initialize == "stochastic") {
+ yummy = ceiling(length(coffeeshop[,1]) / 4)
+ if (yummy > 20) {
+ set.seed(samplingSeed)
+ yummy = sample(1:yummy, 20)
+ } else if (yummy > 0) {
+ yummy = 1:yummy
+ } else {
+ yummy = 1
+ }
}
coffeeshopNord = coffeeshop[yummy,1]
coffeeshopSud = coffeeshop[yummy,2]
+set.seed(samplingSeed)
noise = rnorm(100000, mean=0, sd=0.1)
param = initparam(coffeeshopNord, coffeeshopSud, numdup, counts, cornum, clustnummer, modelnames, noise)
#=======================> DONE!
@@ -828,7 +916,6 @@
scores = lapply(coffeeshop[,1], scorewindow, coffeeshop[,2], numdup, C, bkgd, counts, startlist = coffeeshop[,1])
}
scores = unlist(scores)
- #scores = scores[is.finite(scores)] #just a sanity check
if (windowe == "auto") {
lscores = log(scores)
@@ -867,12 +954,12 @@
# ==============
# Finding Peaks
# ==============
-#if (length(coffeeshop[,1]) > 0) { #any enriched windows left after filtering?, just a sanity check!
+if (length(coffeeshop[,1]) > 0) { #any enriched windows left after filtering?
coffeeshop = coffeeshop[,-3]
if (cornum > 1) {
- peaks = mclapply(coffeeshop[,1], findpeak, coffeeshop[,2], numdup, C, param, bkgd, resol, counts, noise, startlist = coffeeshop[,1], meanAdjust, mc.cores = cornum, mc.preschedule = TRUE)
+ peaks = mclapply(coffeeshop[,1], findpeak, coffeeshop[,2], numdup, C, param, bkgd, resol, counts, noise, startlist = coffeeshop[,1], meanAdjust, clustnummer, mc.cores = cornum, mc.preschedule = TRUE)
} else {
- peaks = lapply(coffeeshop[,1], findpeak, coffeeshop[,2], numdup, C, param, bkgd, resol, counts, noise, startlist = coffeeshop[,1], meanAdjust)
+ peaks = lapply(coffeeshop[,1], findpeak, coffeeshop[,2], numdup, C, param, bkgd, resol, counts, noise, startlist = coffeeshop[,1], meanAdjust, clustnummer)
}
if (!(is.null(peaks))) { #any peaks discovered?
writethis = processPeaks(peaks)
@@ -886,15 +973,16 @@
# Writing Peak Information
# =========================
} else { nothing = TRUE } #no peaks
-#} else { nothing = TRUE } #no enriched windows left after filtering, just a sanity check (at least one window should still be there)
+} else { nothing = TRUE } #no enriched windows left after filtering
} else { nothing = TRUE } #no enriched widnows discovered
+} else { nothing = TRUE; cutthisW = windowe } #no enriched bins discovered
if (isTRUE(nothing)) {
file.create(paste0(out, "/", chromName, ".peaks.bed"))
write(paste(chromName, minpeak, sep = " "), file = paste0(out, "/min.peaksize"), append=TRUE)
if (chrcount == 1) {
- message(paste0("No peaks found! - Window Fold Enrichment: ", round(cutthisW,2)))
+ message(paste0("No peaks found! - Window Fold Enrichment: ", cutthisW, " - Seed: ", samplingSeed))
} else {
message("No peaks found!")
}
@@ -904,7 +992,7 @@
if (chrcount == 1) {
- message(paste0("Done! - Window Fold Enrichment: ", round(cutthisW,2)))
+ message(paste0("Done! - Window Fold Enrichment: ", cutthisW, " - Seed: ", samplingSeed))
} else {
message("Done!")
}
diff -r d42f4d78c85e -r 243f75d0ed6e peakhelper.r
--- a/peakhelper.r Wed Dec 17 10:40:23 2014 -0500
+++ b/peakhelper.r Thu Feb 19 05:39:45 2015 -0500
@@ -1,8 +1,22 @@
-#########################################################################################
-#### Fixes output of peakfinder.r to produce proper narrowPeak format with proper scores
-#### R script
-########################################################################################
-
+########################################################################
+# JAMMv1.0.7rev1 is a peak finder for joint analysis of NGS replicates.
+# Copyright (C) 2014-2015 Mahmoud Ibrahim
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# Contact: mahmoud.ibrahim@mdc-berlin.de
+########################################################################
@@ -38,20 +52,6 @@
-
-# ==============================================
-# Required Libraries check (Source: http://r.789695.n4.nabble.com/test-if-a-package-is-installed-td1750671.htm)
-# ==============================================
-if ((is.element('sqldf', installed.packages()[,1]) == FALSE)) {
- stop("R package 'sqldf' is required. Please install it!")
-}
-suppressPackageStartupMessages(library("sqldf"))
-suppressPackageStartupMessages(library("tcltk"))
-#=======================> DONE!
-
-
-
-
# =================
# Custom Functions
# =================
diff -r d42f4d78c85e -r 243f75d0ed6e readme
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/readme Thu Feb 19 05:39:45 2015 -0500
@@ -0,0 +1,21 @@
+JAMM is a universal peak finder for NGS datasets that can integrate information from biological replicates.
+
+JAMM was developed in the Ohler lab at the Berlin Institute for Medical Systems Biology.
+
+
+Visit this page for documentation: http://code.google.com/p/jamm-peak-finder/
+
+
+To install JAMM:
+- Make sure you are on a Linux system
+- Make sure you have R and Perl installed
+- Install the following R packages: sqldf, signal, mclust
+
+
+To get information about options available to JAMM users, run the JAMM bash script without any parameters.
+
+
+If you encounter errors or have questions, please email mahmoud.ibrahim@mdc-berlin.de or post at the JAMM Google groups page: https://groups.google.com/forum/#!forum/jamm-peak-finder
+
+
+Thank you.
diff -r d42f4d78c85e -r 243f75d0ed6e readshifter.pl
--- a/readshifter.pl Wed Dec 17 10:40:23 2014 -0500
+++ b/readshifter.pl Thu Feb 19 05:39:45 2015 -0500
@@ -1,8 +1,22 @@
-#!/usr/bin/perl
-##############################################
-#### Shifts minus strand reads by read length
-#### Perl script
-#############################################
+########################################################################
+# JAMMv1.0.7rev1 is a peak finder for joint analysis of NGS replicates.
+# Copyright (C) 2014-2015 Mahmoud Ibrahim
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# Contact: mahmoud.ibrahim@mdc-berlin.de
+########################################################################
diff -r d42f4d78c85e -r 243f75d0ed6e signalmaker.r
--- a/signalmaker.r Wed Dec 17 10:40:23 2014 -0500
+++ b/signalmaker.r Thu Feb 19 05:39:45 2015 -0500
@@ -1,9 +1,22 @@
-###############################
-#### Signal Maker for NGS Data
-#### R script
-###############################
-
-
+########################################################################
+# JAMMv1.0.7rev1 is a peak finder for joint analysis of NGS replicates.
+# Copyright (C) 2014-2015 Mahmoud Ibrahim
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# Contact: mahmoud.ibrahim@mdc-berlin.de
+########################################################################
# =======================
@@ -76,10 +89,15 @@
return(reads)
}
-
+#Read in bedpe(st2)
+parseinpe = function(bedfile) {
+ l = read.table(bedfile, header = FALSE)
+ l = cbind((l[[1]] + 1), l[[2]])
+ return(l)
+}
#Produces normalized extended read counts (takes output of parsein(), return a vector of floats)
-countreads = function(bedfile, reads, frag, chromsize, filelist) {
+countreads = function(bedfile, reads, frag, chromsize, filelist, chrcount) {
o = which(filelist == bedfile)
@@ -112,6 +130,36 @@
+#Produces normalized extended read counts (takes output of parsein(), return a vector of floats / PAIRED END)
+countreadspe = function(bedfile, reads, chromsize, filelist, chrcount) {
+
+ o = which(filelist == bedfile)
+
+ counts = vector(mode = "numeric", length = chromsize)
+
+ for (j in 1:length(reads[[o]][,1])) {
+ counts[(reads[[o]][j,1]):(reads[[o]][j,2])] = counts[(reads[[o]][j,1]):(reads[[o]][j,2])] + 1
+ }
+
+ mCount = mean(counts)
+
+ if (chrcount == 1) {
+ counts = counts/mCount
+ write(paste(mCount), file = paste0(out, "/norma.", o, ".info"))
+ } else {
+ meanCounts = mean(as.numeric(read.table(paste0(out, "/norma.", o, ".info"))[[1]]))
+ if ((mCount > (5*meanCounts)) || (mCount < (0.2*meanCounts))) {
+ mCount = meanCounts
+ } else {
+ write(paste(mCount), file = paste0(out, "/norma.", o, ".info"), append = TRUE)
+ }
+ counts = counts/mCount
+ }
+
+ return(counts)
+}
+
+
#generate signal
getsig = function(winStart, winEnds, startList, stepSize, chromName) {
@@ -187,10 +235,8 @@
if (grepl('^-frag=',each.arg)) {
arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
if (! is.na(arg.split[2]) ) {
- frag <- arg.split[2]
- } else {
- stop('No fragment length given')
- }
+ frag <- arg.split[2]
+ }
}
#bakcground files directory
if (grepl('^-bkgd=',each.arg)) {
@@ -252,6 +298,13 @@
chrcount <- as.numeric(arg.split[2])
}
}
+ #alignment type
+ if (grepl('^-t=',each.arg)) {
+ arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]]
+ if (! is.na(arg.split[2]) ) {
+ type <