# Finding SNP/Indel Markers of Interspecific Variation

This is a revised version of my original question, edited to include the portions already answered by @Maximilian Press.

# Desired Outputs:

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

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

# Scratch Materials:

1. Species A reference genome

1. Species B reference genome

1. Illumina short reads for Species A

1. Illumina short reads for Species B

My current plan is to...

1. Map Species B to Species A:
minimap2 -ax asm5 referenceA.fasta referenceB.fasta > BtoA.sam

1. 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

1. 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

1. 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)
2. 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
3. Repeat reciprocally to get alleles exclusive to Species A

Does the logic for this new pipeline check out?

• That seems like a valid approach to me. – Liam McIntyre Jan 31 at 20:21
• 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 Jan 31 at 21:55
• 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 Jan 31 at 21:59
• @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 Feb 3 at 23:23
• 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 Feb 4 at 14:47

## 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:

1. Align reads (individuals of species A to reference A, individuals of species B to reference B).
2. call variants-->get VCF files using preferred workflow for each of A and B independently.
3. compute liftover mapping (chain file, usually) from reference A to reference B.
4. use computed liftover mapping to actually lift VCF over from one genome to another.
5. 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.
6. For each SNP private to each species (output of 5), query the BAMs to ensure that no variation exists.
• 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 Mar 1 at 20:32
• 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 Mar 1 at 21:49
• 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 Mar 1 at 22:31
• @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 Mar 1 at 22:39
• 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 Mar 5 at 2:47

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.