I needed to check the results coming from featureCounts from subread. In fact, samtools outputted less reads then featurecounts did.
I run samtools as follows:
samtools view -c -q 1 -F 4 -L matrix_gene.bed file.bam
2732259
And featureCounts:
featureCounts --minReadOverlap 50 -s 1 -Q 1 -T 12 -a matrix_gene.saf -F SAF -O -o $OUT/total_counts.txt file.bam
featureCounts --minReadOverlap 50 -s 2 -Q 1 -T 12 -a matrix_gene.saf -F SAF -O -o $OUT/total_counts.txt file.bam
(minReadOverlap
= 50 because the length of my reads is 50, -s 1
for stranded and -s 2
for reversely stranded, -Q
mapping quality and -O
for multiple overlap (as samtools outputs alignment multiple times if they overlap multiple regions)
samtools results: 2732259
featureCounts -s 1: 1661976
featureCounts -s 2: 1662917
So, the number of stranded and reverse stranded reads: 3324893
Why do the results differ if I specified the same parameters in both methods?
PS: I did not find what overlap between a read and a region samtools allows, but if the overlap less than the read length, samtools should have more reads than featureCounts.