RNA-seq & gene annotation

[rongronghai]

I have RNA-seq data and need the count of reads at gene level. To do this using BEDTOOLS I need an gene annotation file. I download the hg18 (group: Genes and Gene Prediction Tracks ; track: UCSC Genes) from UCSC in BED format. However, it seems that each line in this file represent an isoform instead of a gene (see the two lines below have the same position but different IDs).

chr7 50625258 50817758 uc003tpk.1 0 - 50628142 50767513 0 19
chr7 50625258 50817758 uc010kzb.1 0 - 50628142 50709814 0 18

What I need actually is each line for a gene. How can I get that? Thanks~

[RockChalkJayhawk]

use the mergeBed feature in BEDtools to give yourself a non-redundant list.

[rongronghai]

use the mergeBed feature in BEDtools to give yourself a non-redundant list.

Thanks for your reply~ I have checked the mergeBed, but I am a little worried about the merged "feature" it gives. It may merge isoforms originated from different genes, right? If so, the resulting "feature" may have no biological meaning.

[RockChalkJayhawk]

Thanks for your reply~ I have checked the mergeBed, but I am a little worried about the merged "feature" it gives. It may merge isoforms originated from different genes, right? If so, the resulting "feature" may have no biological meaning.

Not if you use the strand option.

[epigen]

You could also use htseq-count from the HTSeq package (developer Simon Anders is active on this forum) with a GTF file, e.g. from ENSEMBL. I constructed a fake GTF file from the RefSeq genes with some non-professional Perl scripts that take into account that different genes can overlap (on the same strand!). There may be a professional conversion tool out there or the possibility to download UCSC genes in GTF format.