thank you for your reply, for featureCounts, I use -s 2 to specify the strandness, which is equivalent to the -s reverse in htseq-count
# Strandness
-s <int or string> Perform strand-specific read counting. A single integer
value (applied to all input files) or a string of comma-
separated values (applied to each corresponding input
file) should be provided. Possible values include:
0 (unstranded), 1 (stranded) and 2 (reversely stranded).
Default value is 0 (ie. unstranded read counting carried
out for all input files).
If you are not sure about the best method to count reads or if you want to compare methods, I would pick some gene IDs with small counts differences between methods on a number of countable reads (like ENSG00000284431, one can count 10 reads by eye) and check your bam file on IGV for this specific gene ID to get an idea of how the counting could have been done. Maybe you will see that 1 read is overlapping 2 genes that could be discarded in one method but not the other.
Thank you for your help. I'll definitely try the method you suggested. I was just a bit puzzled because I had initially thought that htseq-count defaults to counting reads rather than read pairs (fragments). However, I found in the Q&A section of htseq-count that it states that for paired-end data, it actually defaults to counting fragments.
https://htseq.readthedocs.io/en/master/htseqcount.html
I choose one of the gene where featureCounts counts as 2 but htseq-count counts as 0. But... I don't know how to read the number of fragments from the screen .... (oh no.....
These tools will account for fragments at the moment you define paired-end reads as input, as the 2 reads of the same fragments are sequencing the same molecule.
Which gene are we talking about ? You have to check where the fragments start and end to understand why one method count them in or out