This question has been asked around the World Wide Web every now and then but there's never a definitive answer because the question is flawed for the following reason:

After you have already made the variant calls, in order to merge multiple variant calls into a single variant call, you need some rule (or rules) about how to do it. What do you do in situations where there are multi-allelic calls at a single position and / or there's a mixture of heterozygous and homozygous calls in your cohort? What do you ultimately report as the final GT (genotype)? Do you just add up the calls and report something like 22/11 as the GT? How do you report AD, DP, and other tags?

Even when you decide on some rules, it gets complicated in implementing it. For example, filtering on individual sample calls in VCF is something that has never really been developed so far. You can produce FILTER conditions with BCFtools on individual samples using index sub-setting, like this, but it's in no way perfect:

bcftools view -f PASS -i 'GT[3]="1/1"'

This looks at the 4th sample in your VCF and will remove that entire variant record unless it's a PASS variant and has 1/1 as genotype. In order to implement this level of filtering for what you want to do, it would require a lot of coding, but it's not impossible.

SnpSift filter also allows you to do this.

-------------------------------------------------------------

If you're literally just interested in which variant calls were made in your samples without much interest as to the individual sample GTs, then you can just strip the VCF of all sample-specific information with:

bcftools view --drop-genotypes

However, prior to doing this, I would split all multi-allelic calls and then recalculate AC/AN in the info field. This way, you would lose sample-specific GT information but you would retain which variant calls were made and also their respective frequencies in your cohort:

bcftools norm -m-any 1000Genomes.bcf -Ob > 1000Genomes.Norm.bcf

bcftools view 1000Genomes.Norm.bcf | /Programs/bcftools-1.3.1/vcfutils.pl fillac > 1000Genomes.Norm.AC.vcf 

awk -F"\t" '/^#CHROM/{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8}; !/^#/{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8}' 1000Genomes.Norm.AC.vcf

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO
1   10177   1:10177:10177:A:AC  A   AC  .   PASS    AC=2130;AN=5008
1   10235   1:10235:10235:T:TA  T   TA  .   PASS    AC=6;AN=5008
1   10352   1:10352:10352:T:TA  T   TA  .   PASS    AC=2191;AN=5008
1   10616   1:10616:10616:CCGCCGTTGCAAAGGCGCGCCG:C  CCGCCGTTGCAAAGGCGCGCCG  C   .   PASS    AC=4973;AN=5008
1   10642   1:10642:10642:G:A   G   A   .   PASS    AC=21;AN=5008
1   11008   1:11008:11008:C:G   C   G   .   PASS    AC=441;AN=5008
1   11012   1:11012:11012:C:G   C   G   .   PASS    AC=441;AN=5008
1   11063   1:11063:11063:T:G   T   G   .   PASS    AC=15;AN=5008
1   13110   1:13110:13110:G:A   G   A   .   PASS    AC=134;AN=5008
1   13116   1:13116:13116:T:G   T   G   .   PASS    AC=486;AN=5008
1   13118   1:13118:13118:A:G   A   G   .   PASS    AC=486;AN=5008
1   13273   1:13273:13273:G:C   G   C   .   PASS    AC=476;AN=5008
1   13284   1:13284:13284:G:A   G   A   .   PASS    AC=7;AN=5008
1   13380   1:13380:13380:C:G   C   G   .   PASS    AC=41;AN=5008
1   13483   1:13483:13483:G:C   G   C   .   PASS    AC=10;AN=5008
1   13494   1:13494:13494:A:G   A   G   .   PASS    AC=7;AN=5008
1   13550   1:13550:13550:G:A   G   A   .   PASS    AC=17;AN=5008
1   14464   1:14464:14464:A:T   A   T   .   PASS    AC=480;AN=5008
1   14599   1:14599:14599:T:A   T   A   .   PASS    AC=739;AN=5008
1   14604   1:14604:14604:A:G   A   G   .   PASS    AC=739;AN=5008
1   14930   1:14930:14930:A:G   A   G   .   PASS    AC=2415;AN=5008
1   14933   1:14933:14933:G:A   G   A   .   PASS    AC=142;AN=5008
1   15211   1:15211:15211:T:G   T   G   .   PASS    AC=3050;AN=5008
1   15245   1:15245:15245:C:T   C   T   .   PASS    AC=21;AN=5008
1   15274   1:15274:15274:A:G   A   G   .   PASS    AC=1739;AN=5008

At least from this you retain the VAF (variant allele frequency) and the variant calls that were made.

Some other related scripts that you may find useful:

• A: calculate Per variant Heterozygosity from VCF file
• A: How to get sample names and genotype for SNP in multi-sample VCF file

------------------------------------------------

The more robust way would be to concatenate your reads together and then re-do the entire pipeline, but you imply that this may not be possible. In such a situation, I would 'sample' reads from the BAM at multiple frequencies and, ironically, call variants on each subset of reads. In your final variant list, you then just take the first FORMAT column's values as the representative set. I do this here: https://github.com/kevinblighe/ClinicalGradeDNAseq

There may very well be other solutions!

Kevin