Annotating a SAM file - best method?

[Cirno]

Hey there,

So far, I have put 76mers through tophat, then did BAM to SAM conversion. I looked at the Tophat page, and was happy to find IGV, etc but this doesn't seem to let me DO anything with the data, just look at it. If there was a way to dump what you see in IGV across the entire genome (i.e. reads annotated to genes and gene features) that would just be too convenient.

So yeah I need to annotate the SAM file. The only option that seems good is HTSeq. The overall goal of the study is levels of alternative splicing; will HTseq be good enough to neatly quantify every type of event? Is there some other method out there?


Thanks in advance.

[colindaven]

I use Bedtools for annotation, I think using BAM files and GFF features. It's nice and quick.

Otherwise Galaxy is good too for those who don't enjoy the command line.

[Cirno]

Thanks for the response.

What part of Bedtools? I don't see any part of it that would do this. :|

What I want to do is find out which genes my reads correspond to. There are many overlapping genes in arabiopsis, and I want reads that can't be assigned to one of these to be excluded - this is why HTseq looked good. I do not want to see which exons they map to really either, because I am looking for reads that are mostly or exclusively in introns.

I have tried with HTSeq, but it almost always returns errors!
I got my GFF file from TAIR, and this happens:

htseq-count -t gene -m union my_reads.sam TAIR10_GFF3_genes.gff
Error occured in line 2 of file TAIR10_GFF3_genes.gff.
Error: Feature AT1G01010 does not contain a 'gene_id' attribute
[Exception type: SystemExit, raised in count.py:55]

I really don't get it. The instructions it gives say to use GFF files, but only GTF files have this....I am really really confused.

[swbarnes2]

If you get a .bed file of features (exons) and coordinates, you can use BEDTools intersect to figure out what reads cross what genes. Or rather, what reads don't cross exons.

But I'm not sure how to get a .bed file of exons. The format is simple; it's just feature name, start, and stop coordinate, tab-delimited.

[Orr Shomroni]

swbarnes2, BEDtools has a simple bamToBed command that converts BAM files to BED format. And I guess to annotate, I need to intersect the converted BED file with a BED file of genes for that genome, right?

P.S.: BEDtools commands can be found here:

http://code.google.com/p/bedtools/wiki/Usage

[dariober]

So yeah I need to annotate the SAM file. The only option that seems good is HTSeq. The overall goal of the study is levels of alternative splicing; will HTseq be good enough to neatly quantify every type of event? Is there some other method out there?

If by annotating a sam file you mean assigning reads to genes to determine the level of expression of each transcript, I think what you can use is the cufflinks suite.

You can write a batch script to dump images from IGV but doing it for the whole genome would result in hundreds of thousands of image files!

Good luck
Dario