Merge pull request #1 from gulki/patch-1

Patch 1

Author: CompEvoMetu

slim_simulation/run_slim.sh

@@ -0,0 +1,13 @@
+#!/bin/bash -l
+#SBATCH -p chimp
+#SBATCH -n 3
+#SBATCH -t 5-00:00:00
+#SBATCH -J slim
+#SBATCH -o slurm-%j-%N-%u.out
+#SBATCH -e slurm-%J-%N-%u.err
+
+simfile=$1
+outfile=$2
+
+#seed=$RANDOM
+/usr/local/sw/SLiM-3.7/build/slim -d seed=$RANDOM ${simfile} > ${outfile}

slim_simulation/ws_model1.txt

@@ -0,0 +1,39 @@
+// Simulations by SLiM - aDNAworkshop2022 - Happy simulations!
+
+// Don't worry, be happy!
+
+// Let's start:
+
+//manual page 86 for more information about initialize
+
+initialize()
+
+{
+        // set the overall mutation rate
+        initializeMutationRate(1e-7);
+
+        // m1 mutation type: neutral
+        initializeMutationType("m1", 0.5, "f", 0.0);
+
+        // g1 genomic element type: uses m1 for all mutations
+        initializeGenomicElementType("g1", m1, 1.0);
+
+        // uniform chromosome of length 100 kb
+        initializeGenomicElement(g1, 0, 99999);
+
+        // uniform recombination along the chromosome
+        initializeRecombinationRate(1e-8);
+
+}
+
+// create a population of 500 individuals
+1
+{
+    sim.addSubpop("p1", 500);
+}
+
+// run to generation 10000
+10000
+{
+    sim.simulationFinished();
+}

slim_simulation/ws_model2.txt

@@ -0,0 +1,31 @@
+// Simulations by SLiM - aDNAworkshop2022 - Happy simulations!
+
+// Let's add output:
+
+//manual page 93 for more information about output and also section 26.1.1 page 657
+
+initialize()
+
+{
+        // set the overall mutation rate
+        initializeMutationRate(1e-7);
+
+        // m1 mutation type: neutral
+        initializeMutationType("m1", 0.5, "f", 0.0);
+
+        // g1 genomic element type: uses m1 for all mutations
+        initializeGenomicElementType("g1", m1, 1.0);
+
+        // uniform chromosome of length 100 kb
+        initializeGenomicElement(g1, 0, 99999);
+
+        // uniform recombination along the chromosome
+        initializeRecombinationRate(1e-8);
+
+}
+
+// create a population of 500 individuals
+1 {sim.addSubpop("p1", 500);}
+
+// output:
+late() { sim.outputFull(); }

slim_simulation/ws_model3.txt

@@ -0,0 +1,47 @@
+// Simulations by SLiM - aDNAworkshop2022 - Happy simulations!
+
+// Let's run a nucleotide-based simulation and a population split scenario
+
+
+
+initialize() {
+        setSeed(seed);
+
+//simulate 100  Kb segment
+        defineConstant("L", 100000); 
+
+//set the mutation rate
+        defineConstant("mu", 1e-8);
+
+// nucleotide-based simulation
+        initializeSLiMOptions(nucleotideBased=T);
+
+//generate random order of ACGTs of length "L" defined above 
+        initializeAncestralNucleotides(randomNucleotides(L)); 
+
+//Introduce nucleotide-based neutral mutations -> dominance coeff is 0.5, DFE was set to fixed fittness effect, selection coeff is zero -> neutral
+        initializeMutationTypeNuc("m1", 0.5, "f", 0.0);
+        m1.convertToSubstitution = F;
+
+        initializeGenomicElementType("g1", m1, 1.0, mmJukesCantor(mu));
+        initializeGenomicElement(g1, 0, L-1);
+
+        //set the recombination rate
+        initializeRecombinationRate(1e-8);
+}
+
+
+// create a population of 10000 individuals
+1 {sim.addSubpop("p1", 10000);}
+
+
+// Population split
+
+1000 { sim.addSubpopSplit("p2", 2000, p1); }
+
+
+// output:
+
+1000 late() {
+        p2.outputVCFSample(5, replace=T, simplifyNucleotides=T);
+}

slim_simulation/ws_model4.txt

@@ -0,0 +1,47 @@
+// Simulations by SLiM - aDNAworkshop2022 - Happy simulations!
+
+// Mutationtype - page 144
+
+
+
+initialize() {
+        setSeed(seed);
+
+//simulate 100  Kb segment
+        defineConstant("L", 100000); 
+
+//set the mutation rate
+        defineConstant("mu", 1e-8);
+
+// nucleotide-based simulation
+        initializeSLiMOptions(nucleotideBased=T);
+
+//generate random order of ACGTs of length "L" defined above 
+        initializeAncestralNucleotides(randomNucleotides(L)); 
+
+//Introduce nucleotide-based neutral mutations -> dominance coeff is 0.5, DFE was set to fixed fittness effect, selection coeff is zero -> neutral
+        initializeMutationTypeNuc("m1", 0.1, "g", -0.03, 0.2);
+        m1.convertToSubstitution = F;
+
+        initializeGenomicElementType("g1", m1, 1.0, mmJukesCantor(mu));
+        initializeGenomicElement(g1, 0, L-1);
+
+        //set the recombination rate
+        initializeRecombinationRate(1e-8);
+}
+
+
+// create a population of 10000 individuals
+1 {sim.addSubpop("p1", 10000);}
+
+
+// Population split
+
+1000 { sim.addSubpopSplit("p2", 2000, p1); }
+
+
+// output:
+
+1000 late() {
+        p2.outputVCFSample(5, replace=T, simplifyNucleotides=T);
+}

slim_simulation/ws_model5.txt

@@ -0,0 +1,119 @@
+// A more complex scenario
+
+// Example paper 1: https://www.nature.com/articles/s41467-018-04191-y 10mb segment
+// Example paper 2: https://genome.cshlp.org/content/24/6/885.long 4 mb segment
+// Multi-population model of ancient Eurasia -> https://stdpopsim.readthedocs.io/en/latest/catalog.html#sec_catalog_HomSap 
+// https://www.tandfonline.com/doi/full/10.1080/01621459.2019.1635482
+
+
+
+// I use mutation rate as 1e-8
+
+// Model-2: Neutral simulation - add mutation to the middle of the sequence
+
+initialize() {
+	setSeed(seed);
+
+//simulate 100 kb segment
+	defineConstant("L", 100000); 
+
+//set the mutation rate
+	defineConstant("mu", 1e-8);
+
+// nucleotide-based simulation
+	initializeSLiMOptions(nucleotideBased=T);
+
+//generate random order of ACGTs of length "L" defined above 
+	initializeAncestralNucleotides(randomNucleotides(L)); 
+
+//Introduce nucleotide-based neutral mutations -> dominance coeff is 0.5, DFE was set to fixed fittness effect, selection coeff is zero -> neutral
+	initializeMutationTypeNuc("m1", 0.5, "f", 0.0);
+	m1.convertToSubstitution = F;
+
+//Non-nucloide based (target) neutral mutation to be introduced at certain time point.
+ 	initializeMutationType("m2", 0.5, "f", 0.0);
+ 	m2.convertToSubstitution = F;
+
+//simulate genomic segment of length of L
+
+	initializeGenomicElementType("g1", m1, 1.0, mmJukesCantor(mu));
+	initializeGenomicElement(g1, 0, L-1);
+
+
+//set the recombination rate
+
+	initializeRecombinationRate(1e-8);
+}
+
+// initial population size 29,100 individuals ancestral population (assume Mbuty and WHG) -> simulate for 58,000 generation -> 1,450,000 years 
+
+	1 { 
+	defineConstant("simID", getSeed());
+	sim.addSubpop("p1", 29100);}
+
+// Time of WHG and Mbuti split at 95,800 years ago (3,832 generations) "p2 is Mbuty". Mbuty pop. size is 17,300
+
+	54164 { sim.addSubpopSplit("p2", 17300, p1); }
+
+// Time of WHG and Basal Eurasian split at 79,800 years ago (3,192 generations), "p3 is Basal Eurasian". Basal Eurasian pop. size is 1,920
+
+	54808 { sim.addSubpopSplit("p3", 1920, p1); }
+
+// Time of WHG and Ust’Ishim split at 51,500 years ago (2,060 generations), "p4 is Ust Ishim". Ust Ishim pop size is 1,920
+
+	55940 { sim.addSubpopSplit("p4", 1920, p1); }
+
+// Time of WHG and Han Chinese split at 50,400 years ago (2,016 generations), "p5 is Han Chinese". Han. Chinese pop size is 6,300
+
+	55984 { sim.addSubpopSplit("p5", 6300, p1); }
+
+// Time of WHG and Mal’ta split at 44,900 years ago (1,796 generation), "p6 is Mal'ta". Mal'ta pop size is 1,920
+
+	56204 { sim.addSubpopSplit("p6", 1920, p1); }
+
+// Time of WHG and AnatoliaNeolithic split at 15,000 years ago (600 generations), "p7 is Anatolia Neolithic". Anatolia Neolithic pop size is 2000
+
+	57400 { sim.addSubpopSplit("p7", 2000, p1); }
+
+
+// Add the target neutral mutation (m2) at generation 57500 to 100 individuals
+// Change 57500 to any time you wish
+57500 late() {
+	// Save the current state of simulation into a tmp file.
+	sim.outputFull("/tmp/slim_" + simID + ".txt");
+
+	// Select 100 random haploid indvidual from p7
+	targetinds = sample(p7.genomes, 100);
+
+	// Introduce m2 mutation into the center of the sequence of these 100 individuals
+	targetinds.addNewDrawnMutation(m2, asInteger(L/2));
+}
+
+
+// Check if the target mutation (m2) is lost
+// Change 5750 to any time you wish
+57500:58000 late() {
+	if (sim.countOfMutationsOfType(m2) == 0) {
+ 		print(simID + ": LOST-RESTARTING ");
+
+		 // if the mutation is lost, restart simulation from the tmp file
+ 		sim.readFromPopulationFile("/tmp/slim_" + simID + ".txt");
+ 		setSeed(getSeed() + 1);
+ 		target = sample(p7.genomes, 100);
+ 		target.addNewDrawnMutation(m2, asInteger(L/2));
+}	}
+
+
+
+57600 late() {
+	p7.outputVCFSample(1000, replace=T, simplifyNucleotides=T);
+}
+
+57800 late() {
+	p7.outputVCFSample(1000, replace=T, simplifyNucleotides=T);
+}
+
+58000 late() {
+	p7.outputVCFSample(1000, replace=T, simplifyNucleotides=T);
+	deleteFile("/tmp/slim_" + simID + ".txt");
+}