Calls peaks using MACS2 software, following from combining samples into condition/replicate aware directory structures (see D0.combine_by_condition.sh. The provided tab-delimited table used in that step can be reused for this script.

E0_call_peaks_macs2.sh

Arguments

COND

condition directory from a repdat folder where sample and inputs" are organized based on condition."

TABLE

tab-delimited table specifying the CONDITION_ID, SAMPLE_ID," and INPUT_ID as columns, and a fourth TYPE column for punctate" vs. broad peak-calling, with any additional columns being" ignored. The first line (header) will be ignored."

Details

This script will output a new folder, macs2, within the condition directory.

This script requires the following:

  • tab separated database specifying condition, sample, and input (NA if none)

  • MACS2 v2.1.0

See also

https://github.com/taoliu/MACS

Examples


## Preset parameters

## 0. Read in Parameters
COND=$(readlink -f $1)
TAB=$(readlink -f $2)

## 0. Echo parameters back
echo ""
echo "Beginning to process your ChIP-seq/ATAC-seq data.."
echo ""
echo "Working on samples in: $COND"
echo ""


#' 1. Routine for all samples within condition

#' Per sample within the sample directory
for I in $(ls -d $COND/samples/*); do 
    PREFIX=$(basename $I)
    echo ""
    echo ".. Calling peaks for: $PREFIX"

    #' Reading database table for parameters
    SAMPLEID=$(grep $PREFIX $TAB | awk '{print $2}')
    INPUTID=$(grep $PREFIX $TAB | awk '{print $3}')
    PEAKTYPE=$(grep $PREFIX $TAB | awk '{print $4}')
    SEQTYPE=$(grep $PREFIX $TAB | awk '{print $5}')

    #' Getting paired sample/input dirs
    SDIR=$(readlink -f $I)
    mkdir -p $SDIR/macs2
    cd $SDIR

    #' Interpreting database params
    if [ $PEAKTYPE == "broad" ]; then
	echo ".. Using broad setting.."
	CALLER="--broad"
    else
	echo ".. Using default punctate setting.."
	CALLER=""
    fi
  
    if [ $SEQTYPE == "paired" ]; then
	echo ".. Using paired-end setting.."
	FORMAT="BAMPE"
    else
	echo ".. Using single-end setting.."
	FORMAT="BAM"
    fi

    #' Interpreting if input or not
    if [ $INPUTID == "NA" ]; then
	echo ".. Calling peaks with no input.."
	echo ""
	macs2 callpeak \
	    -t $SDIR/bam_clean/${PREFIX}.bam \
	    -f $FORMAT -g mm -n ${PREFIX} \
	    --outdir $SDIR/macs2 $CALLER
    else
	IDIR=$(readlink -f $COND/inputs/$INPUTID)
	IPREFIX=$(basename $IDIR)
	echo ".. Calling peaks with input: $(basename $IDIR)"
	echo ""
	macs2 callpeak \
	    -t $SDIR/bam_clean/${PREFIX}.bam \
	    -c $IDIR/bam_clean/${IPREFIX}.bam \
	    -f $FORMAT -g mm -n ${PREFIX} \
	    --outdir $SDIR/macs2 $CALLER
    fi
   
    echo ""
done


echo ""
echo "Peak-calling completed."
echo ""