Example2-DaDi

Compares SFS generation using executable GOFish (GO_Fish and Spectrum) against δaδi.

See subsection δaδi-to-GOFish Parameter Conversion in Getting Started for how to translate parameters from δaδi to GO_Fish in general.

Demography of population 0 (AF) and population 1 (EU)

A complex demographic scenario was chosen as a test case to compare the GO_Fish simulation against an already established SFS method, δaδi. The demographic model is from the YRI-CEU (AF-EU) δaδi example. Using δaδi parameterization to describe the model, the ancestral population, in mutation-selection equilibrium, undergoes an immediate expansion from N_ref to 2N_ref individuals. After time T1 (= 0.005) the population splits into two with a constant, equivalent migration, m_EU-AF (= 1) between the now split populations. The second (EU) population undergoes a severe bottleneck of 0.05N_ref when splitting from the first (AF) population, followed by exponential growth over time T2 (= 0.045) to size 5N_ref. The mutations are weakly deleterious and co-dominant (2N_refs = -2, h =0.5), where 1001 samples were taken of the EU population. The spectrum was then normalized by the number of segregating sites. The corresponding GO Fish parameters for the evolutionary scenario, given a mutation rate of 1x10^-9 per site, 2x10⁹ sites, and an initial population size, N_ref, of 10,000, are: T1 = 0.005*2N_ref = 100 generations, T2 = 900 generations, m_EU-AF = 1/(2N_ref) = 0.00005, 2N_refs = -4, h =0.5, and F = 0. As in δaδi, the population size/time can be scaled together and the simulation will generate the same normalized spectra. Below is the GO_Fish simulation and Spectrum SFS code:

/*
 * run.cu
 *
 *      Author: David Lawrie
 */

#include "go_fish.cuh"
#include "spectrum.h"
#include <vector>

/*
 * population 0 starts of in mutation-selection equilibrium at size N_ind
 * population 0 grows to size 2*N_ind at generation 1
 * at generation 0.01*N_ind, population 1 splits off from population 0
 * population 1's initial size is 0.05*N_ind
 * population 1 grows exponentially to size 5*N_ind for 0.09*N_ind generations
 *
 * migration between populations is at rate 1/(2*N_ind) starting at generation 0.01*N_ind+1
 *
 * selection is weakly deleterious (gamma = -4), mutations are co-dominant (h = 0.5), populations are outbred (F = 0)
 */

void run_validation_test(){
    typedef Sim_Model::demography_constant dem_const;
    typedef Sim_Model::demography_population_specific<dem_const,dem_const> dem_pop_const;
    typedef Sim_Model::demography_piecewise<dem_pop_const,dem_pop_const> init_expansion;
    typedef Sim_Model::demography_exponential_growth exp_growth;
    typedef Sim_Model::demography_population_specific<dem_const,exp_growth> dem_pop_const_exp;

    typedef Sim_Model::migration_constant_equal mig_const;
    typedef Sim_Model::migration_constant_directional<mig_const> mig_dir;
    typedef Sim_Model::migration_constant_directional<mig_dir> mig_split;
    typedef Sim_Model::migration_piecewise<mig_const,mig_split> split_pop0;
    float scale_factor = 1.0f;                                          //entire simulation can be scaled up or down with little to no change in resulting normalized SFS

    GO_Fish::allele_trajectories b;
    b.sim_input_constants.num_populations = 2;                          //number of populations
    b.sim_input_constants.num_generations = scale_factor*pow(10.f,3)+1; //1,000 generations

    Sim_Model::F_mu_h_constant codominant(0.5f);                        //dominance (co-dominant)
    Sim_Model::F_mu_h_constant outbred(0.f);                            //inbreeding (outbred)
    Sim_Model::F_mu_h_constant mutation(pow(10.f,-9)/scale_factor);     //per-site mutation rate 10^-9

    int N_ind = scale_factor*pow(10.f,4);                               //initial number of individuals in population
    dem_const pop0(N_ind);
    dem_const pop1(0);
    dem_pop_const gen0(pop0,pop1,1);                                    //intial population size of N_ind for pop 0 and 0 for pop 1
    dem_const pop0_final(2*N_ind);
    dem_pop_const gen1(pop0_final,pop1,1);
    init_expansion gen_0_1(gen0,gen1,1);                                //population 0 grows to size 2*N_ind
    exp_growth pop1_gen100((log(100.f)/(scale_factor*900.f)),0.05*N_ind,scale_factor*100);
    dem_pop_const_exp gen100(pop0_final,pop1_gen100,1);                 //population 1 grows exponentially from size 0.05*N_ind to 5*N_ind
    Sim_Model::demography_piecewise<init_expansion,dem_pop_const_exp> demography_model(gen_0_1,gen100,scale_factor*100);

    mig_const no_mig_pop0;
    mig_dir no_pop1_gen0(0.f,1,1,no_mig_pop0);
    mig_split create_pop1(1.f,0,1,no_pop1_gen0);                        //individuals from population 0 migrate to form population 1
    split_pop0 migration_split(no_mig_pop0,create_pop1,scale_factor*100);
    float mig = 1.f/(2.f*N_ind);
    mig_const mig_prop(mig,b.sim_input_constants.num_populations);      //constant and equal migration between populations
    Sim_Model::migration_piecewise<split_pop0,mig_const> mig_model(migration_split,mig_prop,scale_factor*100+1);

    float gamma = -4;                                                   //effective selection
    Sim_Model::selection_constant weak_del(gamma,demography_model,outbred);

    b.sim_input_constants.compact_interval = 30;                        //compact interval
    b.sim_input_constants.num_sites = 100*2*pow(10.f,7);                //number of sites
    int sample_size = 1001;                                             //number of samples in SFS

    int num_iter = 50;                                                  //number of iterations
    Spectrum::SFS my_spectra;

    cudaEvent_t start, stop;                                            //CUDA timing functions
    float elapsedTime;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    float avg_num_mutations = 0;
    float avg_num_mutations_sim = 0;
    std::vector<std::vector<float> > results(num_iter);                 //storage for SFS results
    for(int j = 0; j < num_iter; j++){ results[j].reserve(sample_size); }

    for(int j = 0; j < num_iter; j++){
        if(j == num_iter/2){ cudaEventRecord(start, 0); }               //use 2nd half of the simulations to time simulation runs + SFS creation

        b.sim_input_constants.seed1 = 0xbeeff00d + 2*j;                 //random number seeds
        b.sim_input_constants.seed2 = 0xdecafbad - 2*j;
        GO_Fish::run_sim(b, mutation, demography_model, mig_model, weak_del, outbred, codominant, Sim_Model::bool_off(), Sim_Model::bool_off());
        Spectrum::site_frequency_spectrum(my_spectra,b,0,1,sample_size);

        avg_num_mutations += ((float)my_spectra.num_mutations)/num_iter;
        avg_num_mutations_sim += b.maximal_num_mutations()/num_iter;
        for(int i = 0; i < sample_size; i++){ results[j][i] = my_spectra.frequency_spectrum[i]; }
    }

    elapsedTime = 0;
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime, start, stop);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);

    //output SFS simulation results
    std::cout<<"SFS :"<<std::endl<< "allele count\tavg# mutations\tstandard dev\tcoeff of variation (aka relative standard deviation)"<< std::endl;
    for(int i = 1; i < sample_size; i++){
        double avg = 0;
        double std = 0;
        float num_mutations;
        for(int j = 0; j < num_iter; j++){ num_mutations = b.num_sites() - results[j][0]; avg += results[j][i]/(num_iter*num_mutations); }
        for(int j = 0; j < num_iter; j++){ num_mutations = b.num_sites() - results[j][0]; std += 1.0/(num_iter-1)*pow(results[j][i]/num_mutations-avg,2); }
        std = sqrt(std);
        std::cout<<i<<"\t"<<avg<<"\t"<<std<<"\t"<<(std/avg)<<std::endl;
    }

    std::cout<<"\nnumber of sites in simulation: "<< b.num_sites() <<"\ncompact interval: "<< b.last_run_constants().compact_interval;
    std::cout<<"\naverage number of mutations in simulation: "<<avg_num_mutations_sim<<"\naverage number of mutations in SFS: "<<avg_num_mutations<<"\ntime elapsed (ms): "<< 2*elapsedTime/num_iter<<std::endl;
}



int main(int argc, char **argv) { run_validation_test(); }

Below is the corresponding δaδi code:

import time

import numpy
from numpy import array
import dadi

def prior_onegrow_mig((nu1F, nu2B, nu2F, m, Tp, T), (n1,n2), pts):
    """
    Model with growth, split, bottleneck in pop2, exp recovery, migration

    nu1F: The ancestral population size after growth. (Its initial size is
          defined to be 1.)
    nu2B: The bottleneck size for pop2
    nu2F: The final size for pop2
    m: The scaled migration rate
    Tp: The scaled time between ancestral population growth and the split.
    T: The time between the split and present

    n1,n2: Size of fs to generate.
    pts: Number of points to use in grid for evaluation.
    """
    # Define the grid we'll use
    xx = yy = dadi.Numerics.default_grid(pts)
    gamma_model = -2;
    # phi for the equilibrium ancestral population
    phi = dadi.PhiManip.phi_1D(xx, gamma=gamma_model)
    # Now do the population growth event.
    phi = dadi.Integration.one_pop(phi, xx, Tp, nu=nu1F, gamma=gamma_model)

    # The divergence
    phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
    # We need to define a function to describe the non-constant population 2
    # size. lambda is a convenient way to do so.
    nu2_func = lambda t: nu2B*(nu2F/nu2B)**(t/T)
    phi = dadi.Integration.two_pops(phi, xx, T, nu1=nu1F, nu2=nu2_func, m12=m, m21=m, gamma1=gamma_model, gamma2=gamma_model)

    # Finally, calculate the spectrum.
    sfs = dadi.Spectrum.from_phi(phi, (n1,n2), (xx,yy))
    return sfs

def runModel(nu1F, nu2B, nu2F, m, Tp, T):

    ns = (0,1001)
    print 'sample sizes:', ns

    # These are the grid point settings will use for extrapolation.
    pts_l = [110,120,130]

    func = prior_onegrow_mig
    params = array([nu1F, nu2B, nu2F, m, Tp, T])

    # Make the extrapolating version of the demographic model function.
    func_ex = dadi.Numerics.make_extrap_func(func)

    # Calculate the model AFS.
    model = func_ex(params, ns, pts_l)
    model = model.marginalize((0,))
    return model
    

time1 = time.time()

model = runModel(2, 0.05, 5, 1, 0.005, 0.045)

time2 = time.time()
print "total runtime (s): " + str(time2-time1)
model = model/model.S()

for i in range(0,len(model)):
    print model[i]

Below is the makefile to create executable GOFish from examples/example_dadi/run.cu, each line is documented by the top part of the makefile:

Tip: The makefile below compiles machine code explicitly for generation 3.0 and 5.2 GPUs and uses just in time (JIT) compilation for everything else (lowest GPU generation which works for 3P is 3.0). Compilation (and program execution) will be faster if compiling for your specific GPU.

e.g. if running a Tesla K20 or Tesla K40, then the corresponding GPU generation is 3.5: all the current --generate-code arch=##,code=## flags can be deleted and replaced with --generate-code arch=compute_35,code=sm_35.

# Description of Mac/Linux/Unix Makefile for example_dadi. 
#
#############################
# build_path := Where to build program and put executable (note: folder must already exist)
# api_path_source := Location of API source folder
# api_path_include := Location of API include folder
# EXEC_FILE := Name of executable 
#
# NVCC := Compiler path, in this case nvcc is in $PATH
# CFLAGS := Compiler Flags: optimize most, fast math, add API include folder to include search path, equivalent to --relocatable-device-code=true --compile
# CODE := GPU types for which to build explicitly (I have a NVIDIA GTX 780M and 980) https://developer.nvidia.com/cuda-gpus, creates machine code for code=sm_30 (780) and code=sm_52 (980) and virtual architectures for all other generations which can be compiled JIT - code=compute_30 for generations between (3.0,5.0) and code=compute_50 for generations (5.0 and up) http://docs.nvidia.com/cuda/cuda-compiler-driver-nvcc/index.html#options-for-steering-gpu-code-generation
#
# object_names := Objects required for executable
# objects := Prepends build path to object names
#
# all := Command 'make' or 'make all' builds executable file $(EXEC_FILE) in location $(build_path)/
#
# ##### Object Dependencies Lists #####
# If one of the files on the right hand side of the : changes, the corresponding object must be recompiled.
# This is a developer's makefile - it assumes there will be changes to the GOFish API files, if not
# all the non-*.cu (non source file) dependencies are unnecessary.
# ##### End Object Dependencies Lists #####
#
# $(objects) := Make target all objects
#  Compile source code into objects, $< := dependencies from Object Dependencies Lists, $@ := object in $(objects)
#
# $(build_path)/$(EXEC_FILE) := Make target executable EXEC_FILE which depends on all objects
#  Link objects into executable EXEC_FILE, $@ := $(build_path)/$(EXEC_FILE)
#
# .PHONY := Defines 'all' and 'clean' as not true targets (i.e. don't remake executable if can't find files called 'all' or 'clean')
#
# clean := Action to perform for command 'make clean'
#  Remove all objects and EXEC_FILE from build_path
#############################

build_path = ../example_dadi
api_path_source = ../../3P/_internal
api_path_include = ../../3P
EXEC_FILE = GOFish

NVCC = nvcc
CFLAGS = -O3 --use_fast_math -I $(api_path_include)/ -dc
CODE = --generate-code arch=compute_30,code=sm_30 --generate-code arch=compute_52,code=sm_52 --generate-code arch=compute_30,code=compute_30 --generate-code arch=compute_52,code=compute_52

object_names = run.o spectrum.o shared.o go_fish_impl.o
objects = $(addprefix $(build_path)/,$(object_names))

all:$(build_path)/$(EXEC_FILE)

##### OBJECT DEPENDENCIES #####
$(build_path)/run.o: run.cu $(api_path_source)/shared.cuh $(api_path_source)/go_fish_impl.cuh $(api_path_include)/go_fish_data_struct.h $(api_path_include)/go_fish.cuh $(api_path_include)/spectrum.h
$(build_path)/spectrum.o: $(api_path_source)/spectrum.cu $(api_path_include)/spectrum.h $(api_path_source)/shared.cuh $(api_path_include)/go_fish_data_struct.h
$(build_path)/shared.o: $(api_path_source)/shared.cu $(api_path_source)/shared.cuh
$(build_path)/go_fish_impl.o: $(api_path_source)/go_fish_impl.cu $(api_path_source)/go_fish_impl.cuh $(api_path_source)/shared.cuh
##### END OBJECT DEPENDENCIES #####

$(objects):
   $(NVCC) $(CODE) $(CFLAGS) $< -o $@

$(build_path)/$(EXEC_FILE): $(objects)

   $(NVCC) $(CODE) $(objects) -o $@

.PHONY: all clean

clean:
   rm -f $(objects) $(build_path)/$(EXEC_FILE)