samtools mpeileup: extract SNPs between individuals only (discard SNPs betn individuals & reference)

Question

[merodev]

I used samtools after aligning Illumina reads from > 100 samples to a reference genome and used mpileup to call variants. Now, I want to discard the SNPs that compare an individual with the reference and keep the SNPs that are between the individuals only. Is there any ways to get this done? Thanks in advance.

Answer 1

[Jorge Amigo]

you haven't used samtools neither for the alignment nor fhe variant calling, that's for sure. what you've probably done is to merge all your samples bam alignments using samtools merge and to extract non-reference discrepancies using samtools mpileup.

anyway, either if your have individual bam files or a single multisample bam file, the variant calling process rationale is the same: you have to use a variant caller software that reads as input either bam files directly or mpileup files you may have generated. a simple solution is to use the classical variant discovery pipeline (bwa for the alignment, samtools mpileup to get candidate variants, and bcftools to call variants) which, if you already have your bam file/s, it's only reduced to something like this:

samtools mpileup -uf ref.fa aln.bam | bcftools call -mv > var.raw.vcf

[merodev]

Thanks Jorge! I followed this pipeline : bowtie2 > samtools mpileup -uf ref.fa aln1.bam....aln100.bam -v > raw.vcf Now I have all variants and then used bcftools call -mv raw.vcf>called.raw.vcf. This file calls SNPs between my samples and the reference. Now I want to filter the SNPs between the individuals only (without regards to the reference). Is there any tools or filtering criteria to accomplish this?

[Jorge Amigo]

you mean that you want each sample's variants into separate files? then the proper pipeline would have been this one:

for file in aln*.bam; do
  samtools mpileup -uf ref.fa $file \
  | bcftools call -mv > ${file/.bam/.vcf}
done

if you want to extract each sample's variants from a single joined vcf file (as your called.raw.vcf should be), then you would need other solution like this one or this other one.