post-assembly genome analysis workflow question

[tom_mlvs]

Hi All,

We want to extract sequence information (for various genes) from a number of genome assemblies and generate consensus sequences for comparison between genomes representing different experiments.

What we have been doing is using samtools to extract regions from the genomic bam file, then trying to convert those into fasta format using bam2fastq. Everything we've extracted has been groups of overlapping short reads, we have not been successful at obtaining consensus sequences.

Is there an alternative workflow that would be more efficient/better? Are there suggestions for tools we should be using instead of/in addition to samtools and bam2fastq?

(Note: We have tried using the samtools programs (mpileup, bcf view, and vcfutils.pl) to generate a consensus sequence. Unfortunately, the (pipelined and non-pipelined) use of the program ‘bcftools view’ generates the following error: [bcf_sync] incorrect number of fields (0 != 5) at 0:0)).

[swbarnes2]

You are trying to make a vcf with bcftools? That should work, there's probably something off about your input file.

Do you mean genome assemblies, or alignments? I assume alignments, since you have bam files?

I suppose you could try using samtools view to pull sections of your .bam, and then you could put those .bams through velvet, to assemble a consensus for that sample at that region.

But I think getting the vcf files for those regions is the way to go. You just aren't doing it right.

[tom_mlvs]

We've been trying to pull out a chromosome, and go at it that way. Should be be doing it gene by gene instead, or is there something obviously wrong here:

samtools view -b exome_input.bam 12 > chr12.bam
samtools mpileup -f chr12.fa chr12.bam > chr12_pileup
bcftools view -cg chr12_pileup > chr12_vcf


edit:
It seems we're not the only ones to get the following error
[bcf_sync] incorrect number of fields (0 != 5) at 0:0))

Maybe the suggestions here (nohup) will fix it:

http://seqanswers.com/forums/showthread.php?t=15881

Will update when we've tried.

Thanks