<p>Ok, with invaluable aid from @Maximillian Press, I've figured out how to find alleles which are variable between two species but fixed within one species.</p>
<p>Here is how you find SNP/Indel loci which are distinct between species A (hereafter sA) and species B (hereafter sB) but are fixed WITHIN sA:</p>
<pre><code># Maximillian Press's Contribution: Map sA reference to sB reference using minimap2, then convert to bam. 

minimap2 -a sB_ref.fasta sA_ref.fasta &gt; AtoB.sam 
samtools view -bh AtoB.sam &gt; AtoB.bam 


# Sort bam by coordinates, then use bcftools to generate a pileup and call all alleles of sA relative to sB

java -jar <span class="math-container">$EBROOTPICARD/picard.jar SortSam I=AtoB.bam O=sorted_AtoB.bam CREATE_INDEX=true SORT_ORDER=coordinate TMP_DIR=$</span>TMPDIR    
bcftools mpileup -O b -o AtoB.bcf -f sB_ref.fasta sorted_AtoB.bam
bcftools call -m -v -o AtoB.vcf AtoB.bcf 


 # Call all alleles of sA relative to sA. 
# This represents the SNP/Indel loci which are variable WITHIN sA and thus ought to be discarded from the final output.
# &quot;bam_listA.txt&quot; contains names of all the coordinate-sorted bam files for species A
bcftools mpileup -O b -o AtoA.bcf -f sA_ref.fasta -b bam_listA.txt
bcftools call -m -v -o AtoA.vcf AtoA.bcf 

</code></pre>
<p>I then use the &quot;sed&quot; function to replace chromosome names with their analogs (happy to give more details if anyone else needs to do this) and also removed the VCF headers prior to making a C++ program which extracts all the rows found in AtoB.vcf but not in AtoA.vcf. My final output was thus all SNP/Indel loci which vary between sA and sB but NOT within sA.
Happy to share the C++ program if anyone else needs to do this.</p>

CC BY-SA 4.0
2021-03-07T03:07:14.550Z
2021-03-07T03:07:14.550Z
annabelperry
<h2>Problem statement</h2>
<p>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 <em>shared</em> between the two genomes (there is variation in each species), and which ones are <em>private</em> to only one of the species (in the other genome there is no variant at that location).</p>
<p>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 <em>liftover</em> operation.</p>
<h2>Genomes can have different coordinate systems.</h2>
<p>Let's first think about what exactly you are looking at with 2 different reference genomes. For example, here are outputs of <code>samtools faidx</code> for 2 different S. cerevisiae strains, S288C and Kyokai-7:</p>
<pre><code># S288C
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
</code></pre>
<p>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 <em>mapped</em> 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.</p>
<p>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.</p>
<h2>Liftover: a way to solve the mapping problem between 2 genomes.</h2>
<p>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 <a href="https://genome.sph.umich.edu/wiki/LiftOver" rel="nofollow noreferrer">liftovers</a> between different assembly versions of the same genome. I think that you can more or less just <a href="http://genomewiki.ucsc.edu/index.php/LiftOver_Howto#:%7E:text=A%20liftOver%20file%20is%20a,of%20a%20genome%20to%20another." rel="nofollow noreferrer">use the UCSC LiftOver tool to convert SNP coordinates</a>, 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, <a href="http://genomewiki.ucsc.edu/index.php/Minimal_Steps_For_LiftOver" rel="nofollow noreferrer">here is a more minimal version</a> that might be simpler.</p>
<p>The only extra suggestion that I have is that the guide uses BLAT to align the two assemblies, and minimap2 is a <em>much</em> 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.</p>
<h2>Short pipeline:</h2>
<ol>
<li>Align reads (individuals of species A to reference A, individuals of species B to reference B).</li>
<li>call variants--&gt;get VCF files using preferred workflow for each of A and B independently.</li>
<li>compute liftover mapping (chain file, usually) from reference A to reference B.</li>
<li>use computed liftover mapping to actually lift VCF over from one genome to another.</li>
<li>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.</li>
<li>For each SNP private to each species (output of 5), query the BAMs to ensure that no variation exists.</li>
</ol>

CC BY-SA 4.0
2021-03-01T20:32:56.490Z
Thank you for your detailed response! I still am not sure at which point in the pipeline I perform variant calling. Do I use the liftover output in lieu of the Species A and Species B reference genomes in the SNP-calling steps to get the SNPs of each species? Here's an example of what I'm thinking of:
`bcftools mpileup -f liftover  -b bam_listA | bcftools call -o ACompatibleCoordinates.vcf
bcftools mpileup -f liftover -b bam_listB | bcftools call -o BCompatibleCoordinates.vcf`
I would then search the output VCF files for alleles which are unique to each species
annabelperry
CC BY-SA 4.0
2021-03-01T21:49:41.690Z
I have added a short suggested workflow that indicates what I saw as the best option. I am not sure that I totally understand your suggestion, it looks like you are passing a `liftover` to bcftools but I don't totally understand what that is supposed to look like? You could pass a FASTA which is the output of a liftover to bcftools I guess, or just use it as a reference in alignment in the first place. Those would likely also work, though I don't have a great sense for which way is best.
Maximilian Press
CC BY-SA 4.0
2021-03-01T22:31:41.547Z
Ok, thank you! I see a problem with the proposed workflow. In my understanding, to be called as an SNP/Indel, the locus in question must have different alleles in the _reference_ and _sam/bam_. I am trying to find all SNPs/indels which are _exclusive_ to Species A and all SNPs/indels which are _exclusive_ to Species B (i.e. completely _lack_ variation within species). If I call SNPs/Indels of A relative to A, won't the resulting VCF file _not_ include alleles which are 100% fixed in A? I am about to post an answer to my own question which solves this issue, though I would appreciate critique
annabelperry
CC BY-SA 4.0
2021-03-01T22:39:24.077Z
@annabelperry in this case, can you not align the two genomes to call such variants, and then confirm post hoc that they do not have variation in the BAM files? e.g. `minimap2 -a genomea.fa genomeb.fa | samtools view -bh | bcftools` etc., and then scan those against whichever VCF is the reference? that should be quite simple.
Maximilian Press
CC BY-SA 4.0
2021-03-05T02:47:19.907Z
Ah, thank you @Maximilian Press! I will try that and upvote your answer if it works. Else, I will be back here shortly haha
annabelperry
CC BY-SA 4.0
2021-03-01T18:08:33.317Z
2021-03-01T22:42:21.120Z
2021-03-01T22:42:21.120Z
Maximilian Press
<p><a href="https://i.sstatic.net/GQgx9.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/GQgx9.png" alt="enter image description here" /></a>
This is a revised version of my original question, edited to include the portions already answered by @Maximilian Press.</p>
<h1>Goal:</h1>
<h2>Find SNP and indel alleles which are <em>fixed between</em> but <em>lack polymorphism within</em> two closely related species (I'm looking for the BLUE region of the Venn diagram)</h2>
<h1>Desired Outputs:</h1>
<ol>
<li>All SNPs/Indels of Species A which are <em>exclusive to</em> and <em>lack variation within</em> Species A</li>
</ol>
<hr />
<ol start="2">
<li>All SNPs/Indels of Species B which are <em>exclusive to</em> and <em>lack variation within</em> Species B</li>
</ol>
<hr />
<h1>Scratch Materials:</h1>
<ol>
<li>Species A reference genome</li>
</ol>
<hr />
<ol start="2">
<li>Species B reference genome</li>
</ol>
<hr />
<ol start="3">
<li>Illumina short reads for Species A</li>
</ol>
<hr />
<ol start="4">
<li>Illumina short reads for Species B</li>
</ol>
<hr />
<p>My current plan is to...</p>
<ol>
<li>Map Species B to Species A:</li>
</ol>
<pre><code>minimap2 -ax asm5 referenceA.fasta referenceB.fasta &gt; BtoA.sam
</code></pre>
<ol start="2">
<li>Find all SNP/Indel alleles of Species B relative to Species A (Represents the <em>entire</em> circle in the above Venn diagram):</li>
</ol>
<pre><code>bcftools mpileup -f referenceA.fasta -s BtoA.sam | bcftools call -o B_compared_to_A.vcf
</code></pre>
<ol start="3">
<li>Find all SNPs/Indel of Species B reads relative to Species B reference (Represents the <em>red</em> region of the Venn diagram):</li>
</ol>
<pre><code>bcftools mpileup -f referenceB.fasta -b bam_listB | bcftools call -o B_compared_to_B.vcf
</code></pre>
<ol start="4">
<li>Perform a liftover on &quot;B_compared_to_B.vcf&quot; to correct for coordinate differences between Species A and Species B references (I haven't completely figured out how to do this - I want to make sure I have the general pipeline nailed down before I get into the weeds)</li>
<li>Write a program which extracts all loci found in &quot;B_compared_to_A.vcf&quot; but NOT in &quot;B_compared_to_B_liftedOver.vcf&quot; (This should yield the BLUE region of the Venn Diagram</li>
<li>Repeat reciprocally to get alleles exclusive to Species A</li>
</ol>
<p>Does the logic for this new pipeline check out?</p>

CC BY-SA 4.0
2021-01-31T20:21:16.590Z
That seems like a valid approach to me.
Liam McIntyre
CC BY-SA 4.0
2021-01-31T21:55:18.633Z
I might be missing something, but does this approach not require that the coordinate systems of `referenceA.fasta` and `referenceB.fasta` be the same? I would expect a liftover step to map variant coordinates between the two. For example, the BAM headers in `bam_listA` will not match the fasta index of `referenceB.fasta`. Maybe these tools don't care about that, I don't know that off the top of my head. but I worry about coordinate frameshifting.
Maximilian Press
CC BY-SA 4.0
2021-01-31T21:59:23.063Z
At the risk of  sounding dumb, is there a reason that you don't use reciprocal whole-genome alignment via e.g. minimap2 to identify the interspecific differences, which you can then check with reads in each species separately? i don't have a specific workflow fleshed out, but i would hope that the reference assemblies were good enough to get the (vast) majority of things that would be picked up by reciprocal shotgun read mapping.
Maximilian Press
CC BY-SA 4.0
2021-02-03T23:23:37.127Z
@MaximilianPress could you link me to an article describing coordinate frameshifting? I am relatively new to bioinformatics and haven't encountered this concept before. I will also give minimap2 a look, I haven't heard of this software before.
annabelperry
CC BY-SA 4.0
2021-02-04T14:47:47.407Z
not really an article-level concept (and others would use different term probably), it's just that the two genomes probably don't have the exact same coordinates (e.g. chr1 is 10.5 Mbp in one genome, and 10.6 Mbp in the other genome). If this is true, then when you call a variant in one genome, those coordinates will usually differ from the coordinates for _the homologous variant_ in the other genome. Stated differently: it is very unlikely that two species' genomes can be aligned without gaps. Wherever there is a gap, variant coordinates will shift.
Maximilian Press
CC BY-SA 4.0
2021-02-28T01:14:06.040Z
@MaximilianPress I apologize for dredging up a nearly month-old conversation, but I am still struggling with this problem. When you mentioned "reciprocal whole-genome mapping" with minimap2, what did you intend for me to do with the reads once I have mapped the Species A reference genome to the Species B reference genome? I understand you didn't have a pipeline worked out, but I'm trying to figure out the logic of the approach.
annabelperry
CC BY-SA 4.0
2021-03-01T17:39:49.397Z
I will post as a short answer, easier to format etc.
Maximilian Press
<vcf><genomics><variant-calling><software-recommendation>
Finding SNP/Indel Markers of Interspecific Variation
# Finding SNP/Indel Markers of Interspecific Variation

(Asked by: username_0 on 2021-01-30)

This is a revised version of my original question, edited to include the portions already answered by @Maximilian Press.
Goal:
Find SNP and indel alleles which are fixed between but lack polymorphism within two closely related species (I'm looking for the BLUE region of the Venn diagram)
Desired Outputs:

All SNPs/Indels of Species A which are exclusive to and lack variation within Species A



All SNPs/Indels of Species B which are exclusive to and lack variation within Species B


Scratch Materials:

Species A reference genome



Species B reference genome



Illumina short reads for Species A



Illumina short reads for Species B


My current plan is to...

Map Species B to Species A:

minimap2 -ax asm5 referenceA.fasta referenceB.fasta > BtoA.sam


Find all SNP/Indel alleles of Species B relative to Species A (Represents the entire circle in the above Venn diagram):

bcftools mpileup -f referenceA.fasta -s BtoA.sam | bcftools call -o B_compared_to_A.vcf


Find all SNPs/Indel of Species B reads relative to Species B reference (Represents the red region of the Venn diagram):

bcftools mpileup -f referenceB.fasta -b bam_listB | bcftools call -o B_compared_to_B.vcf


Perform a liftover on "B_compared_to_B.vcf" to correct for coordinate differences between Species A and Species B references (I haven't completely figured out how to do this - I want to make sure I have the general pipeline nailed down before I get into the weeds)
Write a program which extracts all loci found in "B_compared_to_A.vcf" but NOT in "B_compared_to_B_liftedOver.vcf" (This should yield the BLUE region of the Venn Diagram
Repeat reciprocally to get alleles exclusive to Species A

Does the logic for this new pipeline check out?

**Tags:** vcf, genomics, variant-calling, software-recommendation
> **Comments:**
> **username_1**: That seems like a valid approach to me.  
> **username_2**: I might be missing something, but does this approach not require that the coordinate systems of `referenceA.fasta` and `referenceB.fasta` be the same? I would expect a liftover step to map variant coordinates between the two. For example, the BAM headers in `bam_listA` will not match the fasta index of `referenceB.fasta`. Maybe these tools don't care about that, I don't know that off the top of my head. but I worry about coordinate frameshifting.  
> **username_2**: At the risk of  sounding dumb, is there a reason that you don't use reciprocal whole-genome alignment via e.g. minimap2 to identify the interspecific differences, which you can then check with reads in each species separately? i don't have a specific workflow fleshed out, but i would hope that the reference assemblies were good enough to get the (vast) majority of things that would be picked up by reciprocal shotgun read mapping.  
> **username_0**: @MaximilianPress could you link me to an article describing coordinate frameshifting? I am relatively new to bioinformatics and haven't encountered this concept before. I will also give minimap2 a look, I haven't heard of this software before.  
> **username_2**: not really an article-level concept (and others would use different term probably), it's just that the two genomes probably don't have the exact same coordinates (e.g. chr1 is 10.5 Mbp in one genome, and 10.6 Mbp in the other genome). If this is true, then when you call a variant in one genome, those coordinates will usually differ from the coordinates for _the homologous variant_ in the other genome. Stated differently: it is very unlikely that two species' genomes can be aligned without gaps. Wherever there is a gap, variant coordinates will shift.  
> **username_0**: @MaximilianPress I apologize for dredging up a nearly month-old conversation, but I am still struggling with this problem. When you mentioned "reciprocal whole-genome mapping" with minimap2, what did you intend for me to do with the reads once I have mapped the Species A reference genome to the Species B reference genome? I understand you didn't have a pipeline worked out, but I'm trying to figure out the logic of the approach.  
> **username_2**: I will post as a short answer, easier to format etc.  

---

## 2 answers

### Answer by username_0 on 2021-03-07. Score: 1 | Accepted answer

Ok, with invaluable aid from @Maximillian Press, I've figured out how to find alleles which are variable between two species but fixed within one species.
Here is how you find SNP/Indel loci which are distinct between species A (hereafter sA) and species B (hereafter sB) but are fixed WITHIN sA:
# Maximillian Press's Contribution: Map sA reference to sB reference using minimap2, then convert to bam. 

minimap2 -a sB_ref.fasta sA_ref.fasta > AtoB.sam 
samtools view -bh AtoB.sam > AtoB.bam 


# Sort bam by coordinates, then use bcftools to generate a pileup and call all alleles of sA relative to sB

java -jar $EBROOTPICARD/picard.jar SortSam I=AtoB.bam O=sorted_AtoB.bam CREATE_INDEX=true SORT_ORDER=coordinate TMP_DIR=$TMPDIR    
bcftools mpileup -O b -o AtoB.bcf -f sB_ref.fasta sorted_AtoB.bam
bcftools call -m -v -o AtoB.vcf AtoB.bcf 


 # Call all alleles of sA relative to sA. 
# This represents the SNP/Indel loci which are variable WITHIN sA and thus ought to be discarded from the final output.
# "bam_listA.txt" contains names of all the coordinate-sorted bam files for species A
bcftools mpileup -O b -o AtoA.bcf -f sA_ref.fasta -b bam_listA.txt
bcftools call -m -v -o AtoA.vcf AtoA.bcf 


I then use the "sed" function to replace chromosome names with their analogs (happy to give more details if anyone else needs to do this) and also removed the VCF headers prior to making a C++ program which extracts all the rows found in AtoB.vcf but not in AtoA.vcf. My final output was thus all SNP/Indel loci which vary between sA and sB but NOT within sA.
Happy to share the C++ program if anyone else needs to do this.


---

### Answer by username_2 on 2021-03-01. Score: 2

Problem statement
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:
# S288C
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.
Short pipeline:

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.

> **Comments:**
> **username_0**: Thank you for your detailed response! I still am not sure at which point in the pipeline I perform variant calling. Do I use the liftover output in lieu of the Species A and Species B reference genomes in the SNP-calling steps to get the SNPs of each species? Here's an example of what I'm thinking of: `bcftools mpileup -f liftover  -b bam_listA | bcftools call -o ACompatibleCoordinates.vcf bcftools mpileup -f liftover -b bam_listB | bcftools call -o BCompatibleCoordinates.vcf` I would then search the output VCF files for alleles which are unique to each species  
> **username_2**: I have added a short suggested workflow that indicates what I saw as the best option. I am not sure that I totally understand your suggestion, it looks like you are passing a `liftover` to bcftools but I don't totally understand what that is supposed to look like? You could pass a FASTA which is the output of a liftover to bcftools I guess, or just use it as a reference in alignment in the first place. Those would likely also work, though I don't have a great sense for which way is best.  
> **username_0**: Ok, thank you! I see a problem with the proposed workflow. In my understanding, to be called as an SNP/Indel, the locus in question must have different alleles in the _reference_ and _sam/bam_. I am trying to find all SNPs/indels which are _exclusive_ to Species A and all SNPs/indels which are _exclusive_ to Species B (i.e. completely _lack_ variation within species). If I call SNPs/Indels of A relative to A, won't the resulting VCF file _not_ include alleles which are 100% fixed in A? I am about to post an answer to my own question which solves this issue, though I would appreciate critique  
> **username_2**: @annabelperry in this case, can you not align the two genomes to call such variants, and then confirm post hoc that they do not have variation in the BAM files? e.g. `minimap2 -a genomea.fa genomeb.fa | samtools view -bh | bcftools` etc., and then scan those against whichever VCF is the reference? that should be quite simple.  
> **username_0**: Ah, thank you @Maximilian Press! I will try that and upvote your answer if it works. Else, I will be back here shortly haha  

---

