If pileup is deprecated, so how do I get consensus?

[kentk]

I wanted to get the consensus sequence from a bam file and came across this

http://sourceforge.net/apps/mediawiki/samtools/index.php?title=SAM_FAQ

When I typed it, it cant run anymore because pileup has apparently been removed. It recommends using mpileup but theres no -c switch for consensus. Any way to do this using mpileup or should I just get an older version of samtools?

[jfk]

I wanted to get the consensus sequence from a bam file and came across this

http://sourceforge.net/apps/mediawiki/samtools/index.php?title=SAM_FAQ

When I typed it, it cant run anymore because pileup has apparently been removed. It recommends using mpileup but theres no -c switch for consensus. Any way to do this using mpileup or should I just get an older version of samtools?

I ran into this problem too as the local cluster installation of samtools is v 0.1.18. The work I am repeating/checking requires pileup and the ability to generate a consensus sequence.

Locally I have a desktop running BioLinux and installed v 0.1.17 of samtools on that. So I can do my mapping/alignment on the cluster and then generated a consensus sequence locally.

Feel free to PM me.

In the end I just installed a local

[chadn737]

Did you read the Samtools page about calling SNPs and IN/DELs with mpileup?

Since r865, it is possible to generate the consensus sequence with

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

http://samtools.sourceforge.net/mpileup.shtml

[jfk]

I would be interested in generating a consensus sequence - can mpileup do this as I'm not sure from the documentation that this is possible?

mpileup is fine for variant calling and SNP detection - the option for the -c flag in pileup being preserved in the further releases of samtools would have been useful.

[lh3]

Hasn't chadn737 already told you the answer?

[jflowers]

I think chadn737 may be asking if it is possible to use the mpileup | bcftools view | vcfutils.pl to get the consensus for MULTIPLE genomes. So far, in my reading, I only see examples of consensus sequences being generated when there is one input bam.

Is it possible to get the consensus for multiple input bams?

[chadn737]

I think chadn737 may be asking if it is possible to use the mpileup | bcftools view | vcfutils.pl to get the consensus for MULTIPLE genomes. So far, in my reading, I only see examples of consensus sequences being generated when there is one input bam.

Is it possible to get the consensus for multiple input bams?

I'm not asking anything other than whether or not the original posters bothered to read the instructions. I am pointing out to those asking the question that it says in the samtools mpileup instructions how to go about getting a consensus sequence. Pileup only worked on one bam file, mpileup allows multiple. The pipeline that is described I would think would work for multiple bam inputs given that mpileup is designed to do exactly that.

[jflowers]

Just wanted to rephrase my previous question. If say multiple different samples (e.g., humans) were sequenced and we want a separate consensus for each human, how would we generate multiple consensus sequences (one for each genome) using mpileup | bcftools view | vcfutils.pl?

[chadn737]

Just wanted to rephrase my previous question. If say multiple different samples (e.g., humans) were sequenced and we want a separate consensus for each human, how would we generate multiple consensus sequences (one for each genome) using mpileup | bcftools view | vcfutils.pl?

mpileup works equally well on one input bam as it does on multiple input bams. Thats how I typically use it, on a single bam.

So just run the pipeline for each genome/bam individually and you should get exactly what you want.

[jflowers]

Yes, mpileup works just fine on individual genomes. However, for my application, I would like to use the mpileup to call snps for many genomes simultaneously and then obtain a consensus sequence for EACH genome.

For example, I would like to run something like mpileup [options] sample1.bam sample2.bam ... sampleN.bam | bcftools view | vcfutils.pl and obtain N consensus sequences (where N is the number of input BAMs representing different sequenced individuals)

I have tried and the desired output doesnt seem to be implemented in this workflow.

I would appreciate any suggestions.

[swbarnes2]

Yes, mpileup works just fine on individual genomes. However, for my application, I would like to use the mpileup to call snps for many genomes simultaneously and then obtain a consensus sequence for EACH genome.

For example, I would like to run something like mpileup [options] sample1.bam sample2.bam ... sampleN.bam | bcftools view | vcfutils.pl and obtain N consensus sequences (where N is the number of input BAMs representing different sequenced individuals)

I have tried and the desired output doesnt seem to be implemented in this workflow.

I would appreciate any suggestions.

So you run mpileup with all the samples (though I'm not sure that doing that helps much), generate a vcf with a column for every sample. Why can't you just make individual vcfs from that to put through vcfutils?

[kentk]

Sorry chad didnt see that mpileup documentation, I read the wiki page FAQs instead which seems dated. Anyway I tried it and it worked. Thanks!

[jflowers]

swbarnes2,
Thanks for the suggestion, that just might work

[jflowers]

Hi swbarnes2,

Just wanted to follow up on the issue of generating multiple consensus sequences from mpileup run with multiple BAMs in the event that someone happend upon this thread with the same question.

I looked more carefully into your suggestion of splitting the output VCF into many VCFs (one for each genotype) and running it through vcfutils.pl

So you run mpileup with all the samples (though I'm not sure that doing that helps much), generate a vcf with a column for every sample. Why can't you just make individual vcfs from that to put through vcfutils?

A good idea, but it won't work. Looking at the vcfutils.pl vcf2fq routine, this script doesn't use any information from VCF column 10 (the column with genotypes, phred likelihoods, genotype qualities; e.g., 0/1:17,0,81:20).

All vcf2fq does is use information from the FQ and AF1 tags to call genotypes and these are set by bcftools for the set of genotypes.

Looks like creating multiple consensus sequences from multiple bams would require a currently unavailable utility.