Calculate Number of Reads per class code (cufflinks)

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

cuffcompare

• 2,405 views • link

updated 3.9 years ago by Ram

• written 10.9 years ago by guisa.santos •

1 answer

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

• 0 views • link

written 10.9 years ago by michael.ante

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?

• 0 views • link

updated 3.9 years ago by Ram

• written 10.9 years ago by guisa.santos •

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

• 0 views • link

written 10.9 years ago by michael.ante

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

• 0 views • link

written 10.9 years ago by guisa.santos •

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

• 1 views • link

written 10.9 years ago by michael.ante

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

• 1 views • link

written 10.9 years ago by guisa.santos •