Finding SNP/Indel Markers of Interspecific Variation

Question

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:

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

---

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

---

Scratch Materials:

1. Species A reference genome

---

2. Species B reference genome

---

3. Illumina short reads for Species A

---

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

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

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

4. 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)
5. 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
6. Repeat reciprocally to get alleles exclusive to Species A

Does the logic for this new pipeline check out?

Answer 1

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.

Answer 2

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.