You want to find a set of SNPs in populations of two species, for each of which you have a distinct reference genome. You then want to know which of those SNPs are shared between the two genomes (there is variation in each species), and which ones are private to only one of the species (in the other genome there is no variant at that location).
To accomplish this task, you need to first find SNP variation in each genome, and then figure out which SNPs are the same between the two genomes. In other words, line up two VCF files for different genomes and figure out which SNPs are the same between them. In comparing genome assemblies, this is sometimes called a liftover operation.
Genomes can have different coordinate systems.
Let's first think about what exactly you are looking at with 2 different reference genomes. For example, here are outputs of
samtools faidx for 2 different S. cerevisiae strains, S288C and Kyokai-7:
ref|NC_001133| 230203 95 60 61
CHR2.19980913 813142 234173 60 61
CHR3.19980521 315339 1060907 60 61
CHR4.19990210 1531929 1381540 60 61
CHR5.19970727 576870 2939039 60 61
CHR6.19970727 270148 3525562 60 61
CHR7.19970703 1090936 3800252 60 61
CHR8.19970727 562638 4909411 60 61
ref|NC_001141| 439885 5481523 60 61
CHR10.19970727 745440 5928778 60 61
CHR11.19970727 666448 6686681 60 61
CHR12.19970730 1078172 7364277 60 61
ref|NC_001145| 924430 8460517 60 61
ref|NC_001146| 784328 9400452 60 61
CHR15.19970811 1091283 10197892 60 61
ref|NC_001148| 948061 11307461 60 61
ref|NC_001224| 85779 12271444 60 61
# Kyokai first 10 lines
gi|347729946|dbj|BABQ01000001.1| 5491 137 80 81
gi|347729945|dbj|BABQ01000002.1| 4704 5834 80 81
gi|347729944|dbj|BABQ01000003.1| 203046 10734 80 81
gi|347729943|dbj|BABQ01000004.1| 1036 216456 80 81
gi|347729942|dbj|BABQ01000005.1| 126714 217642 80 81
gi|347729941|dbj|BABQ01000006.1| 121900 346077 80 81
gi|347729940|dbj|BABQ01000007.1| 20385 469638 80 81
gi|347729939|dbj|BABQ01000008.1| 243458 490415 80 81
gi|347729938|dbj|BABQ01000009.1| 189876 737054 80 81
gi|347729937|dbj|BABQ01000010.1| 81973 929441 80 81
# goes on for many thousands of lines
What you will notice is that there are significant differences between the 2 assemblies. The first is almost all in chromosomes, the second is largely mapped to various chromosomes but exists in thousands of tiny contigs. There is no way to just map reads to them, call variants, and assume that you (or bcftools) can figure out what is going on.
This is possibly an extreme example, but even in the case of two genomes assembled to chromosome pseudomolecules, there is basically no chance that those two chromosomes are on the same coordinate system or have similar metadata. That means again that bcftools probably won't be able to figure it out, because a SNP at position 112929 in genome A will be homologous to position 110201 in genome B.
Liftover: a way to solve the mapping problem between 2 genomes.
You can create a mapping between the two genomes. This is 100% a solved problem in bioinformatics, and it is more or less the same as doing liftovers between different assembly versions of the same genome. I think that you can more or less just use the UCSC LiftOver tool to convert SNP coordinates, unless I am missing something. That how-to guide that I linked is going to give you much better instructions than I can give you, here is a more minimal version that might be simpler.
The only extra suggestion that I have is that the guide uses BLAT to align the two assemblies, and minimap2 is a much faster tool. You may need to fight the input formatting a little bit, but if you run into runtime issues you might make that change. This is why I suggested it in the comment, though I wasn't aware that the liftover docs were so complete.
- Align reads (individuals of species A to reference A, individuals of species B to reference B).
- call variants-->get VCF files using preferred workflow for each of A and B independently.
- compute liftover mapping (chain file, usually) from reference A to reference B.
- use computed liftover mapping to actually lift VCF over from one genome to another.
- align A and B references, call SNPs between them (possibly can be done from (3) output). Ensure that you are using as reference the liftover target genome so that the VCF is informative.
- For each SNP private to each species (output of 5), query the BAMs to ensure that no variation exists.