Skip to content

Commit 4249354

Browse files
bhallerbhaller
committed
rewrite recipe 13.3 to use multiple chromosomes
1 parent 0ddc806 commit 4249354

File tree

2 files changed

+46
-74
lines changed

2 files changed

+46
-74
lines changed

‎QtSLiM/recipes/Recipe 13.3 - A model of discrete QTL effects with a structured chromosome.txt‎

Lines changed: 23 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,6 @@
11
// Keywords: migration, dispersal, QTL, quantitative trait loci
22

33
initialize() {
4-
initializeMutationRate(1e-6);
5-
64
// neutral mutations in non-coding regions
75
initializeMutationType("m1", 0.5, "f", 0.0);
86
initializeGenomicElementType("g1", m1, 1.0);
@@ -15,30 +13,20 @@ initialize() {
1513
m2.mutationStackPolicy = "l";
1614

1715
// set up our chromosome: 10 QTLs, surrounded by neutral regions
18-
defineConstant("C", 10); // number of QTLs
16+
defineConstant("C", 10); // number of QTLs / chromosomes
1917
defineConstant("W", 1000); // size of neutral buffer on each side
20-
pos = 0;
21-
q = NULL;
2218

2319
for (i in 1:C)
2420
{
25-
initializeGenomicElement(g1, pos, pos + W-1);
26-
pos = pos + W;
21+
initializeChromosome(i, 0, W*2 + 1);
2722

28-
initializeGenomicElement(g2, pos, pos);
29-
q = c(q, pos);
30-
pos = pos + 1;
23+
initializeGenomicElement(g1, 0, W-1);
24+
initializeGenomicElement(g2, W, W);
25+
initializeGenomicElement(g1, W+1, W+1 + W-1);
3126

32-
initializeGenomicElement(g1, pos, pos + W-1);
33-
pos = pos + W;
27+
initializeMutationRate(1e-6);
28+
initializeRecombinationRate(1e-8);
3429
}
35-
36-
defineConstant("Q", q); // remember our QTL positions
37-
38-
// we want the QTLs to be unlinked; build a recombination map for that
39-
rates = c(rep(c(1e-8, 0.5), C-1), 1e-8);
40-
ends = (repEach(Q + W, 2) + rep(c(0,1), C))[0:(C*2 - 2)];
41-
initializeRecombinationRate(rates, ends);
4230
}
4331
1 late() {
4432
sim.addSubpop("p1", 500);
@@ -48,17 +36,15 @@ initialize() {
4836
p1.setMigrationRates(p2, 0.01);
4937
p2.setMigrationRates(p1, 0.01);
5038

51-
// replicate the s1 event to output in tick 2 also
52-
community.registerEarlyEvent("s2", s1.source, 2, 2);
53-
5439
// optional: give m2 mutations to everyone, as standing variation
55-
g = sim.subpopulations.haplosomes;
40+
individuals = sim.subpopulations.individuals;
5641

57-
for (q in Q)
42+
for (i in 1:C)
5843
{
44+
g = sim.subpopulations.individuals.haplosomesForChromosomes(i);
5945
isPlus = asLogical(rbinom(size(g), 1, 0.5));
60-
g[isPlus].addNewMutation(m2, 1.0, q);
61-
g[!isPlus].addNewMutation(m2, -1.0, q);
46+
g[isPlus].addNewMutation(m2, 1.0, W);
47+
g[!isPlus].addNewMutation(m2, -1.0, W);
6248
}
6349
}
6450
mutationEffect(m2) { return 1.0; }
@@ -78,7 +64,7 @@ mateChoice() {
7864
others = sourceSubpop.individuals.tagF;
7965
return weights * dnorm(others, phenotype, 5.0);
8066
}
81-
s1 2001 early() {
67+
c(2,2001) early() {
8268
cat("-------------------------------\n");
8369
cat("Output for end of cycle " + (sim.cycle - 1) + ":\n\n");
8470

@@ -96,18 +82,18 @@ s1 2001 early() {
9682
minus = muts[muts.selectionCoeff == -1.0];
9783

9884
cat("\nOverall frequencies:\n\n");
99-
for (q in Q)
85+
for (i in 1:C)
10086
{
101-
qPlus = plus[plus.position == q];
102-
qMinus = minus[minus.position == q];
103-
pf = sum(sim.mutationFrequencies(NULL, qPlus));
104-
mf = sum(sim.mutationFrequencies(NULL, qMinus));
105-
pf1 = sum(sim.mutationFrequencies(p1, qPlus));
106-
mf1 = sum(sim.mutationFrequencies(p1, qMinus));
107-
pf2 = sum(sim.mutationFrequencies(p2, qPlus));
108-
mf2 = sum(sim.mutationFrequencies(p2, qMinus));
87+
iPlus = plus[plus.chromosome.id == i];
88+
iMinus = minus[minus.chromosome.id == i];
89+
pf = sum(sim.mutationFrequencies(NULL, iPlus));
90+
mf = sum(sim.mutationFrequencies(NULL, iMinus));
91+
pf1 = sum(sim.mutationFrequencies(p1, iPlus));
92+
mf1 = sum(sim.mutationFrequencies(p1, iMinus));
93+
pf2 = sum(sim.mutationFrequencies(p2, iPlus));
94+
mf2 = sum(sim.mutationFrequencies(p2, iMinus));
10995

110-
cat(" QTL " + q + ": f(+) == " + pf + ", f(-) == " + mf + "\n");
96+
cat(" QTL " + i + ": f(+) == " + pf + ", f(-) == " + mf + "\n");
11197
cat(" in p1: f(+) == " + pf1 + ", f(-) == " + mf1 + "\n");
11298
cat(" in p2: f(+) == " + pf2 + ", f(-) == " + mf2 + "\n\n");
11399
}

‎SLiMgui/Recipes/Recipe 13.3 - A model of discrete QTL effects with a structured chromosome.txt‎

Lines changed: 23 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,6 @@
11
// Keywords: migration, dispersal, QTL, quantitative trait loci
22

33
initialize() {
4-
initializeMutationRate(1e-6);
5-
64
// neutral mutations in non-coding regions
75
initializeMutationType("m1", 0.5, "f", 0.0);
86
initializeGenomicElementType("g1", m1, 1.0);
@@ -15,30 +13,20 @@ initialize() {
1513
m2.mutationStackPolicy = "l";
1614

1715
// set up our chromosome: 10 QTLs, surrounded by neutral regions
18-
defineConstant("C", 10); // number of QTLs
16+
defineConstant("C", 10); // number of QTLs / chromosomes
1917
defineConstant("W", 1000); // size of neutral buffer on each side
20-
pos = 0;
21-
q = NULL;
2218

2319
for (i in 1:C)
2420
{
25-
initializeGenomicElement(g1, pos, pos + W-1);
26-
pos = pos + W;
21+
initializeChromosome(i, 0, W*2 + 1);
2722

28-
initializeGenomicElement(g2, pos, pos);
29-
q = c(q, pos);
30-
pos = pos + 1;
23+
initializeGenomicElement(g1, 0, W-1);
24+
initializeGenomicElement(g2, W, W);
25+
initializeGenomicElement(g1, W+1, W+1 + W-1);
3126

32-
initializeGenomicElement(g1, pos, pos + W-1);
33-
pos = pos + W;
27+
initializeMutationRate(1e-6);
28+
initializeRecombinationRate(1e-8);
3429
}
35-
36-
defineConstant("Q", q); // remember our QTL positions
37-
38-
// we want the QTLs to be unlinked; build a recombination map for that
39-
rates = c(rep(c(1e-8, 0.5), C-1), 1e-8);
40-
ends = (repEach(Q + W, 2) + rep(c(0,1), C))[0:(C*2 - 2)];
41-
initializeRecombinationRate(rates, ends);
4230
}
4331
1 late() {
4432
sim.addSubpop("p1", 500);
@@ -48,17 +36,15 @@ initialize() {
4836
p1.setMigrationRates(p2, 0.01);
4937
p2.setMigrationRates(p1, 0.01);
5038

51-
// replicate the s1 event to output in tick 2 also
52-
community.registerEarlyEvent("s2", s1.source, 2, 2);
53-
5439
// optional: give m2 mutations to everyone, as standing variation
55-
g = sim.subpopulations.haplosomes;
40+
individuals = sim.subpopulations.individuals;
5641

57-
for (q in Q)
42+
for (i in 1:C)
5843
{
44+
g = sim.subpopulations.individuals.haplosomesForChromosomes(i);
5945
isPlus = asLogical(rbinom(size(g), 1, 0.5));
60-
g[isPlus].addNewMutation(m2, 1.0, q);
61-
g[!isPlus].addNewMutation(m2, -1.0, q);
46+
g[isPlus].addNewMutation(m2, 1.0, W);
47+
g[!isPlus].addNewMutation(m2, -1.0, W);
6248
}
6349
}
6450
mutationEffect(m2) { return 1.0; }
@@ -78,7 +64,7 @@ mateChoice() {
7864
others = sourceSubpop.individuals.tagF;
7965
return weights * dnorm(others, phenotype, 5.0);
8066
}
81-
s1 2001 early() {
67+
c(2,2001) early() {
8268
cat("-------------------------------\n");
8369
cat("Output for end of cycle " + (sim.cycle - 1) + ":\n\n");
8470

@@ -96,18 +82,18 @@ s1 2001 early() {
9682
minus = muts[muts.selectionCoeff == -1.0];
9783

9884
cat("\nOverall frequencies:\n\n");
99-
for (q in Q)
85+
for (i in 1:C)
10086
{
101-
qPlus = plus[plus.position == q];
102-
qMinus = minus[minus.position == q];
103-
pf = sum(sim.mutationFrequencies(NULL, qPlus));
104-
mf = sum(sim.mutationFrequencies(NULL, qMinus));
105-
pf1 = sum(sim.mutationFrequencies(p1, qPlus));
106-
mf1 = sum(sim.mutationFrequencies(p1, qMinus));
107-
pf2 = sum(sim.mutationFrequencies(p2, qPlus));
108-
mf2 = sum(sim.mutationFrequencies(p2, qMinus));
87+
iPlus = plus[plus.chromosome.id == i];
88+
iMinus = minus[minus.chromosome.id == i];
89+
pf = sum(sim.mutationFrequencies(NULL, iPlus));
90+
mf = sum(sim.mutationFrequencies(NULL, iMinus));
91+
pf1 = sum(sim.mutationFrequencies(p1, iPlus));
92+
mf1 = sum(sim.mutationFrequencies(p1, iMinus));
93+
pf2 = sum(sim.mutationFrequencies(p2, iPlus));
94+
mf2 = sum(sim.mutationFrequencies(p2, iMinus));
10995

110-
cat(" QTL " + q + ": f(+) == " + pf + ", f(-) == " + mf + "\n");
96+
cat(" QTL " + i + ": f(+) == " + pf + ", f(-) == " + mf + "\n");
11197
cat(" in p1: f(+) == " + pf1 + ", f(-) == " + mf1 + "\n");
11298
cat(" in p2: f(+) == " + pf2 + ", f(-) == " + mf2 + "\n\n");
11399
}

0 commit comments

Comments
 (0)