----------------------------
ESTIMATING SNP HERITABILITY-
----------------------------

At this point, the dataset is stored in data.speed, data.bim and data.fam, while prune.keep provides a list of unrelated (and population homogeneous) individuals and covar.covar contains covariates. If any highly-associated SNPs were identified, these are listed in top.in. Suppose phenotypes are stored in phen.pheno. Again, always keep an eye on the screen output from LDAK.

------------------------------------------------------------------------------

#Prepare to compute SNP weightings
ldak --cut-weights sections --speed data --keep prune.keep
#The details of each section will be stored in sections/section.details

#Compute weightings for each section - for this example, suppose there are 200 sections
#The exact number of sections is provided in sections/section.number
#When data come from multiple cohorts (here 3), using Subset Options guards against genotype errors
#For this, we create the files keep1, keep2 & keep3, containing individuals from Cohorts A, B & C
awk '(NR==FNR){arr[$1];next}($1 in arr){print $1, $2}' prune.keep cohortA/phase1.sample > keep1
awk '(NR==FNR){arr[$1];next}($1 in arr){print $1, $2}' prune.keep cohortB/phase1.sample > keep2
awk '(NR==FNR){arr[$1];next}($1 in arr){print $1, $2}' prune.keep cohortC/phase1.sample > keep3
for j in {1..200}; do
ldak --calc-weights sections --speed data --subset-number 3 --subset-prefix keep --section $j
done
#Typically sections complete in under 2 hours (most within 10 minutes)

#Join weightings across sections
ldak --join-weights sections --speed data
#Weights will be stored in sections/weightsALL

------------------------------------------------------------------------------

We will compute kinships for each chromosome separately; this will allow us to incorporate a test for inflation of h2SNP due to cryptic relatedness. It is normally not necessary here to include highly-associated SNPs as covariates, but if you do, make sure not to double count their contribution. If the trait is quantitative and individuals span multiple cohorts, we can also test for inflation due to genotyping errors.

------------------------------------------------------------------------------

#Get a list of SNPs on each chromosome
awk '{print $4 > "chr"$1""}' data.bim

#Compute kinships for each chromosome
ldak --cut-kins kinships --speed data --by-chr YES
for j in {1..22}; do
ldak --calc-kins kinships --speed data --partition $j --weights sections/weights.all --power -0.25
done

#Join to get genome-wide kinships
ldak --join-kins kinships
#Genome-wide kinships will be saved with prefix kinships/kinship.all

#Construct kinships for Chromosomes 1-3, 4-7, 8-11 & 12-22
rm listA; for j in {1..3}; do echo "kinships/kinships.$j" >> listA; done
rm listB; for j in {4..7}; do echo "kinships/kinships.$j" >> listB; done
rm listC; for j in {8..11}; do echo "kinships/kinships.$j" >> listC; done
rm listD; for j in {12..22}; do echo "kinships/kinships.$j" >> listD; done
for j in {A,B,C,D}; do
ldak --add-grm kins$j --mgrm list$j
done

#Test each quarter separately and combined
for j in {A,B,C,D}; do
ldak --reml quad$j --grm kins$j --pheno phen.pheno --covar covar.covar --keep prune.keep \
--top-preds top.in --speed data
done
echo -e "kinsA\nkinsB\nkinsC\nkinsD" > listABCD
ldak --reml quadALL --mgrm listABCD --pheno phen.pheno --covar covar.covar --keep prune.keep \
--top-preds top.in --speed data

#The results are stored in quadA.reml, quadB.reml, quadC.reml, quadD.reml and quadALL.reml
#The estimated inflation is (h2A + h2B + h2C + h2D - h2ALL)/3 - useful to view as a percentage
grep Her_K quad{A,B,C,D,ALL}.reml | awk '(NR<=4){sum+=$2}(NR>4){sum2+=$2}END{I=(sum-sum2)/3;print "Inflation:", I, "=", I/sum2*100,"%"}'

#Test for inflation due to genotyping errors
ldak --he inflation --grm kinships/kinships.all --pheno phen.pheno --covar covar.covar \
--keep prune.keep --subset-number 3 --subset-prefix keep --top-preds top.in --speed data
#The results are stored in inflation.he and inflation.he.compare

------------------------------------------------------------------------------

If satisfied that there is minimal inflation due to cryptic relatedness (and genotyping error if tested), we are ready to estimate SNP heritability. This command assumes single-SNP regression identified highly-associated SNPs; if not, ignore the final two arguments.

------------------------------------------------------------------------------

#Estimate SNP heritability
ldak --reml final --grm kinships/kinships.all --pheno phen.pheno --covar covar.covar \
--keep prune.keep --top-preds top.in --speed data