This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Re-generating rnaSTAR reported values for "Number of input reads" and "Uniquely mapped reads number"

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

rna-seq star
• 2,571 views • link
updated by Ram Ram • written by noushin.farnoud noushin.farnoud

1 answer

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.

• 0 views • link
updated by Ram Ram • written by Devon Ryan Devon Ryan

Log in to answer this question.

More posts like this