coverage strand count not accurate compared to samtools for paired end alignment

[mheskett]

Using coverage -c -s where -a is a bed file with strand info, and -b is a strand specific BAM alignment (RNA-seq in my case), I get very different values compared to using the standard samtools commands to extract stranded. This is what I do in samtools for PE library

## get plus strand
samtools view -b -f 131 -F 16 $input > $out_dir$filename.fwd1.bam
samtools index $out_dir$filename.fwd1.bam

# now get first in pair mapping to reverse strand
samtools view -b -f 83 $input > $out_dir$filename.fwd2.bam
samtools index $out_dir$filename.fwd2.bam

# now combine, this should contail all plus strand genes
samtools merge -f $out_dir$filename.plus.bam $out_dir$filename.fwd1.bam $out_dir$filename.fwd2.bam
samtools index $out_dir$filename.plus.bam

and the converse for minus strand