What Better Way To Get Consensus Sequence Of Mapping Against Reference Genome Using Samtools Mpileup?

Question

[thiagomafra]

Hi people,

I have a set of reads mapped against a reference genome and need take consensus sequence using samtools mpileup. I used the command:

samtools mpileup -AuDS -Q 10 -q 10 -d 10000 -f <reference.fasta> <file.sorted.bam> | bcftools view -gc - | vcfutils.pl vcf2fq | perl fastq2fasta.pl cns.fasta &

Someone uses other parameters in samtools mpileup?

[Rm]

check this similar post: How to generate a consensus fasta sequence from SAM tools pileup?

Answer 1

[Jorge Amigo]

you could also try not using samtools but GATK instead, which allows selecting better calling algorithms than samtools (UnifiedGenotyper and HaplotypeCaller). the underlying idea would be first to perform the variant calling, and then to generate the fasta file:

java -jar GenomeAnalysisTK.jar -R hg19.fa -T HaplotypeCaller -I sample.bam -o sample.vcf
java -jar GenomeAnalysisTK.jar -R hg19.fa -T FastaAlternateReferenceMaker -o sample.fa --variant sample.vcf

Answer 2

[Joseph Hughes]

I use the following command in samtools Version: 0.1.18:

samtools mpileup -E -uf <reference.fasta> <file.sorted.bam> | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fastq

The -E option helps with sensitivity especially if you have sites with multiple nucleotide polymorphisms.