featuresCount -s 2 option

Question

[noeD]

Hello! :)

I have paired-end RNA-seq data aligned by STAR, I used featuresCount in order to get read summarization.(I know that STAR can give the count, but I would like to try with featuresCount). I don't understand why if I use the "-s 2" option I obtained a lower number of assigned read than without it.

This is the output of STAR:

 Mapping speed, Million of reads per hour | 66.37

                      Number of input reads |   20999464
                  Average input read length |   202
                                UNIQUE READS:
               Uniquely mapped reads number |   18063283
                    Uniquely mapped reads % |   86.02%
                      Average mapped length |   200.30
                   Number of splices: Total |   4054756
        Number of splices: Annotated (sjdb) |   0
                   Number of splices: GT/AG |   4012566
                   Number of splices: GC/AG |   28175
                   Number of splices: AT/AC |   1578
           Number of splices: Non-canonical |   12437
                  Mismatch rate per base, % |   0.31%
                     Deletion rate per base |   0.02%
                    Deletion average length |   1.36
                    Insertion rate per base |   0.01%
                   Insertion average length |   1.34
                         MULTI-MAPPING READS:
    Number of reads mapped to multiple loci |   1876550
         % of reads mapped to multiple loci |   8.94%
    Number of reads mapped to too many loci |   40430
         % of reads mapped to too many loci |   0.19%
                              UNMAPPED READS:
   % of reads unmapped: too many mismatches |   0.00%
             % of reads unmapped: too short |   4.79%
                 % of reads unmapped: other |   0.07%
                              CHIMERIC READS:
                   Number of chimeric reads |   0
                        % of chimeric reads |   0.00%

This is the output of featuresCount with the following command:

featureCounts-T 5 -t exon -g gene_name -a /Homo_sapiens.GRCh38.89.chr.gtf -s 2 -o output align.bam

Assigned    8798522
Unassigned_Unmapped 0
Unassigned_MappingQuality   0
Unassigned_Chimera  0
Unassigned_FragmentLength   0
Unassigned_Duplicate    0
Unassigned_MultiMapping 11493680
Unassigned_Secondary    0
Unassigned_Nonjunction  0
Unassigned_NoFeatures   26603481
Unassigned_Overlapping_Length   0
Unassigned_Ambiguity    724563

And this is the output without using the -s 2 option:

Assigned    16492987
Unassigned_Unmapped 0
Unassigned_MappingQuality   0
Unassigned_Chimera  0
Unassigned_FragmentLength   0
Unassigned_Duplicate    0
Unassigned_MultiMapping 11493680
Unassigned_Secondary    0
Unassigned_Nonjunction  0
Unassigned_NoFeatures   17636623
Unassigned_Overlapping_Length   0
Unassigned_Ambiguity    1996956

Thank you in advance!

Answer 1

[Devon Ryan]

Apparently you have a non-directional library.

[noeD]

Thank you for your reply. How can I confirm it? Furthermore, is it normal to obtain that percentage of assigned reads? Thank you in advance for your help!

[Devon Ryan]

Ask whoever did the library prep., they can tell you whether it was directional or not. Your percentage assigned is rather low, but there can be many reasons for this, such as sequencing total RNA rather than doing polyA-enrichment.

[noeD]

I apologize for my stupid question.. but if I have paired-end reads, it doesn't mean that the library is directional?

Thank you in advance

[GenoMax]

No it does not. dUTP based method is described here.