Calculate Number of Reads per class code (cufflinks)

Question

[guisa.santos]

Dear all,

Does anybody know how to calculate the total number of reads "used" by each transcript class code (=, c, j, i, etc) from cuffcompare?

Thank you in advance.

Belisa

Answer 1

[michael.ante]

Hi Guisa,

If you need the number of reads I would avoid cufflinks. Anyway you can use awk's associative array to make an approximation:

awk 'NR>1{class[$3]+=($10*$11)}END{for(i in class){print i"\t"class[i]}}' cuffcompare.tmap

You add per line the average read-depth times the length to the value of each class.

Cheers,

Michael

[guisa.santos]

Thank you Michael. This is exactly what I have done (also using awk :) So now I know that my reasoning is correct.

(?) My issue with this approach is that the sum of the reads estimated this way add-up to more than the total number of mapped reads... So I do do not know how to discuss this difference: my total number of aligned pair reads is 1,259,591,756 and the sum of all estimated reads (from all class codes from all tmap files) is 3,732,338,824.4 - these values are quite different). Do you have any clue why this is?

[michael.ante]

Argh; I think you have to divide by the read length / alignment length.

[guisa.santos]

Yes, i know, these values are already divided by the read length (in my case 101 bp)...

[michael.ante]

You can use htseq-count with the combined gtf as reference counting the uniquely mapping reads per class_code:

htseq-count ... -t class_code my.bam cufflinks.combined.gtf

[guisa.santos]

Great suggestion, thanks! I will give it a try (and then post the results :)