Get consensus sequence from blasr alignment

[guibar]

Hello,

I have some pacbio reads which I have aligned to a reference sequence using blasr. I can see the resulting alignment in igv with the help of samtools and it works fine.
What I would like is a simple way of recovering the consensus sequence that results from this alignment. I can see any simple way of doing this. Is there one?

Thanks a lot for your help?

[GenoMax]

This should do what you need: https://github.com/PacificBiosciences/pbdagcon

[flxlex]

I'm not sure the pbdagcon program is the right one. It is meant for aligning raw pacbio reads to themselves for correction. I would suggest the Quiver program: https://github.com/PacificBiosciences/GenomicConsensus. The more coverage, the better the consensus will be, though...

[guibar]

Thanks Genomax and flxlex for your replies. I will have a look at those tools.

I also realized that my problem is not specific to Pacbio, it can be rephrased as extracting the consensus from a bam alignement (regardless of the tool used to produce that alignment) and I have seen that there are a few ways of doing this.

samtools mpileup is one apparently.

http://www.researchgate.net/post/How_to_extract_the_final_sequence_after_sequence_reads_are_aligned

Does that make sense?

[rhall]

For PacBio data the highest quality consensus requires the rich QV information that is not stored in a BAM file. The workflow would be to use pbalign (which calls blasr) to generate a cmp.h5 file, then GenomicConsensus.
https://github.com/PacificBiosciences/pbalign

pbdagcon is optimized for speed, and is really intended as a tool for use in developing algorithms.

pbalign (blasr) -> GenomicConsensus (quiver) will give the most accurate consensus.

Due to the predominant indel error in PacBio data, samtools and mpileup won't give as accurate a consensus as those tools designed for PacBio data.