how to extract information about mapped genes from a genome-mapping bam file

[capricy]

Hello, there,

I mapped the reads to the genome scaffold using tophat and got a bam file. I have gff file for the genome sequences. I wonder if there is an easy way to extract the information about the mapped genes and the corresponding reads.

Thank you for your help.

Capricy

[GenoMax]

Are you looking for coverage information (http://bedtools.readthedocs.org/en/l.../coverage.html) or actually looking to extract reads that are mapped to the genes (http://seqanswers.com/forums/showthread.php?t=50390)?

[capricy]

I am looking for which gene was mapped by the reads. For example:

My alignment file:
----
HWI-M01439:125:000000000-A7P33:1:1110:21257:22290 99 A_Cont998 55849 50 106M95N138M98N6M = 55887 481 ATTTTGAAGATATCGGAGTATTAGACCTCGACGCCTCACGTGAGCCAATGAGGGCTTTAGTTTGACTTCGTGTGACCTTCACCGCAGGATCAGTTGTGGAGAGGAACAGTTCCGTCACTGTGTTCTTATGCGTAGGATCAAATAACTTTTTCAATTCGCCAGATGCAGCAGCCACTTCAGCGGCCGTCTGCCCATAAAAGACGTCATCCTCCTGCAGTTCCCGAGGTTTAAGGCCAGTTTTATCATCTCT CDCEEFDFFFFFGGGGGGGGGGHHHHHHHGGGGGGHGHHHHHHGHHHHHHHHHGGHHHHHHHHHHHHHHGHGHHHHHHHHHHHGGGGGGHHHHHHHHGGHHHGGHHHHHHHHHGHHGGHHHHHHHHHHHHHGGGGGHHHHHHHHHHHHHHHHHHHHHGGGGGHHHHHHHHHHHHHHHHHHHGGGGGGGHGHGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFFFFFFFFFFFFFFFHHFFFHFFHFHF MD:Z:250 XG:i:0 NH:i:1 NM:i:0 XM:i:0 XO:i:0 AS:i:0 XS:A:-
HWI-M01439:125:000000000-A7P33:1:1110:21257:22290 147 A_Cont998 55887 50 68M95N138M98N44M = 55849 -481 CGTGAGCCAATGAGGGCTTTAGTTTGACTTCGTGTGACCTTCACCGCAGGATCAGTTGTGGAGAGGAACAGTTCCGTCACTGTGTTCTTATGCGTAGGATCAAATAACTTTTTCAATTCGCCAGATGCAGCAGCCACTTCAGCGGCCGTCTGCCCATAAAAGACGTCATCCTCCTGCAGTTCCCGAGGTTTAAGGCCAGTTTTATCATCTCTAGTAACTATTTCCGAAACGTACTCCCAACGTGGGCCTC EFBFFFFFFFFFFFFFFFFFFFFFFFEFFFFFFFFBBGFBA9.FAGFGGGGFFGGGGFBFFFGGGHHGHHEEECGGGHHHHHHHHHHFGGGGGHHHHHHHHHHHHHGGHHHGGHHGC<CHHFHHHGHHHGGHGHHGHHGGGGGGGGGHGGHHHHHHHHHGHGGGHHHHFGGGHHHHHHGGGGGGHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGGGGGHHHGGGGGGGGGGFCCDDDDDDDDD MD:Z:250 XG:i:0 NH:i:1 NM:i:0 XM:i:0 XO:i:0 AS:i:0 XS:A:-
----

my gff file indicate that: B gene is located on scaffold A_Cont998, between 558000 - 578000

I would like to have an output: HWI-M01439:125:000000000-A7P33:1:1110:21257:22290 "B gene"

I don't need to count, just match reads to genes, based on bam file and gff file

Any idea? Thanks a lot!

[cascoamarillo]

Why would you want to do that? Probably you have your reasons. One idea: I think using htseq-count (you need sam file as input) have the option of giving another sam output (--samout) which have some extra gene alignment information on it (along the flags). Then you can parse this sam output with your desire info. Good luck.

[capricy]

Ok, I will go with this idea!

Thanks