----------------------------
MERGING---------------------
----------------------------

Having imputed each cohort separately, we next merged imputed cohort data to form datasets. We stored data in Speed Format, a format which accommodates non-integer values and thus allows us to analyze expected allele counts. An alternative is to convert data to hard genotypes and store in PLINK binary format, which will have negligible impact if using only high-quality SNPs. This guide uses strict quality control, matching those we used for the main analyses. Specifically, we retained only bialellic autosomal SNPs which in all cohorts satisfied (expected) MAF>0.01 and r j > 0.99 (if imputed) or r2 type0 > 0.99 (if directly genotyped), where r j is our information score (Supplementary Figure 20) and r2 type0 is a metric computed by IMPUTE2. Suppose we are combining three cohorts, stored in the folders cohortA, cohortB and cohortC.

These instructions are designed to process the output from IMPUTE2, however, they should be easily modifiable for other imputation software (see www.ldak.org/file-formats for ways to incorporate different genotype formats). Always keep an eye on the screen output from LDAK to understand how datasets are being processed and for suggested options.

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

#Create badgeno.snps, listing SNPs which in any cohort were directly genotyped and had r2_type0<0.99
#Column 9 of the the info file indicates imputed (0) or genotyped (2); Column 12 provides r2_type0
rm badgeno.A badgeno.B badgeno.C
for j in {1..518}; do
awk < cohortA/chunk${j}_info '($9==2 && $12<.99){print $2}' >> badgeno.A;
awk < cohortB/chunk${j}_info '($9==2 && $12<.99){print $2}' >> badgeno.B;
awk < cohortC/chunk${j}_info '($9==2 && $12<.99){print $2}' >> badgeno.C;
done
cat badgeno.A badgeno.B badgeno.C | sort | uniq > badgeno.snps

#Merge cohorts by region, excluding SNPs in badgeno.snps, and filtering based on MAF and r_j
#'--common-preds YES' restricts to variants present in all (three) cohorts
#'--exclude-odd YES' excludes variants with alleles other than A, C, G, T (i.e., none SNPs)
#'--exclude-dups YES' excludes variants with matching positions (i.e., none diallelic SNPs)
#'--pass-all YES' ensures SNPs must pass QC in all cohorts
#Replace --make-speed with --make-bed to save data in Binary PLINK format

for j in {1..518}; do
echo -e "cohortA/chunk$j cohortA/phase1.sample" > list$j
echo -e "cohortB/chunk$j cohortB/phase1.sample" >> list$j
echo -e "cohortC/chunk$j cohortC/phase1.sample" >> list$j
chr=`awk -v j=$j '(NR==j){print $1}' allreg.txt`
ldak --make-speed merge$j --mgen list$j --oxford-single-chr $chr 
--exclude badgeno.snps --exclude-odd YES --exclude-dups YES 
--min-maf 0.01 --min-info .99 --pass-all YES --common-preds YES
done
#Merged data for Region j will be stored in speed format with prefix merge${j}

#Combine across regions
#Speed format files can be combined using cat; for PLINK data files, use --merge-list in PLINK
cat merge{1..518}.speed > data.speed
cat merge{1..518}.bim > data.bim
cp merge1.fam data.fam
head -n 1 merge1.stats > data.stats
for j in {1..518}; do
tail -q -n +2 merge$j.stats >> data.stats
done
awk < data.stats '(NR>1){print $1, $6}' > data.infos 
#The dataset is now saved in data.speed, data.bim and data.stats
#data.infos contains info scores (although these are not used for our main analysis)

#If hard genotypes are required, allele counts can be rounded to the nearest integer
ldak --make-bed hard --speed data --threshold 0.5
#Genotypes will be saved to hard.bed, hard.bim and hard.fam

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

To reliably estimate h2SNP , individuals must be unrelated, which we achieve by filtering based on alleleic correlations computed from a pruned subset of SNPs. Next we make the covariate file covar.covar, which has 33 columns: two IDs and sex (Columns 1, 2, 5 of data.fam), 20 dataset principal axes (Columns 3 to 22 of prune.vect) and 10 population axes (Columns 5, 7, ..., 23) of tg.profile. To obtain population axes, we use 1000 Genomes data; suppose these are stored in tg.bed, tg.bim and tg.fam. We additionally perform single-SNP analysis; this used as a sanity check, and also to identify any highly-associated SNPs (P < 10 −20 ) which we subsequently include (after pruning) as additional fixed-effect covariates.

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

#Thin SNPs, then use these to compute allelic correlations
#We ensure no SNPs within 1Mb have correlation squared >0.2
#Could also exclude SNPs in high-LD regions (e.g., the MHC)
ldak --thin prune --speed data --window-kb 1000 --window-prune 0.2
ldak --calc-kins-direct prune --speed data --extract prune.in --ignore-weights YES --power -0.25
#The list of pruned SNPs are contained in prune.in; the kinship has stem prune

#Filter relatedness (if you have multiple cohorts, for each cohort, then for all together), then obtain top 20 eigen-vectors
#By default, LDAK filters using the threshold c, where -c is smallest observed kinship
#To instead, filter based on a threshold of (say) 0.05, add --max-rel 0.05

#If you have just one cohort, use 
ldak --filter prune --grm prune
ldak --pca prune --grm prune --axes 20 --keep prune.keep

#If you have multiple cohorts, use
ldak --filter keepA --grm prune --keep cohortA/phase1.sample
ldak --filter keepB --grm prune --keep cohortB/phase1.sample
ldak --filter keepC --grm prune --keep cohortC/phase1.sample
cat keepA.keep keepB.keep keepC.keep > keepABC.keep
ldak --filter prune --grm prune --keep keepABC.keep
ldak --pca prune --grm prune --axes 20 --keep prune.keep

#Compute 10 population axes from 1000 Genomes data
ldak --calc-kins-direct tg --bfile tg --extract prune.in --ignore-weights YES --power -0.25
ldak --pca tg --grm tg --axes 10
ldak --calc-pca-loads tg --bfile tg --grm tg --pcastem tg
ldak --calc-scores tg --speed data --scorefile tg.load --allow-flips YES --keep prune.keep --power 0

#Identify highly-associated SNPs from single-SNP analysis (we define as P<1e-20)
ldak --linear linear --speed data --covar covar.covar --keep prune.keep --pheno phen.pheno
awk '(NR>1&&$7<1e-20){print $2}' linear.assoc > top.snps            #if using LDAK
#OR
plink --logistic hide-covar --out logistic --bfile hard --pheno phen.pheno \
--covar covar.covar --keep prune.keep
awk '(NR>1&&$9<1e-20){print $2}' logistic.assoc.logistic > top.snps  #if using PLINK

#Prune these - we used a correlation-squared threshold of 0.5 and window of 10Mb
ldak --thin top --speed data --extract top.snps --keep prune.keep \
--window-prune 0.5 --window-kb 10000

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