These scripts explain how to estimate functional enrichments and genetic correlations

#####################

Estimate functional enrichments. Now we need to create a (multi-part) tagfile, which contains $q_j I_{jk} + \sum_{l \in N_j} q_j I_{lk} r^2_{jl}$ for each SNP and category. To work out which SNPs are in each of the 24 annotations, we use the genefiles in \verb|annotations.zip|, available at \blue[http://ldak.org/annotations] (these are modified versions of the bedfiles provided at \blue[https://data.broadinstitute.org/alkesgroup/LDSCORE]).


#First work out which (reference panel) SNPs are in each annotation
#Will save these in ann_snps.1, ..., ann_snps.24 (the SNP lists must have consecutive suffixes)
for j in {1..24}; do
./ldak5.linux --cut-genes ann$j --bfile ref --genefile ann$j.genefile --ignore-weights YES
mv ann$j/genes.predictors.used ann_snps.$j
done

#Now compute a tagfile - as the categories overlap, we use --annotation-number and --annotation-prefix
./ldak5.linux --calc-tagging alz_ldak.ann --bfile ref --extract alz.txt --weights alz.weights \
--power -.25 --window-cm 1 --annotation-number 24 --annotation-prefix ann_snps.

#Perform the regression (adding --exclude, --check-sums NO, --genomic-control YES as required)
./ldak5.linux --sum-hers alz_ldak.ann --tagfile alz_ldak.ann.tagging --summary alz.txt --exclude alz.excl

#The category enrichments and SDs are in Columns 5 and 6 of alz_ldak.ann.enrich
#Suppose we have results for alz, years and height; to compute inverse-variance weighted averages enrichments use
for i in {1..24}; do 
grep "Enrich_A$i " {alz,years,height}_ldak.ann.enrich | awk -v i=$i '{a+=$5/$6/$6;b+=1/$6/$6}END\
{print "Annotation:", i, "Av. Enr.:", a/b, "SD:", 1/sqrt(b), "Num. Traits:", NR}'
done

#For the PRS, we must save the average heritability fractions of each annotation and the base (use the .share files)
rm props.ldak
for i in {1..24}; do 
grep "Share_A$i " {alz,years,height}_ldak.ann.share | awk -v i=$i '{a+=$2/$3/$3;b+=1/$3/$3}END\
{print a/b, 1/sqrt(b)}' >> props.ldak
done
grep Share_Base {alz,years,height}_ldak.ann.share | awk -v i=$i '{a+=$2/$3/$3;b+=1/$3/$3}END\
{print a/b, 1//sqrt(b)}' >> props.ldak

#To instead use the GCTA Model, use --ignore-weights YES and --power -1 when computing the tagfile
#Note that to copy Finucane et al, you must use all 52 annotations (not just the first 24)
./ldak5.linux --calc-tagging alz_ldak.ann --bfile ref --extract alz.txt --ignore-weights YES --power -1 --window-cm 1 --annotation-number 24 --annotation-prefix ann_snps.

#If categories partition the genome (LDAK allows some overlap), should use --partition-number and --partition-prefix
#In this example, we construct lists part.1 and part.2 containing coding (Annotation 1) and non-coding SNPs
cp ann_snps.1 part.1
awk '(NR==FNR){arr[$1];next}!($2 in arr){print $2}' part.1 ref.bim > part.2
./ldak5.linux --calc-tagging alz_ldak.coding --bfile ref --extract alz.txt --weights alz.weights \
--power -.25 --window-cm 1 --partition-number 2 --partition-prefix part.

#Note that we can achieve the same using
./ldak5.linux --calc-tagging alz_ldak.coding --bfile ref --extract alz.txt --weights alz.weights \
--power -.25 --window-cm 1 --partition-number 1 --partition-prefix part. --background YES

#####################

Estimate genetic correlations. We calculate a (1-part) tagfile, using only SNPs common to the two traits (here alzheimers and years education). Now the direction of effects matters, so we recommend using only non-ambiguous SNPs (alleles A \& C, A \& G or C \& T), unless very confident the alleles are consistently aligned (note that if you include ambiguous SNPs, you must add \verb|--allow-ambiguous YES| when performing the regression).


#Find the overlap of non-ambiguous SNPs
awk < alz.txt '(($2=="A"&&$3=="C") || ($2=="A"&&$3=="G") || ($2=="C"&&$3=="T") || ($2=="C"&&$3=="A") \
||($2=="G"&&$3=="A") || ($2=="T"&&$3=="C")){print $1}' > alz.nonamb
awk < years.txt '(($2=="A"&&$3=="C") || ($2=="A"&&$3=="G") || ($2=="C"&&$3=="T") || ($2=="C"&&$3=="A") \
||($2=="G"&&$3=="A") || ($2=="T"&&$3=="C")){print $1}' > years.nonamb
awk '(NR==FNR){arr[$1];next}($1 in arr){print $1}' alz.nonamb years.nonamb > overlap.nonamb

#Will exclude SNPs in the exclusion lists for either trait
cat alz.excl years.excl | sort | uniq > alz_years.excl

#Calculate LDAK weights for the overlap
for j in {1..22}; do
./ldak5.linux --cut-weights alz_years_chr$j --bfile ref --extract overlap.nonamb --chr $j
./ldak5.linux --calc-weights-all alz_years_chr$j --bfile ref --extract overlap.nonamb --chr $j
done
cat alz_years_chr{1..22}/weights.short > alz_years.weights

#Calculate LDAK tagging file
./ldak5.linux --calc-tagging alz_years_ldak --bfile ref --extract overlap.nonamb --weights alz_years.weights \
--power -.25 --window-cm 1

#Perform the regression (adding --exclude, --check-sums NO, --genomic-control YES as required)
./ldak5.linux --sum-cors alz_years_ldak --tagfile alz_years_ldak.tagging --summary alz.txt --summary2 years.txt \
--exclude alz_years.excl

#The estimate of genetic correlation is in alz_years_ldak.cors

#####################

Comparing heritability models. The easiest way to compare two heritability models is by likelihood. By default, LDAK excludes SNPs with zero weight when computing tagfiles. However, when comparing heritability models, it is important that each tagfile contain the same SNPs. Therefore, when computing an LDAK tagfile, add \verb|--reduce NO|. The second way to compare two heritability models is by creating a merged tagfile, then LDAK will estimate the share of heritability allocated to each.


#Comparing the GCTA and LDAK Models by likelihood

#Compute (standard) GCTA tagfile and the (unreduced) LDAK tagfile (as before, could divide by chromosome)
./ldak5.linux --calc-tagging alz_gcta --bfile ref --extract alz.txt --ignore-weights YES --power -1 --window-cm 1
./ldak5.linux --calc-tagging alz_ldak.non$j --bfile ref --extract alz.txt --weights alz.weights \
--power -.25 --window-cm 1 --reduce NO

#Perform the regression for each (adding --exclude, --check-sums NO, --genomic-control YES as required)
./ldak5.linux --sum-hers alz_gcta --tagfile alz_gcta.tagging --summary alz.txt --exclude alz.excl
./ldak5.linux --sum-hers alz_ldak.non --tagfile alz_ldak.non.tagging --summary alz.txt --exclude alz.excl

#The likelihoods and likelihood ratio test statistics are in alz_gcta.extra and alz_ldak.non.extra

#Comparing the GCTA and LDAK Models by creating a merged tagfile
echo "alz_gcta.tagging
alz_ldak.non.tagging" > list.txt
./ldak5.linux --merge-tagging alz_both --taglist list.txt
./ldak5.linux --sum-hers alz_both --tagfile alz_both.tagging --summary alz.txt --exclude alz.excl

#The share of heritability allocated to the LDAK Model (p) is the value Share_P2 in alz_both.share