Using population-level sequencing data in SLiM for neutral evolution simulations

[B-Demirbas]

Hi! Following a question I posted on slim-discuss (which was answered), Ben Haller suggested I raise an issue on the GitHub on the specific problem I'm facing.

My initial goal for a project is to model neutral evolution in a population of ~500 individuals over N generations using standard Wright-Fisher simulations. The project is supposed to complement experimental work, and hence I want to use real sequencing data in the form of .vcf files containing information on SNPs.

The problem is that I only have population sequence data (for a population of ~500 individuals). It does contain information on average allele frequencies, but there is no individual-level data. As far as I can tell from the SLiM manual there is no straightforward way to load in and use population-level .vcf files directly in SLiM (it seems it always assumes it is individual-level data), is this correct?

As a workaround I now use the allelic depth data to estimate the frequency of individuals with a given SNP in a population. I use that to generate artificial individuals in a new .vcf (in Python) which I then load into SLiM. This seems to work, but it would be more convenient if there was a way to do this directly in SLiM. I attached an (abridged) example of a .vcf file that I would like to use.

example VCF population single chromosome.txt

[bhaller]

Hi @B-Demirbas! Thanks for filing this, it's very helpful. I'm looking at the VCF now. It's interesting; I've never seen a VCF without individual-level calls before! :-O Here is a call line (the first one):

NC_003279.8 620 . C A 100.008 PASS DP=37;VDB=0.000145619;SGB=-0.662043;RPBZ=-0.548832;MQBZ=-0.55636;MQSBZ=-0.0199689;BQBZ=0.169721;SCBZ=1.23205;MQ0F=0;AF1=0.5;AC1=1;DP4=18,6,5,4;MQ=49;FQ=103.016;PV4=0.399859,0.480028,0.236702,0.265776 GT:PL:AD 0/1:130,0,253:24,9

Do you know which of these fields is the frequency of the mutation in the population that I would want to use? The headers are rather unclear to me – I don't work with empirical data. I'm guessing maybe these are relevant:

##INFO=<ID=AF1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele frequency (assuming HWE)">
##INFO=<ID=AF2,Number=1,Type=Float,Description="Max-likelihood estimate of the first and second group ALT allele frequency (assuming HWE)">
##INFO=<ID=AC1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele count (no HWE assumption)">

but I'm not sure which one I'd want to use, etc. And are these fields standard at all, or are they just what your particular pipeline produces? It looks like for this file AF1 is always 0.5 or 1, and AC is correspondingly either 1 or 2, so it looks like maybe data that was originally obtained from a single diploid individual; am I understanding correctly?

I'm trying to wrap my brain around how SLiM might possibly support this. I guess a method that adds mutations to a set of individuals from a population VCF file? Like:

individuals.addMutationsFromVCF(string$ filePath, [string$ vcfFrequencyKey = "AF1"], [iso$ mutType = "m1"])

or something?

[bhaller]

Anyhow, don't hold your breath on this – I'm in the middle of a rather large project right now, so I won't have cycles for this for probably several months. But it's a good idea that I'd like to do; it would be useful for lots of people, I imagine. Thanks!

[B-Demirbas]

Hi Ben, I'm glad to hear that you find the problem interesting, at least!

Do you know which of these fields is the frequency of the mutation in the population that I would want to use? The headers are rather unclear to me – I don't work with empirical data. I'm guessing maybe these are relevant:

I use the allelic depth field (AD in GT:PL:AD 0/1:130,0,253:24,9 in the example) to calculate allelic frequency. So in the example, total_count = 24 + 9, and alt_count = 9. Then, allelic frequency = alt_count / total_count.

And are these fields standard at all, or are they just what your particular pipeline produces? It looks like for this file AF1 is always 0.5 or 1, and AC is correspondingly either 1 or 2, so it looks like maybe data that was originally obtained from a single diploid individual; am I understanding correctly?

I don’t have much experience working with this particular data format either, unfortunately. As far as I understand it’s fairly standard, but I have asked my experimental collaborator to chip in. She will try to reply to this thread tomorrow!

individuals.addMutationsFromVCF(string$ filePath, [string$ vcfFrequencyKey = "AF1"], [iso$ mutType = "m1"])

Yes, something like this could work I think, if it automatically calculates allelic frequency from e.g. the AD field.

Anyhow, don't hold your breath on this – I'm in the middle of a rather large project right now, so I won't have cycles for this for probably several months. But it's a good idea that I'd like to do; it would be useful for lots of people, I imagine. Thanks!

No worries, for now I think my workaround works fine. It would be cool if this could help other people integrate their data in SLiM models!