sashagusev

# ---
# Code to estimate and simulate selection on healthy individuals (controls) for a threshold trait
# Updated with alternative simulation setup (producing the same result) and averaging to get smoother curves
# ---

library("VGAM")

# --- Estimate using the Breeder's Equation for threshold traits (Walsh & Lynch 2018; Dempster & Lerner 1950)
# We will select on being "healthy" where 1-prev of the population is healthy
ThresholdSelection <- function( FractionHealthy, Heritability ) {

sashagusev

library("RColorBrewer")
library("umap")
library("Rtsne")

# Number of markers
M = 2e3
# Number of people
N = 1e3

# frequencies in the ancestral population

sashagusev

par(mfrow=c(1,5))
par(oma=c(0.5,0.5,0.5,0.5))
par(mar=c(0.5,0.5,1,0.5))

# fst parameters
# frequency in one population
p_b = 0.2
# scan of frequencies in the second population
p_a = seq(0.01,0.99,by=0.01)
# estimate using Hudson's Fst

sashagusev

par(mfrow=c(1,2))
seeds = 10
x = seq(1,10,by=0.001)
all_u_coefs = seq(-5,5,by=0.1)

# lin reg
all_bs = rep(NA,length(all_u_coefs))
for ( i in 1:length(all_u_coefs) ) {
  u_coef = all_u_coefs[i]
  cur_r = rep(NA,seeds)

sashagusev

# simple factor analysis example
# ---

# lavaan for factor analysis
library('lavaan')

# --- parameters
# number of people
N = 1e3
# number of factors

sashagusev

# --- admixture analysis

# sample size
N = 5e3

# mate correlation
AM_cor = 0.75

# variance in ancestry transmission
genetic_sd = 0.025

sashagusev

# number of samples in each study
N = 2e3
# number of studies to run
seeds = 50

# vectors for storing results
est_a = rep(NA,seeds)
est_c = rep(NA,seeds)
est_e = rep(NA,seeds)
est_y1 = rep(NA,seeds)

sashagusev

set.seed(42)

# sample size
N = 10e3

# generate two groups
group = rbinom( N , 1 , 0.5 )

# generate group-specific environments
shared_env = rnorm(N,group*-15,1)

sashagusev

args = as.numeric(commandArgs(trailingOnly=TRUE))
# N,M,var_direct,var_indirect

# code for implementation of REML:
# https://github.com/sashagusev/SKK-REML-sim/blob/master/func_reml.R

source("reml_func.R")

# Num individuals
N = 1e3

sashagusev

# Num individuals
N = 10e3
# Num markers
M = 10

# Assortative Mating
AM_cor = 0.75
# Minor allele frequency
maf = 0.5
# Heritability