Reference Mapping To Consensus Fasta

Question

[Juliofdiaz]

I have mapped illumina reads to a reference fasta. I can visualize the map using the tview tool of samtools or CLC genomics. I already have the bam file of the mapping and it is sorted as well. I also have the mpileup file. Is there a way to get the "assembled" reads in a fasta format? For example if my reference and the mapping looks like this

    ref ATGCATGC
        ...A
         ..A...
           A....
             ...

and position 4 is recognized as a SNP. Can I get a fasta that looks like this:

>result
ATGAATGC

Also will it take gaps into account? Thanks

[Zev.Kronenberg]

Why do you want to do this? What analyses would you use this for? Going from SNVs to FASTA requires some assumptions that could bias subsequent analyses.

[Juliofdiaz]

I wanted to BLAST MLST genes to the consensus sequence. Do you think its better to map my reads to the sequence of the MLST genes?

[bioinfo8]

@Juliofdiaz Did you find the solution? I have a similar problem posted here.

Answer 1

[David Quigley]

As far as I can see, samtools has you covered. From the man page, under "Generate the consensus sequence for one diploid individual:"

samtools mpileup -uf ref.fa aln.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq