Question about featureCounts strand options

Question

[lorenzo]

Hello, I am new to bioinformatics and I am using featureCounts to quantify RNA-Seq reads but I don't grasp the behaviour and meaning of the strand specific options. I have bam files obtain from pair-end reads.

In particular, according to the official documentation:

• using the "-p" option: fragments will be counted instead of reads
• using the "-s 0" option: non strand specific counting is perfomed
• using the "-s 1" option: strand specific counting is performed
• using the "-s 2" option: strand specific counting is performed, but the strands are reversed.

My questions are:

• What is the actual behaviour of featureCounts with these options?
• With a bam file obtained from pair-end reads, should I always use the "-p" option? If not, why should I use the other ones?
• With my data, I observed that using the "-p" option I get about 70% of assigned alignments and I get about the same amount with the "-s 0" option, while I get 0% with "-s 1 / 2". Is this reasonable?

I checked: https://chipster.csc.fi/manual/library-type-summary.html and https://littlebitofdata.com/en/2017/08/strandness_in_rnaseq/ but it didn't help.

Thanks!

Answer 1

[h.mon]

With paired-end reads, you should always use -p.

It is really odd to get 70% assignment with -s 0 and 0% with -s 1 and -s 2. There is something funny either with your mapping or with your annotation. Is the annotation extremely compact, with all feature (genes) overlapping each other?

[lorenzo]

Using "-s 1" or "-s 2" the program warns me that "Paired-end reads were found and excluded", as you can see in the output snippet below. Can this be the cause?

Thanks

||              Paired-end : no                                               ||
||      Multimapping reads : not counted                                      ||
|| Multi-overlapping reads : not counted                                      ||
||   Min overlapping bases : 1                                                ||
||                                                                            ||
\\============================================================================//

//================================= Running ==================================\\
||                                                                            ||
|| Load annotation file gencode.v35lift37.basic.annotation.gtf ...            ||
||    Features : 62475                                                        ||
||    Meta-features : 62475                                                   ||
||    Chromosomes/contigs : 25                                                ||
||                                                                            ||
|| Process BAM file bc1-1.bam...                                              ||
||    Strand specific : stranded                                              ||
||    WARNING: Paired-end reads were found and excluded.                      ||
||    Total alignments : 44908324                                             ||
||    Successfully assigned alignments : 0 (0.0%)                             ||
||    Running time : 0.21 minutes                                             ||
||                                                                            ||

[GenoMax]

I have bam files obtain from pair-end reads.

Then you have to use -p as pointed out by @h.mon. Based on your log above that option must be missing. The stand option (-s) is exclusive of paired-end option.

Just to be sure, you did align the paired-end reads together when you generated the BAM file, correct?

Edit (Aug 2022): New versions of featureCounts require in addition to -p.

--countReadPairs    Count read pairs (fragments) instead of reads. This option
                      is only applicable for paired-end reads.

[lorenzo]

Thanks for the information! Yes, I used hisat2 with 2 read files as input.

[DareDevil]

Paired-end : no