This script follows from 01.SE/PE_chip_atac_process.sh to shift reads per the ATAC-seq protocol, producing new read-shifted bam files, and also removes contaminating mitochondrial reads.

B1_atac_shift.sh

Arguments

PREFIX

sample directory where aligned bam files are located in a" samtools directory with *-bt2-st.bam suffix"

THREADS

allocated threads for the job; note that memory amt" proportional to number of threads / 20 * 128G (Ruddle)" which ends up being about 6GB per thread"

Details

This script requires the following:

  • bedtools v2.22.1

  • samtools 1.3.1

Examples


#' Preset parameters
MEMPERTHREAD=6                 # GB; 128 GB of RAM / 20 cores per node

#' 0. Read in Parameters
SDIR=$(readlink -f $1)            # Sample directory
PREFIX=$(basename $SDIR)          # prefix for files 
THREADS=$2


#' 0. Echo parameters back
echo ""
echo "Beginning to process your paired/single-end ATAC-seq data.."
echo ".. Memory per thread is: $MEMPERTHREAD GB"
echo ".. Using for Chromsome sizes: $CHRSIZES"
echo ""
echo "Working on: $SDIR"
echo ".. Using prefix: $PREFIX"
echo ".. Using $THREADS threads"
echo ""


#' 1. Remove Mitochondrial reads and shift reads
echo "1. Removing mitochondrial reads, and shifting by +4/-5 for +/- strands.."
#' Note: uses bitwise flags for paired end:
#'   0  and  16 are +/- for (unpaired) single-end
#'   99 and 147 are +, with 99 being mate1 (left)
#'   83 and 163 are -, with 83 being mate1 (left)
mkdir -p $SDIR/shift
cd $SDIR/shift

samtools view -h $SDIR/samtools/${PREFIX}_bt2-st.bam | \
    awk '{ FS = "\t"; OFS = "\t" } { 
  if (substr($1,1,1) == "@") { print $0 } else
  if ($3 == "MT" || $3 == "M") {        } else
  if ($2 == 0 || $2 == 99 || $2 == 147)  { $4 = $4 + 4; print $0 } else 
  if ($2 == 16 || $2 == 83 || $2 == 163) { $4 = $4 - 5; print $0 } 
  else { print $0 }
}' | samtools view -b - > tmp_shift.bam

#' 2. Sort and index resulting shifted file
echo "2. Sorting bam file and indexing.."
samtools sort -@ $THREADS -m $(($MEMPERTHREAD * $THREADS))G \
    -T tmp -O bam tmp_shift.bam > ${PREFIX}_shift.bam
samtools index ${PREFIX}_shift.bam

#' 3. Clean up
echo "3. Cleaning up.."
rm tmp_shift.bam

#' 3. Success?
echo ""
echo "Shifting completed."