reading SNPs from a VCF file

[andrea mrnjavac]

Hi all,

I am trying to read SNPs from a VCF file to set up the initial population state in a simple neutral WF simulation, however, slim omits about half of the SNPs from the VCF file. I ran the simulation 100 times, and slim always imports the same set of SNPs from a VCF file. I have checked some of the omitted SNPs and did not find anything unusual about them.

Does anybody have an idea what might be the problem?

This is my code:

"

initialize() {

initializeMutationRate(0);

initializeMutationType("m1", 0.5, "f", 0.0);

m1.convertToSubstitution = F;

initializeGenomicElementType("g1", m1, 1.0);

initializeGenomicElement(g1, 0, 23513711); //chr 2L

//read recombination map for chr 2L

lines = readFile("~/coo6_2L.tbl");

rates = NULL;

ends = NULL;

for (line in lines)

{

components = strsplit(line, "\t");

ends = c(ends, asInteger(components[0]));

rates = c(rates, asFloat(components[1]));

}

ends = c(ends[1:(size(ends)-1)] - 2, 23513711);

rates = rates * 1e-8;

initializeRecombinationRate(rates, ends);

}

1 late() {

sim.addSubpop("p1", 8);

p1.genomes.readFromVCF("~/line10_2L.recode.vcf", m1);

}

1 late(){ sim.outputFull("~/test_slim_gen1.out"); }

6 late(){ sim.outputFull("~/test_slim_gen6.out"); }

"

Cheers,

Andrea

[Zuxi Cui]

Hi Andrea,

I did something similar recently. Have you checked the minor allele count for those SNPs omitted? You may have lost them if there's no minor.

Best,

Terry

--
SLiM forward genetic simulation: http://messerlab.org/slim/
---
You received this message because you are subscribed to the Google Groups "slim-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to slim-discuss...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/slim-discuss/921b37f0-fc29-44a2-bb0e-e3b63d5cb85fn%40googlegroups.com.

[Ben Haller]

Hi Andrea.  I'm not sure offhand what the issue might be.  Can you send me everything needed to reproduce the problem, off-list?  Please make the script runnable in SLiMgui (i.e., not needing any parameters supplied on the command line).  Also please simplify the script as much as possible; i.e., if the recombination rate map can be eliminated while still reproducing the bug, please make those sorts of simplifying changes.  Thanks!

Cheers,

-B.

Benjamin C. Haller

Messer Lab

Cornell University

[andrea mrnjavac]

Hi Terry and Ben,

Thank you for your input. After checking my data more thoroughly, I've realized that most of the missing SNPs can be explained with p(alternative allele)=0.

Also, Ben, thanks for making a 700 pg manual and debugging everybody's problems every day!

Best,

Andrea

[Ben Haller]

No worries; it's my job and I enjoy it.  :->

Let me know if it seems like there's any remaining issue here (you do say "most", not "all"!).

Cheers,

-B.

Benjamin C. Haller

Messer Lab

Cornell University