This is a test version of Biostars. For the public version, visit https://www.biostars.org.
What Better Way To Get Consensus Sequence Of Mapping Against Reference Genome Using Samtools Mpileup?

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?

consensus mpileup reference
• 13,377 views • link
updated by Jorge Amigo Jorge Amigo • written by thiagomafra thiagomafra

2 answers

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
• 0 views • link
written by Jorge Amigo Jorge Amigo

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.

• 0 views • link
written by Joseph Hughes Joseph Hughes

Log in to answer this question.

More posts like this