Annotating bed/bam from GFF

[Pratap, Abhishek]

Hi Aaron

 

I have superficially followed the toolkit until now.  Before I ask my question I will admit I haven’t read the manual in detail but  want to know if doing  the following is possible at all.

 

I have a bed/bam  file/s as  result of 454 reads alignment. The  annotation file  corresponding to this sample  is in GFF format.  My question will it be possible to annotate the bed/bam   based on GFF. I know this could be done with a script but before I go that route I wanted to check if there already exists a slick way of doing it in this toolkit.

 

Thanks!

-Abhi

[Aaron Quinlan]

Hi Abhi,

   

intersectBed has an "-abam" option that allows one to find overlaps between alignments in BAM format and annotations in either BED or GFF format.  However, only the GFF annotation will be reported in GFF format.  The alignments will be in BED format.

Here's an example with some 454 data I have:

$ intersectBed -abam 454.bam -b genes.gff -wa -wb -bed | head

chr16 73795015 73795502 A12-F4JXVRN01EV87A 57 + chr16 hg18_knownGene exon 73795495 73795721 0.000000 - . gene_id "uc002fdr.1"; transcript_id "uc002fdr.1"; 

chr16 73795015 73795504 A06-F51QE6Q01CWCKM 57 - chr16 hg18_knownGene exon 73795495 73795721 0.000000 - . gene_id "uc002fdr.1"; transcript_id "uc002fdr.1"; 

chr16 73795015 73795503 A06-F51QE6Q01EQTVD 57 - chr16 hg18_knownGene exon 73795495 73795721 0.000000 - . gene_id "uc002fdr.1"; transcript_id "uc002fdr.1"; 

chr16 73795015 73795500 A08-F4JXVRN01A75HJ 58 + chr16 hg18_knownGene exon 73795495 73795721 0.000000 - . gene_id "uc002fdr.1"; transcript_id "uc002fdr.1"; 

chr16 73795017 73795500 A01-F4JXVRN01EMQZE 57 + chr16 hg18_knownGene exon 73795495 73795721 0.000000 - . gene_id "uc002fdr.1"; transcript_id "uc002fdr.1"; 

chr16 73795019 73795504 A07-F51QE6Q01DAQEA 54 + chr16 hg18_knownGene exon 73795495 73795721 0.000000 - . gene_id "uc002fdr.1"; transcript_id "uc002fdr.1"; 

chr16 73795044 73795518 A05-F4JXVRN01B6IDU 57 - chr16 hg18_knownGene exon 73795495 73795721 0.000000 - . gene_id "uc002fdr.1"; transcript_id "uc002fdr.1"; 

chr16 73795044 73795502 A06-F51QE6Q01C8C0J 55 - chr16 hg18_knownGene exon 73795495 73795721 0.000000 - . gene_id "uc002fdr.1"; transcript_id "uc002fdr.1"; 

chr16 73795046 73795532 A08-F51QE6Q01D5TLX 58 - chr16 hg18_knownGene exon 73795495 73795721 0.000000 - . gene_id "uc002fdr.1"; transcript_id "uc002fdr.1"; 

chr16 73795047 73795522 A06-F4JXVRN01C2WVV 50 - chr16 hg18_knownGene exon 73795495 73795721 0.000000 - . gene_id "uc002fdr.1"; transcript_id "uc002fdr.1"; 

The first six columns of each line are the BAM alignment in BED format and the remaining fields are the GFF annotation that the alignment overlapped.  I do not have an option to also write the BAM alignment in GFF format.  Note that you must use the "-bed" optio if you want to write BED output when using BAM input.  The default is to write BAM.  That is, the default is to write BAM output for those alignments that overlap your annotations.  This allows you to create subsets of your alignments based on annotations.

I hope this helps.

Aaron

[Pratap, Abhishek]

Thanks for a detailed reply Aaron. I tried the same thing on my dataset and it is hard for me to believe that there is no overlap in the reads to gene in gff or bed file. If you have some time let me know, I will be happy to share the dataset with you offline.

Thanks!
-Abhi
________________________________________
From: bedtools...@googlegroups.com [bedtools...@googlegroups.com] On Behalf Of Aaron Quinlan [aaronq...@gmail.com]
Sent: Tuesday, February 02, 2010 9:43 PM
To: bedtools...@googlegroups.com
Subject: Re: [bedtools-discuss] Annotating bed/bam from GFF

[Aaron Quinlan]

Hi Abhi,
Please make sure you annotation file is tab-delimited and that the chromosome labels in the annotation file match those in the BAM file (e.g., chr1 in both, not "1" and "chr1"). If this holds up, I'd be happy to look at a small dataset for you.

Aaron

[Pratap, Abhishek]

Hi Aaron

You got my error. It was due to the mismatch in chr label. Works fine now.

Thanks for releasing what looks like a robust toolset. I will learning more about it in the coming weeks.

Cheers,
-Abhi