Using the recent versions of SAMtools and BCFtools, you can do it like this:

1, call variants

samtools mpileup --redo-BAQ --min-BQ 30 --per-sample-mF --output-tags DP,AD -f "${Ref_FASTA}" \
--BCF BAM1.bam BAM2.bam BAM3.bam BAM4.bam | \
bcftools call --multiallelic-caller --variants-only -Ob > Out.bcf ;

Note that here I have 4 replicates, which samtools mpileup accepts.

For less stringent calling (outputs pretty much anything that can possibly be reported, then pipe it into:

bcftools call --consensus-caller --variants-only --pval-threshold 1.0 -Ob

2, split multi-alleles (1st pipe) and left-align indels (2nd pipe)

bcftools norm -m-any Out.bcf | bcftools norm -Ob --check-ref w -f human_g1k_v37.fasta > Out.norm.bcf ;
bcftools index Out.norm.bcf ;

3, add all available tags to the BCF

See my answer here: A: How to use bcftools to calculate AF INFO field from AC and AN in VCF?

4, further filter the variants

#Options:
#   -Q INT  minimum RMS mapping quality for SNPs [10]
#   -d INT  minimum read depth [2]
#   -D INT  maximum read depth [10000000]
#   -a INT  minimum number of alternate bases [2]
#   -w INT  SNP within INT bp around a gap to be filtered [3]
#   -W INT  window size for filtering adjacent gaps [10]
#   -1 FLOAT    min P-value for strand bias (given PV4) [0.0001]
#   -2 FLOAT    min P-value for baseQ bias [1e-100]
#   -3 FLOAT    min P-value for mapQ bias [0]
#   -4 FLOAT    min P-value for end distance bias [0.0001]
#   -e FLOAT    min P-value for HWE (plus F<0) [0.0001]
#   -p      print filtered variants
echo "Applying further filtering to called variants..." ;
bcftools view Out.norm.bcf | \
misc/vcfutils.pl varFilter -d "${ReadDepthCutoff} \
-w 1 -W 3 -a 1 -1 0.05 -2 0.05 -3 0.05 -4 0.05 -e 0.05 -p \
> Out.norm.filt.vcf  ;

echo "Variants filtered out:"

Note that, in the final part (#4), I use a perl executable (vcfutils.pl) that comes bundled with BCFtools in the misc directory.

Kevin