View Single Post
Old 06-30-2015, 10:33 PM   #2
SNPsaurus
Registered Vendor
 
Location: Eugene, OR

Join Date: May 2013
Posts: 521
Default

I do think the mpileup command is what you want. I am not sure what your pipeline is making, but it sounds like the pileup format (. or , for the reference base with ATCG with alternates).
You can output a vcf with genotype compared to the reference for each position:
samtools mpileup -gu -f ref.fasta sorted.bam | bcftools call -c - > all.vcf
The genotype will be 0/0 for homozygous reference, 0/1 for heterozygous with the alternate listed allele, and 1/1 for homozygous alternate.

bcftools call -cv would output the variants only.

You might want to use bcftools consensus to create a modified fasta file with your variants:
http://www.htslib.org/doc/bcftools.html#consensus
__________________
Providing nextRAD genotyping and PacBio sequencing services. http://snpsaurus.com
SNPsaurus is offline   Reply With Quote