Count all reads, nucleosomal-free reads and mononucleosomal reads in a bam file.

samtools.count.sh

#!/bin/bash
#$ -N samtools.count
#$ -cwd
#$ -q long-sl7,mem_512_12h,mem_512
#$ -l virtual_free=50G,h_rt=12:00:00,disk=50G
#$ -M anamaria.elek@crg.eu
#$ -o samtools.count.out
#$ -pe smp 16
#$ -j y

SAMTOOLS_PATH="/users/asebe/aelek/bin/samtools-1.12"

echo $(date)

samtools_count () {
  local BAM=$1
  local COUNTS=${BAM%%bam}counts.txt
  local OUT=${BAM%%bam}counts.tmp
  local ROWS=$( dirname $BAM )/tmp.txt
  printf "FILE\nREADS\nNFR\nMONO\nOTHER" > $ROWS
  echo $(date): Counting reads in $BAM
  READS=$( ${SAMTOOLS_PATH}/samtools view -c $BAM )
  NFR=$( ${SAMTOOLS_PATH}/samtools view $BAM | awk -v LEN=147 '{if ( $9 <= LEN && $9 >= -(LEN) && $9 != 0) print $0}' | wc -l )
  MONO=$( ${SAMTOOLS_PATH}/samtools view $BAM | awk -v LENMIN=147 -v LENMAX=240 '{if ( $9 > LENMIN && $9 < LENMAX || $9 < -(LENMIN) && $9 > -(LENMAX) ) print $0}' | wc -l )
  OTHER=$( ${SAMTOOLS_PATH}/samtools view $BAM | awk -v LEN=240 '{if ( $9 >= LEN || $9 <= -(LEN) || $9 == 0) print $0}' | wc -l )
  printf "$BAM\n$READS\n$NFR\n$MONO\n$OTHER" > $OUT
  paste $ROWS $OUT > $COUNTS
  rm $OUT # $ROWS
  echo $(date): Done counting reads in $BAM
}



BAMS=$( ls /users/asebe/aelek/proj/scATAC_nvec*/scATAC_pro/output*/mapping_result/*positionsort.bam )
for BAM in $BAMS; do samtools_count "$BAM" & done
wait
echo "All done."