I am analyzing an rna-seq data set, which consists of matched tumor/normal samples from two patients. Alignment of the trimmed fastq files were done using bowtie2 with default parameters and iGenomes bowtie2 hg19 index. The alignement BAM files were sorted by name and coordiantes using samtools. Now I am trying to get the counts for further downstream analysis. I have tried several approaches, including htseq-count, countOverlaps, and summarizeOverlaps. Each of them gives slightly different results - i.e., the counts for the same gene and sample are different. My question is that is there a way to extract a subset of alignments from the BAM files and create a "subsetted" BAM file that can be analyzed in depth with different methods so as to get an idea why there are differences. All the methods for getting the counts are really "black boxes" for me, and I am not sure what is going on under the covers to produce different counts. I thought it would be a good idea to just try out different methods with a set of alignments that have mapped to a particular gene (or even exon) and try different methods with that subset so it is more tractable and faster to try out and compare various methods.