Wrapper for htseq-count

htseq_count_em.sh

#!/bin/bash
samdir="./"
annotation="my.gtf"
htseqdir="./htseq-count"
htseqcount=`which htseq-count`
samtools=`which samtools`
feature="gene"
idattr="gene_name"
ss="reverse"

if [ -d "$htseqdir" ];
 then   
        echo "$htseqdir available"
else    
        mkdir $htseqdir
fi


for BAM in $(ls $samdir/*.bam)
do
bn=$(basename $BAM .uniq.bam)
echo "processing $bn"
set -x
$samtools sort -n -o $BAM $bn | $samtools view -h - | $htseqcount -t $feature  -i $idattr -s $ss - $annotation > $htseqdir/$bn.htseq.count 2> $htseqdir/$bn.htseq.count.err
set +x
done