How To Convert My Denovos To Gff/Bed

Question

[upendrakumar.devisetty]

Hi, I have recently done denovo assembly and now i want to convert the denovos into either gff or bed format. This is to both to visualize them on Genome browser and also for quantification purposes (like cufflinks). I found that exonerate model est2genome can do the job and so i have used that but unfortunately the exonerate is not the correct job it seems.

Is there any other alternative exists besides exonerate?

Thanks Upendra

Answer 1

[Alex Reynolds]

If you can convert your assembly into one of the formats supported in BEDOPS conversion scripts, such as BAM or SAM, you could use those scripts to turn converted assembly data into BED.

Update

I get more than six columns when running bam2bed on your BAM file:

$ bam2bed < velvet_trinity_tophat_comb_novel_final_transcripts.fasta.cap.contigs.transdecoder.bam | head
A01    158188    158533    Contig545|m.1503    0    +    229    123M3D66M3I153M    *    0    0    ATGGAGGTGATTGATTTTTCGGAGACGGAGATGGAAGCAGCGGAGCAGCTAGTGCAATTGAGCGAGGATGATACGTTGAGCTGTAGCAGCGGTACAGGCTTGAGCGTTAGCGGTTGTGAAGGCGGATTCCATAGGAAGAGACACAGCGATGTTATTAGTGATGAGGTGCAAAACGACGGCGTTGTACGTACGACGATGAATAACAATGCTACAGACGCTCAATGTTTTGTGAAGGCAATCACGGAGACGAACATAATAAGGAGGAGGTATAAGAAGAAAAAGTTTCGGTCTCTTGCGAGTCTTTACAGGGCAACGAAGGAGATGACGACGGACGAGTTGATCTAG    *    AS:i:276    XS:i:0    XF:i:0    XE:i:7    NM:i:17
A01    160347    160578    Contig1727|m.4383    16    -    138    231M    *    0    0    TTAATCCATTTGTTTCACTTTGTTAGCAGTCTCAATTGGTTTCTTGCTGATCATTTCATGATGGGATATATCTGTGTGGAATCTTTGTTTTGGTTTGACTTTTTGGTCGATGATGATTGACGTTAAGGCTTGTTTCTCAAAACATTCGCAGTGGCCGCCTAAGGCTTTCGCCTTTCGCCGGGCTTGAGAGCTTATGATGTTTTTATCCTTCCTTGGATGTCTGAGACCAAT    *    AS:i:231    XS:i:107    XF:i:0    XE:i:6    NM:i:0
A01    888838    889347    Contig1477|m.3679    16    -    92    509M94S    *    0    0    TTATAAAAAAGCAATACTATCATGGCGTTGGGGCGAGAAAGAGGCGCAACACATGCCTGTGAGTCCATGACGGCTCCAGAACTCTTCAGTATCAGCGACTCCGTGGATCCTGACTCTTCTCGTGCTCTTCGTGTCTGTGTTGAAGTAGAAAACATAGAAGGGCTGGGCAGGGCCTGGCAGAGATGTTGGGGAGAAAATGATTTCACCGGCCTTGTTGGTTCCTTGGGAGATCATACGCTGACTCCTCCGGAAAGCGAAGGGAAGAGACAAGGAAAGCTGAAATGATTGTCTGGACCACTCATGTTTGGTCACATCCTCCAGTATCCAAAGATCAAAACGGTCGAAAGAAGAACCAGTTCTCAAACAGTTTGTTGGTACAAGGACATCTAGCTTCCCTTTGTATTCTATCAATGATGTATAAACCTTCCAAACCAGGACCTCCTTAGGCGTAGTGATAAAACTTAGTATCCTCTCATATCTTACATCAAAGCACACAAACACAGGAGGAGTTTGATTTGGAACAGAAGCAGCATAATAGATGAAACCATTGATGCTTAGTCCCACAGTCATAGGGCAATAATACGGGGAGGTTCCTTGCGTTTC    *    AS:i:449    XS:i:293    XF:i:0    XE:i:8    NM:i:15
A01    888864    889123    Contig1477|m.3680    0    +    68    5S259M    *    0    0    ATGGCGTTGGGGCGAGAAAGAGGCGCAACACATGCCTGTGAGTCCATGACGGCTCCAGAACTCTTCAGTATCAGCGACTCCGTGGATCCTGACTCTTCTCGTGCTCTTCGTGTCTGTGTTGAAGTAGAAAACATAGAAGGGCTGGGCAGGGCCTGGCAGAGATGTTGGGGAGAAAATGATTTCACCGGCCTTGTTGGTTCCTTGGGAGATCATACGCTGACTCCTCCGGAAAGCGAAGGGAAGAGACAAGGAAAGCTGAAATGA    *    AS:i:227    XS:i:167    XF:i:3    XE:i:4    NM:i:8 
A01    1167009    1167386    Contig1598|m.4011    16    -    236    377M4S    *    0    0    CGCCGTGATAGGAATGATCGCAGCAGTGGAGGGCTTCGACTGTAGCTTGGTCTCATGATTGAGTAGCTTCTCGTGAACCTCCGTGATCGAGGGAGGATGATCGCGGCCTTCAATCTGATCAACGATCTGTTTGTATTCCTCCGGCAGGCCCTCCAGTATGTATTCAATCTGATCTTCAATGTCATATGGTTTCCCAAGGAGAGCCAACTGATCAAATCGAGTGGTAAACCCCTGAACGTAGGCATCGATAGTCATCGTGCCTTTGGACCAGTTTTTGACTTGCTGCCGAAGTTGTTGGATATGACCTCGACTCGGTTTTGCATACGTGCTGGACAAGATACTCCAGATCTGGGCCGAAGATGTTGCAGTGGAGAGGACTCT    *    AS:i:377    XS:i:0    XF:i:0    XE:i:9    NM:i:0
...

I'm not sure what's wrong with your setup. Maybe check the version of samtools you are using, and upgrade, if need be. Here is the version I am running locally:

$ samtools --version
samtools 0.2.0-rc5-27-g1b408f7
Using htslib 0.2.0-rc5-43-g650202c
Copyright (C) 2013 Genome Research Ltd.

[upendrakumar.devisetty]

I have seen them before but i haven't tried any one those. I will try the bam2bed now....Thanks for help anyway...

[upendrakumar.devisetty]

Hi Alex, I have managed to convert my contigs to bam first and then with the help of BEDOPS tool bam2bed i have converted it to bed. But i have only got 6 columns in the bed file. How can i get the rest of the columns so that i can get the exon and intron structure? Any ideas of how to do that?

[Alex Reynolds]

I'm not sure how you got six columns. The bam2bed script should give you (at least) 13 columns, assuming that your BAM file has those records; see the column-mapping part of the bam2bed documentation. If you think you've run into a bug, feel free to post your BAM file somewhere public and I'll take a look.

[upendrakumar.devisetty]

Just sent the link to the file through PM....

[Alex Reynolds]

Please see the update to my answer.

[upendrakumar.devisetty]

Thanks Alex. You think problem might have occurred during converting sam to bam then? Also the latest samtools version is samtools-0.1.19 and how can you get 0.20? This is what i did before trying bam2bed

bwa bwasw -t 12 ../db/B.rapa_genome_sequence_0830 velvet_trinity_tophat_comb_novel_final_transcripts.fasta.cap.contigs.transdecoder.cds > velvet_trinity_tophat_comb_novel_final_transcripts.fasta.cap.contigs.transdecoder.sam
samtools view -bS velvet_trinity_tophat_comb_novel_final_transcripts.fasta.cap.contigs.transdecoder.sam > velvet_trinity_tophat_comb_novel_final_transcripts.fasta.cap.contigs.transdecoder.bam 
bam2bed < velvet_trinity_tophat_comb_novel_final_transcripts.fasta.cap.contigs.transdecoder.bam > velvet_trinity_tophat_comb_novel_final_transcripts.fasta.cap.contigs.transdecoder.bed

Am I doing anything wrong. Thanks

Upendra

[Alex Reynolds]

I don't think conversion from SAM to BAM would affect anything, but it is a superfluous step as you can just run sam2bed on the SAM output directly. There's no need to convert from SAM to BAM unless you need a version of these reads in BAM.

To get the latest version of samtools, you might do something like this:

$ git clone https://github.com/samtools/samtools.git
$ git clone https://github.com/samtools/htslib.git
$ cd samtools
$ make
...
$ cp /usr/local/bin/samtools /usr/local/bin/samtools.backup
$ cp samtools /usr/local/bin/samtools

This would depend on your setup and where you keep your copy of samtools.

Documentation for sam2bed is available here.

You might also make sure you're running BEDOPS 2.4 and all updated components. I don't think any older versions of conversion scripts were ever limited in number of output columns (we were always "non-lossy" in conversion), but it can't hurt to check that you're running the latest stuff.

If you need a compressed version of your reads, consider converting to the BEDOPS Starch format with the sam2starch script, which stores data in a fraction of the size of a BAM file and offers random-access features. BEDOPS tools operate on Starch archives, as well as common BED.