The Differences Of The Number Of Reads From Command Line And Rna-Seqc

Question

[mad.cichlids]

Hi, all

I am trying to use RNA-SeQC to evaluate my RNA seq results. In its output from my pair end sequences:

Sample    Note     Total Purity Filtered Reads Sequenced    Alternative Aligments        
jmp         qc               10,846,940                             8,148,101

Which gave me a total number of reads = 10,846,940 + 8,148,101 = 18,995,041

However, when I use command lines to count:

zcat sub592_1.fastq.gz | echo $((`wc -l`/4)) #5924954


zcat sub592_2.fastq.gz | echo $((`wc -l`/4))#5924954

which totals up to 5924954x2 = 11 849 908

The differences equals 18,995,041 - 11,849,908 = 7 145 133

Could anyone please explain what caused this discrepant result ? Thank you very much !

[P]

I have the same question! Did you find an answer for this?

~Rini

Answer 1

[madk00k]

Hi!

RNAseq-QC does not count the reads that were not mapped. One can output the number of mapped reads in BAM file using samtools flagstat command or Qualimap tool. Additionaly, as I suppose, one has to exclude number of reported rRNA reads to get the number reported by RNASeq-QC.