Commit ab808a5a authored by nchernia's avatar nchernia
Browse files

added LSF and CPU beta versions. modified the mega script in AWS and LSF....

added LSF and CPU beta versions. modified the mega script in AWS and LSF. Juicer 1.6 will reconcile these different issues
parent ada6bde8
Loading
Loading
Loading
Loading
+70 −0
Original line number Diff line number Diff line
#!/bin/bash
##########
#The MIT License (MIT)
#
# Copyright (c) 2015 Aiden Lab
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
#  THE SOFTWARE.
##########
# Juicer postprocessing script.
# This will find contact domains with Arrowhead.

## Read arguments
usageHelp="Usage: ${0} [-h] -j <juicebox_file_path> -i <hic_file_path>"

printHelpAndExit() {
    echo "$usageHelp"
    exit $1
}

#set defaults
genomeID="hg19"
hic_file_path="$(pwd)/aligned/inter_30.hic"
juiceboxpath="/opt/juicer/scripts/juicebox"

while getopts "h:j:i:" opt; do
    case $opt in
	h) printHelpAndExit 0;;
	j) juiceboxpath=$OPTARG ;;
	i) hic_file_path=$OPTARG ;;
	[?]) printHelpAndExit 1;;
    esac
done

## Check that juicebox exists 
if [ ! -e "${juiceboxpath}" ]; then
  echo "***! Can't find juicebox in ${juiceboxpath}";
  exit 100;
fi

## Check that hic file exists    
if [ ! -e "${hic_file_path}" ]; then
  echo "***! Can't find inter.hic in ${hic_file_path}";
  exit 100;
fi

echo -e "${juiceboxpath} is post-processing Hi-C for ${genomeID}\nData read from ${hic_file_path}.\nMotifs read from ${bed_file_dir}\n"
echo -e "ARROWHEAD:\n"
${juiceboxpath} arrowhead ${hic_file_path} ${hic_file_path%.*}"_contact_domains.txt"
if [ $? -ne 0 ]; then
    echo "***! Problem while running Arrowhead";
    exit 100
else
    echo -e "\n(-: Arrowhead Postprocessing successfully completed (-:"
fi
+90 −0
Original line number Diff line number Diff line
#!/bin/bash
##########
#The MIT License (MIT)
#
# Copyright (c) 2015 Aiden Lab
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
#  THE SOFTWARE.
##########
# Juicer postprocessing script.
# This will run the major post-processing on the HiC file, including finding
# loops with HiCCUPS and inding motifs of these loops with MotifFinder.

## Read arguments
usageHelp="Usage: ${0} [-h] -j <juicebox_file_path> -i <hic_file_path> -m <bed_file_dir> -g <genome ID>"

printHelpAndExit() {
    echo "$usageHelp"
    exit $1
}

#set defaults
genomeID="hg19"
hic_file_path="$(pwd)/aligned/inter_30.hic"
juiceboxpath="/opt/juicer/scripts/juicebox"
bed_file_dir="/opt/juicer/references/motif"

while getopts "h:g:j:i:m:" opt; do
    case $opt in
	h) printHelpAndExit 0;;
	j) juiceboxpath=$OPTARG ;;
	i) hic_file_path=$OPTARG ;;
	m) bed_file_dir=$OPTARG ;; 
	g) genomeID=$OPTARG ;;
	[?]) printHelpAndExit 1;;
    esac
done

## Check that juicebox exists 
if [ ! -e "${juiceboxpath}" ]; then
  echo "***! Can't find juicebox in ${juiceboxpath}";
  exit 100;
fi

## Check that hic file exists    
if [ ! -e "${hic_file_path}" ]; then
  echo "***! Can't find inter.hic in ${hic_file_path}";
  exit 100;
fi

## Check that bed folder exists    
if [ ! -e "${bed_file_dir}" ]; then
  echo "***! Can't find folder ${bed_file_dir}";
  exit 100;
fi

echo -e "\nHiCCUPS:\n"
${juiceboxpath} hiccups ${hic_file_path} ${hic_file_path%.*}"_loops.txt"
if [ $? -ne 0 ]; then
    echo "***! Problem while running HiCCUPS";
    exit 100
fi

if [ -f ${hic_file_path%.*}"_loops.txt" ]
then
    echo -e "\nAPA:\n"
    ${juiceboxpath} apa ${hic_file_path} ${hic_file_path%.*}"_loops.txt" "apa_results"
    echo -e "\nMOTIF FINDER:\n"
    ${juiceboxpath} motifs ${genomeID} ${bed_file_dir} ${hic_file_path%.*}"_loops.txt"
    echo -e "\n(-: Feature annotation successfully completed (-:"
else
    # if loop lists do not exist but Juicebox didn't return an error, likely 
    # too sparse
    echo -e "\n(-: Postprocessing successfully completed, maps too sparse to annotate (-:"
fi
+112 −60
Original line number Diff line number Diff line
@@ -43,10 +43,10 @@ juicer_version="1.5"
#load_gpu=""
# Juicer directory, contains scripts/ and restriction_sites/
juiceDir="/opt/juicer"
# default queue, can also be set in options
# default queue and time
queue="Stat"
# default long queue, can also be set in options
long_queue="Merge"
queue_time="12:00"

# unique name for jobs in this run
groupname="a$(date +%s)"

@@ -55,15 +55,16 @@ groupname="a$(date +%s)"
# top level directory, can also be set in options
topDir=$(pwd)
# restriction enzyme, can also be set in options
site="DpnII"
site="MboI"
# genome ID, default to human, can also be set in options
genomeID="hg19"

## Read arguments
usageHelp="Usage: ${0##*/} -g genomeID [-d topDir] [-s site] [-hx]"
usageHelp="Usage: ${0##*/} -g genomeID [-d topDir] [-s site] [-r resolutions] [-hx]"
genomeHelp="   genomeID must be defined in the script, e.g. \"hg19\" or \"mm10\" (default \"$genomeID\")"
dirHelp="   [topDir] is the top level directory (default \"$topDir\") and must contain links to all merged_nodups files underneath it"
siteHelp="   [site] must be defined in the script, e.g.  \"HindIII\" or \"MboI\" (default \"$site\"); alternatively, this can be the restriction site file"
resolutionsHelp="   [resolutions] is a comma-delimited list of resolutions, such as 10000,5000,1000,5f (default is 2.5M,1M,500K,250K,100K,50K,25K,10K,5K in base pair and 500f,250f,100f,50f,25f,10f,5f,2f,1f)"
excludeHelp="   -x: exclude fragment-delimited maps from Hi-C mega map (will run much faster)"
helpHelp="   -h: print this help and exit"

@@ -72,18 +73,20 @@ printHelpAndExit() {
    echo "$genomeHelp"
    echo "$dirHelp"
    echo "$siteHelp"
    echo "$resolutionsHelp"
    echo "$excludeHelp"
    echo "$helpHelp"
    exit "$1"
}

while getopts "d:g:hxs:" opt; do
while getopts "d:g:r:hxs:" opt; do
    case $opt in
	g) genomeID=$OPTARG ;;
	h) printHelpAndExit 0;;
	d) topDir=$OPTARG ;;
	s) site=$OPTARG ;;
	x) exclude=1 ;;
	r) resolutions=$OPTARG ;;
	[?]) printHelpAndExit 1;;
    esac
done
@@ -117,9 +120,10 @@ fi
megadir=${topDir}"/mega"
outputdir=${megadir}"/aligned"
tmpdir=${megadir}"/HIC_tmp"
export TMPDIR=${megadir}"/HIC_tmp"
outfile=${megadir}/lsf.out
touchfile1=${megadir}/touch1
# set global tmpdir so no problems with /var/tmp
export TMPDIR=${tmpdir}
#output messages
logdir=${megadir}/debug

## Check for existing merge_nodups files:

@@ -130,9 +134,27 @@ then
	exit 100
fi

merged_names=$(find -L ${topDir} | grep merged_nodups.txt | tr '\n' ' ')
merged_names1=$(find -L ${topDir} | grep merged_nodups.txt)
merged_names=$(echo $merged_names1 | tr '\n' ' ')
inter_names=$(find -L ${topDir} | grep inter.txt | tr '\n' ' ')

if [[ $merged_names == *".txt.gz"* ]]
then
    gzipped=1
    echo "***! Mega map of gzipped files not yet supported, please unzip before running."
    exit 100
    # we need to unzip here
    for i in $merged_names1
    do
	if [[ $i != *".txt.gz"* ]]
	then
	    echo "***! Mixture of gzipped and unzipped merged_nodups files"
	    echo "Ensure that the merged_nodups are all either unzipped or gzipped then rerun"
	    echo "Files: $merged_names"
	    exit 100
	fi
    done
fi
## Create output directory, exit if already exists
if [[ -d "${outputdir}" ]] 
then
@@ -148,17 +170,25 @@ if [ ! -d "$tmpdir" ]; then
    chmod 777 $tmpdir
fi

#Prepare an empty log file:
touch ${outfile}
chmod 777 ${outfile}
## Create log directory
if [ ! -d "$logdir" ]; then
    mkdir $logdir
    chmod 777 $logdir
fi

if [ -n "$resolutions" ]; then
    resolutions="-r $resolutions"
fi

## Arguments have been checked and directories created. Now begins
## the real work of the pipeline

#source $usePath
#$load_cluster

jid1=`bsub -o ${megadir}/lsf.out -q "${queue}" -J "${groupname}_topstats" <<-TOPSTATS
bsub -o ${logdir}/header.out -q "${queue}" -J ${groupname}cmd <<-EOF
  date
  echo "Juicer version:$juicer_version"
  echo "$0 $@"
EOF
touchfile1=${megadir}/touch1
jid1=`bsub -o ${logdir}/topstats.out -e ${logdir}/topstats.err -q "${queue}" -J "${groupname}_topstats" <<-TOPSTATS

export LC_ALL=C
if ! awk -f ${juiceDir}/scripts/makemega_addstats.awk ${inter_names} > ${outputdir}/inter.txt
@@ -173,15 +203,10 @@ TOPSTATS`

touchfile2=${megadir}/touch2
# Merge all merged_nodups.txt files found under current dir
jid2=`bsub -o ${megadir}/lsf.out -q "${long_queue}" -J ${groupname}_merge -R "rusage[mem=16000]" -w "done(${groupname}_topstats)" <<- MRGSRT
if [ ! -f "${touchfile1}" ]
then
   echo "***! Top stats job failed, type bjobs -l $jid1 to see what happened."
   exit 100;
fi
jid2=`bsub -o ${logdir}/merge.out -e ${logdir}/merge.err -q "${queue}" -J ${groupname}_merge -R "rusage[mem=16000]" <<- MRGSRT
if ! sort -T ${tmpdir} -m -k2,2d -k6,6d ${merged_names} > ${outputdir}/merged_nodups.txt
then 
echo "***! Some problems occurred somewhere in creating sorted merged_nodups files."
echo "***! Some problems occurred somewhere in merging sorted merged_nodups files."
    exit 100
else
echo "(-: Finished sorting all merged_nodups files into a single merge."
@@ -193,7 +218,12 @@ MRGSRT`
touchfile3=${megadir}/touch3  

# Create statistics files for MQ > 0
jid3=`bsub -o ${megadir}/lsf.out -q "${queue}" -J "${groupname}_inter0" -w "done(${groupname}_merge)" <<- INTER0
jid3=`bsub -o ${logdir}/inter0.out -e ${logdir}/inter0.err -q "${queue}" -J "${groupname}_inter0" -w "done(${groupname}_merge) && done(${groupname}_topstats)" <<- INTER0
if [ ! -f "${touchfile1}" ]
then
   echo "***! Top stats job failed, type bjobs -l $jid1 to see what happened."
   exit 100;
fi
if [ ! -f "${touchfile2}" ]
then
   echo "***! Sort job failed."
@@ -207,7 +237,13 @@ INTER0`
touchfile4=${megadir}/touch4

# Create statistics files for MQ > 30
jid4=`bsub -o ${megadir}/lsf.out -q "${queue}" -J "${groupname}_inter30" -w "done(${groupname}_merge)" <<- INTER30
jid4=`bsub -o ${logdir}/inter30.out -e ${logdir}/inter30.err -q "${queue}" -J "${groupname}_inter30" -w "done(${groupname}_merge) && done(${groupname}_topstats)" <<- INTER30
if [ ! -f "${touchfile1}" ]
then
   echo "***! Top stats job failed, type bjobs -l $jid1 to see what happened."
   exit 100;
fi

if [ ! -f "${touchfile2}" ]
then
   echo "***! Sort job failed."
@@ -220,24 +256,24 @@ INTER30`
touchfile5=${megadir}/touch5

# Create HIC maps file for MQ > 0
jid5=`bsub -o ${megadir}/lsf.out -q "${queue}" -J "${groupname}_hic0" -w "done(${groupname}_inter30)" -R "rusage[mem=16000]" <<- HIC0
#source $usePath
#$load_java
jid5=`bsub -o ${logdir}/hic0.out -e ${logdir}/hic0.err -q "${queue}" -J "${groupname}_hic0" -w "done(${groupname}_inter30)" -R "rusage[mem=16000]" <<- HIC0
if [ ! -f "${touchfile3}" ]
then
   echo "***! Statistics q=1 job failed."
   echo "***! Statistics q=1 job failed, type bjobs -l $jid3 to see what happened"
   exit 100;
fi
exitcode=-999
myexitcode=1
if [ -z "$exclude" ]
then
	${juiceDir}/scripts/juicebox pre -f ${site_file} -s ${outputdir}/inter.txt -g ${outputdir}/inter_hists.m -q 1 ${outputdir}/merged_nodups.txt ${outputdir}/inter.hic ${genomeID}
    exitcode=\$?
  echo "Launching ${juiceDir}/scripts/juicebox pre $resolutions -f ${site_file} -s ${outputdir}/inter.txt -g ${outputdir}/inter_hists.m -q 1 ${outputdir}/merged_nodups.txt ${outputdir}/inter.hic ${genomeID}"
  ${juiceDir}/scripts/juicebox pre $resolutions -f ${site_file} -s ${outputdir}/inter.txt -g ${outputdir}/inter_hists.m -q 1 ${outputdir}/merged_nodups.txt ${outputdir}/inter.hic ${genomeID}
  myexitcode=\$?
else
	${juiceDir}/scripts/juicebox pre -s ${outputdir}/inter.txt -g ${outputdir}/inter_hists.m -q 1 ${outputdir}/merged_nodups.txt ${outputdir}/inter.hic ${genomeID}
    exitcode=\$?
  echo "Launching ${juiceDir}/scripts/juicebox pre $resolutions -s ${outputdir}/inter.txt -g ${outputdir}/inter_hists.m -q 1 ${outputdir}/merged_nodups.txt ${outputdir}/inter.hic ${genomeID}"
  ${juiceDir}/scripts/juicebox pre $resolutions -s ${outputdir}/inter.txt -g ${outputdir}/inter_hists.m -q 1 ${outputdir}/merged_nodups.txt ${outputdir}/inter.hic ${genomeID}
  myexitcode=\$?
fi
if [ "\${exitcode}" -eq 0 ]
if [ "\${myexitcode}" -eq 0 ]
then
    touch $touchfile5
fi
@@ -245,24 +281,24 @@ HIC0`

touchfile6=${megadir}/touch6
# Create HIC maps file for MQ > 30
jid6=`bsub -o ${megadir}/lsf.out -q "${queue}" -J "${groupname}_hic30" -w "done(${groupname}_inter30)" -R "rusage[mem=16000]" <<- HIC30
#source $usePath
#$load_java	
jid6=`bsub -o ${logdir}/hic30.out -e ${logdir}/hic30.err -q "${queue}" -J "${groupname}_hic30" -w "done(${groupname}_inter30)" -R "rusage[mem=16000]" <<- HIC30
if [ ! -f "${touchfile4}" ]
then
   echo "***! Statistics q=30 job failed."
   echo "***! Statistics q=30 job failed, type bjobs -l $jid4 to see what happened."
   exit 100;
fi
exitcode=-999
myexitcode=1
if [ -z "${exclude}" ]
then
	${juiceDir}/scripts/juicebox pre -f ${site_file} -s ${outputdir}/inter_30.txt -g ${outputdir}/inter_30_hists.m -q 30 ${outputdir}/merged_nodups.txt ${outputdir}/inter_30.hic ${genomeID}
   exitcode=\$?
    echo "Launching ${juiceDir}/scripts/juicebox pre $resolutions -f ${site_file} -s ${outputdir}/inter_30.txt -g ${outputdir}/inter_30_hists.m -q 30 ${outputdir}/merged_nodups.txt ${outputdir}/inter_30.hic ${genomeID}"
    ${juiceDir}/scripts/juicebox pre $resolutions -f ${site_file} -s ${outputdir}/inter_30.txt -g ${outputdir}/inter_30_hists.m -q 30 ${outputdir}/merged_nodups.txt ${outputdir}/inter_30.hic ${genomeID}
   myexitcode=\$?
else
	${juiceDir}/scripts/juicebox pre -s ${outputdir}/inter_30.txt -g ${outputdir}/inter_30_hists.m -q 30 ${outputdir}/merged_nodups.txt ${outputdir}/inter_30.hic ${genomeID}
   exitcode=\$?
   echo "Launching ${juiceDir}/scripts/juicebox pre $resolutions -s ${outputdir}/inter_30.txt -g ${outputdir}/inter_30_hists.m -q 30 ${outputdir}/merged_nodups.txt ${outputdir}/inter_30.hic ${genomeID}"
   ${juiceDir}/scripts/juicebox pre $resolutions -s ${outputdir}/inter_30.txt -g ${outputdir}/inter_30_hists.m -q 30 ${outputdir}/merged_nodups.txt ${outputdir}/inter_30.hic ${genomeID}
   myexitcode=\$?
fi
if [ "\${exitcode}" -eq 0 ]
if [ "\${myexitcode}" -eq 0 ]
then
touch $touchfile6
fi
@@ -270,39 +306,55 @@ HIC30`

touchfile7=${megadir}/touch7
# Create loop and domain lists file for MQ > 30
jid7=`bsub  -o ${megadir}/lsf.out -q "${queue}" -J "${groupname}_postproc" -w "done(${groupname}_hic30)" <<- POSTPROC
#source $usePath
#$load_java	
jid7=`bsub  -o ${logdir}/hiccups.out -e ${logdir}/hiccups.err -q "${queue}" -J "${groupname}_hiccups" -w "done(${groupname}_hic30)" <<- HICCUPS
export _JAVA_OPTIONS=-Xmx16384m;
export LC_ALL=C
if [ ! -f "${touchfile6}" ]
then
   echo "***! Failed to create domain and loop lists"
   echo "***! Failed to make inter_30.hic, type bjobs -l $jid6 to see what happened."
   exit 100;
fi
${juiceDir}/scripts/juicer_postprocessing.sh -j ${juiceDir}/scripts/juicebox -i ${outputdir}/inter_30.hic -m ${juiceDir}/references/motif -g ${genomeID}
${juiceDir}/scripts/juicer_hiccups.sh -j ${juiceDir}/scripts/juicebox -i ${outputdir}/inter_30.hic -m ${juiceDir}/references/motif -g ${genomeID}
touch $touchfile7
POSTPROC`
HICCUPS`

jid8=`bsub -o ${megadir}/lsf.out -j y -q "${queue}" -J "${groupname}_done" -w "done(${groupname}_postproc)" <<- FINAL
touchfile8=${megadir}/touch8
jid8=`bsub  -o ${logdir}/arrowhead.out -e ${logdir}/arrowhead.err -q "${queue}" -J "${groupname}_hiccups" -w "done(${groupname}_hic30)" <<- ARROWHEAD
export _JAVA_OPTIONS=-Xmx16384m;
export LC_ALL=C
if [ ! -f "${touchfile6}" ]
then
   echo "***! Failed to make inter_30.hic, type bjobs -l $jid6 to see what happened."
   exit 100;
fi
${juiceDir}/scripts/juicer_arrowhead.sh -j ${juiceDir}/scripts/juicebox -i ${outputdir}/inter_30.hic
touch $touchfile8
ARROWHEAD`

bsub -o ${megadir}/done.out -e ${megadir}/done.err -q "${queue}" -J "${groupname}_done" -w "done(${groupname}_hiccups) && done(${groupname}_arrowhead)" <<- FINAL
if [ ! -f "${touchfile5}" ]
then
   echo "***! Failed to make inter.hic."   
   echo "***! Failed to make inter.hic, type bjobs -l $jid5 to see what happened."   
   exit 100;
fi
if [ ! -f "${touchfile6}" ]
then
   echo "***! Failed to make inter_30.hic."   
   echo "***! Failed to make inter_30.hic, type bjobs -l $jid6 to see what happened."   
   exit 100;
fi
if [ ! -f "${touchfile7}" ]
then
   echo "***! Failed in postprocessing"   
   echo "***! Failed to create loop lists, type bjobs -l $jid7 to see what happened."   
   exit 100;
fi

rm $touchfile1 $touchfile2 $touchfile3 $touchfile4 $touchfile5 $touchfile6 $touchfile7
if [ ! -f "${touchfile8}" ]
then
   echo "***! Failed to create contact domain lists, type bjobs -l $jid8 to see what happened."   
   exit 100;
fi
rm -r ${tmpdir}
rm $touchfile1 $touchfile2 $touchfile3 $touchfile4 $touchfile5 $touchfile6 $touchfile7 $touchfile8
echo "(-: Successfully completed making mega map. Done. :-)"
FINAL`
FINAL

echo "(-: Finished adding all jobs... please wait while processing."

CPU/README.md

0 → 100644
+7 −0
Original line number Diff line number Diff line
# Juicer single CPU version BETA

**BETA VERSION** 

Report bugs using the Issues tab. This script has not been thoroughly tested.

Make sure you've installed BWA and Java-1.7 before proceeding.
+4.55 KiB

File added.

No diff preview for this file type.

Loading