Re-generating rnaSTAR reported values for "Number of input reads" and "Uniquely mapped reads number"

Question

[noushin.farnoud]

Hello,

I am trying to query of the rnaSTAR Aligned.sam output using the SAM flags to re-generate Number of input reads and Uniquely mapped reads number that are reported in Log.final.out for paired end RNASeq data. But so far, I had no success in getting the exact STAR reported numbers.

As for "Number of input reads", I am using:

samtools view -S -f 2 -F 8 -c Aligned.out.sam

and then dividing this number by two, which seems to be still more than what Log.final.out has reported. I assume if I can get how this number is calculated, I can then use grep by NH:i:1 tag to get the number of unique reads that STAR reports (unless you know otherwise).

I greatly appreciate your input.

Noushin

Answer 1

[Devon Ryan]

I wouldn't necessarily expect that samtools command to match either of the numbers you mentioned. If you want to match the total number of input reads, you need to:

1. Ensure that unmapped reads are included in the output
2. Ensure that multimappers are
   1. counted only once or...
   2. reported only once

Then a simple samtools view -Sc foo.sam should produce the total number of input reads (divide by 2 for read pairs).

For the "uniquely mapped" reads, a simple samtools view -Scf 2 foo.sam will work if and only if STAR doesn't produce singleton alignments. Whether it will do so is dependent upon the settings you give.