#!/bin/bash

### MAKE SURE YOU RUN IN BASH:
bash


### Make a few replicates of random phenotypes from common or rare SNPs:
### 1. Randomly choose 1000 causal variants, either rare (MAF<0.05) or common (MAF>=0.05) from the MAF distribution, which I already obtained using plink.
### Doing this 4 times (to keep it short). How many replicates might you want to simulate? 
### This is one way to randomly choose some variants, and you could do the same thing in lots of ways. You can change things as well - how might you change the number of CVs? Why might you want to vary that number?
for ii in {1..4}
do
    awk '$5>=0.05' simulation_2/5K.unrel0.05.freq.frq | perl -MList::Util -e 'print List::Util::shuffle <>' | head -1000 | awk '{print $2}' > simulation_2/common.rep"$ii".txt 
    awk '$5<0.05 && 5>0' simulation_2/5K.unrel0.05.freq.frq | perl -MList::Util -e 'print List::Util::shuffle <>' | head -1000 | awk '{print $2}' > simulation_2/rare.rep"$ii".txt
done


### 2. Use GCTA to create phenotypes from those causal SNPs within your sample.
### When no effect size is specified, allelic effects will be drawn from a standard normal distribution in GCTA. You might want to test something else (we will try this in the third simulation exercise)
for ii in {1..4}
do
    gcta64 --bfile simulation_2/5K.unrel0.05 --simu-qt --simu-hsq 0.5 --simu-causal-loci simulation_2/common.rep"$ii".txt --out simulation_2/pheno.common.rep"$ii".txt &
    gcta64 --bfile simulation_2/5K.unrel0.05 --simu-qt --simu-hsq 0.5 --simu-causal-loci simulation_2/rare.rep"$ii".txt --out simulation_2/pheno.rare.rep"$ii".txt & 
done
wait


### 3. Use GCTA to test the influence of causal variant MAF (which you've simulated two different conditions) on SNP-heritability
ii=1
gcta64 --grm simulation_2/5K.unrel0.05 --pheno simulation_2/pheno.common.rep"$ii".txt.phen --reml --reml-no-lrt --out simulation_2/hsq.common.rep"$ii" --thread-num 4
gcta64 --grm simulation_2/5K.unrel0.05 --pheno simulation_2/pheno.rare.rep"$ii".txt.phen --reml --reml-no-lrt --reml-priors 0.2 0.8 --out simulation_2/hsq.rare.rep"$ii" --thread-num 4
# how would you do this efficiently for every simulated phenotype?


### 4. How might you test the influence of sample size?  Try varying to 2.5K instead of 5K:
perl -MList::Util -e 'print List::Util::shuffle <>' simulation_2/5K.unrel0.05.grm.id | head -2500 > simulation_2/2.5K.unrel0.05.id
gcta64 --grm simulation_2/5K.unrel0.05 --keep simulation_2/2.5K.unrel0.05.id --pheno simulation_2/pheno.common.rep"$ii".txt.phen --reml --reml-no-lrt --out simulation_2/hsq.common.rep"$ii".2.5K --thread-num 4                         



grep -P 'V\(G\)\/Vp' simulation_2/*hsq  ### How could you test if there is a difference in hsq estimates based on MAF of the (simulated) causal variants?

