Walkthrough

Introduction

This Walkthrough demonstrates how MimcrEE can be used to model experimental evolution. We will test the performance with a ROC graph for strongly selected SNPs (s=0.1) and weakly selected SNPs (s=0.025). The whole analysis will be performed in replicates of 10 in order to obtain error bars for the resulting ROC graphs.

Requirments

• data for the Walkthrough
• Software for the Walktrough
• MimicrEESummary
• PoPoolation2
• ROCR

Walkthrough

randomly pick SNPs

We first randomly pick 10 SNPs with the selection coefficients s=0.1 and s=0.025; Documentation for the scripts can be displayed with the --help option. The file Dmel-example.mimhap contains a small example of a natural population of D. melanogaster. In more detail the file contains 8000 SNPs from the four major chromosomes '2L, 2R, 3L, 3R' from a population of size N=250 (a small population was used in order to speed up the example).

gzip -dc Dmel-example.mimhap.gz>  Dmel-example.mimhap
mkdir snps
for i in {1..10}; do python <path>/pick-random-addtive-snps.py -n 10 -s 0.1 -e 0.5 -m 0.8 --loci-count 4000 --input Dmel-example.mimhap > snps/s01-n$i.txt; done
for i in {1..10}; do python <path>/pick-random-addtive-snps.py -n 10 -s 0.025 -e 0.5 -m 0.8 --loci-count 4000 --input Dmel-example.mimhap > snps/s0025-n$i.txt; done
rm Dmel-example.mimhap

Simulations with MimicrEE

We perform forward simulations of EE with MimicrEE using 60 generations and 3 replicates.

mkdir sim
for i in {1..10}; do nohup java -Xmx1g -jar  <path>/MimicrEESummary.jar --haplotypes-g0 Dmel-example.mimhap.gz --recombination-rate dmel.rr.txt --output-mode 60 --replicate-runs 3 --output-format sync --threads 6 --additive snps/s01-n$i.txt --output-file sim/s01-n$i.sync; done
for i in {1..10}; do nohup java -Xmx1g -jar <path>/MimicrEESummary.jar --haplotypes-g0 Dmel-example.mimhap.gz --recombination-rate dmel.rr.txt --output-mode 60 --replicate-runs 3 --output-format sync --threads 6 --additive snps/s0025-n$i.txt --output-file sim/s0025-n$i.sync; done

Run the cmh-test

mkdir cmh
for i in {1..10}; do perl <path>/popoolation2/cmh-test.pl --min-count 1 --min-coverage 1 --max-coverage 100000 --min-logpvalue 0.0 --population 1-2,3-4,5-6 --input sim/s01-n$i.sync --output cmh/s01-n$i.cmh; done &
for i in {1..10}; do perl <path>/popoolation2/cmh-test.pl --min-count 1 --min-coverage 1 --max-coverage 100000 --min-logpvalue 0.0 --population 1-2,3-4,5-6 --input sim/s0025-n$i.sync --output cmh/s0025-n$i.cmh; done &

Create labels and predictions for ROCR

mkdir simres
python <path>/rocr-generate-labellist.py sim/s01-n1.sync snps/s01-n{1..10}.txt > simres/s01.labels
python <path>/rocr-generate-labellist.py sim/s01-n1.sync snps/s0025-n{1..10}.txt > simres/s0025.labels
python <path>/rocr-generate-predictionlist.py cmh/s01-n{1..10}.cmh > simres/s01.predictions
python <path>/rocr-generate-predictionlist.py cmh/s0025-n{1..10}.cmh > simres/s0025.predictions

Visualize the results with ROCR

R --vanilla --args simres/s01.labels simres/s01.predictions simres/s0025.labels simres/s0025.predictions roc < ../toysoft/toyroc.R

The result

https://docs.google.com/file/d/0BweZxn-dqeZDNWY4V1pXaTVpeTA/edit?usp=sharing