LDAK is able to perform set-based association analysis. We typically use this feature for gene-based testing. However, an alternative is to perform a chunk-based association analysis (e.g., where each chunk is a 100kb section of the genome)
To perform a set-based analysis involves three steps: the first divides predictors into sets, the second tests each set for association with the phenotype, while the third joins results across sets. If your sets correspond to genes, the first step requires gene annotations. On Resources, you can download gene annotations for the human genome. The page also explains how I created these files, which might be useful if you wish to construct annotations yourself (e.g., for non-human species).
Please note that the commands on this page are designed for use with individual-level data; if you instead plan to use summary statistics, you should instead read the instructions for LDAK-GBAT.
Always read the screen output, which suggests arguments and estimates memory usage.
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Divide predictors into sets:
The main argument is --cut-genes <folder>.
This requires the options
--bfile/--gen/--sp/--speed <datastem> or --bgen <datafile> - to specify the genetic data files (see File Formats).
--genefile <genefile>, --chunks <float> or --chunks-bp <integer> - to specify the sets of predictors (see below).
You can use --extract <extractfile> and/or --exclude <excludefile> to restrict to a subset of predictors (for more details, see Data Filtering).
There are three ways to specify predictor sets. The first is to provide gene annotations using --genefile <genefile>. The first four columns of <genefile> should specify the name, chromosome, start and end basepairs of each gene. If there is a fifth column, its value should be either + or -, indicating whether the gene is on the forward or backwards strand. Note that the basepairs should provide 0-start, half-open positions (this is the format used by Browser Extensible Data). For example, if Gene ABC is on Chromosome 7 and contains basepairs 1-10, inclusive, the corresponding row of <genefile> would be "ABC 7 0 10" (if providing strand orientation, the row would be either "ABC 7 0 10 +" or "ABC 7 0 10 -").
By default, LDAK will only consider predictors inside genes. To change this, you can use --gene-buffer <integer>, --up-buffer <integer> or --down-buffer <integer>. For example, if you use --gene-buffer 500, each set will contain predictors inside or within 500bp of a gene. Note that you can only use --up-buffer and --down-buffer if <genefile> provides the orientation of the gene (i.e., has five columns). If you add --overlap NO, LDAK will ensure that each predictor is only allocated to the nearest gene.
Instead of providing gene annotations, you can use --chunks <float> to divide predictors into chunks of weighting approximately <float> (you can use --weights <weightsfile> to provide predictor weightings, else each predictor will get weighting one). Alternatively, you can use --chunks-bp <integer> to divide predictors into chunks of length <integer> base pairs. By default, LDAK will construct overlapping chunks. For example, if you use --chunk-bp 1000, the first three sets on a chromosome will span basepairs 1-1000, basepairs 501-1500 and basepairs 1001-2000. However, if you add --overlap NO, then the first three sets on a chromosome will instead span basepairs 1-1000, basepairs 1001-2000 and basepairs 2001-3000.
If you add --partition-length or --by-chr YES, LDAK will divide the genome into partitions, which can be analysed separately. You can exclude small genes using --min-weights (for example, using --min-weights 5 will, by default, exclude genes containing fewer than 5 predictors).
If you have predictor p-values (e.g., from Single-Predictor Analysis), you can provide these using the option --pvalues <pvalues>. The file <pvalues> should have one row per predictor and two columns, which provide predictor names then p-values. LDAK will then report the minimum p-value for each gene/chunk, and also Bonferroni-corrected p-values.
The main output file is <folder>/genes.details, which details which predictors are in each set (and perhaps the minimum and Bonferroni-corrected p-value for each set). The file <folder>/genes.predictors lists which predictors are within one or more sets. If you provided gene annotations, then the distance from each predictor to the nearest gene will be reported in the file <folder>/genes.distance (negative values indicate the predictor is within a gene).
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Test each set for association:
The main argument is --calc-genes-reml <folder>.
This requires the options
--pheno <phenofile> - to specify phenotypes (in PLINK format). Samples without a phenotype will be excluded. If <phenofile> contains more than one phenotype, specify which should be used with --mpheno <integer> or --pheno-name <string> (the latter requires that <phenofile> has a header row), or use --mpheno ALL to analyse all phenotypes.
--bfile/--gen/--sp/--speed <datastem> or --bgen <datafile> - to specify the genetic data--exclude <excludefile> files (see File Formats).
--power <float> - to specify how predictors are scaled. In general, we recommend using --power -0.25.
If you used --partition-length or --by-chr YES when dividing predictors into sets, you must run this step multiple times, using --partition <integer> to specify which partition to analyse.
You can use --keep <keepfile> and/or --remove <removefile> to restrict to a subset of samples, and --extract <extractfile> and/or --exclude <excludefile> to restrict to a subset of predictors (for more details, see Data Filtering). Please note that the predictor filterings used in this step must be match those used when dividing predictors into sets.
You can use --covar <covarfile> or --factors <factorfile> to provide quantitative or categorical covariates (in PLINK format) as fixed effects in the regression (if <covarfile> contains multiple covariates, you can specify a subset using either --covar-names <string> or --covar-numbers <string>). Furthermore, if the data contains predictors that are definitely associated with the phenotype, you can include these as fixed effects using --top-preds <toppredsfile>.
You can use --grm <kinfile> to provide a kinship matrix. LDAK will then include this when performing REML (in addition to the kinship matrix calculated across the set of predictors being tested). We recommend doing this when your dataset contains related samples and/or population structure (it is an effective way to guard against confouding due to cryptic relatedness). To construct the kinship matrix, we suggest Calculating Kinships assuming a thinned version of the Uniform Model, that restricts to SNPs in approximate linkage equilibrium (see below).
You can add --permute YES in order to shuffle the predictor values. This option is useful if wishing to perform permutation analysis to determine the distribution of p-values or test statistics when there is no true signal.
To reduce runtime and increase computational stability, LDAK prunes the predictors in each set so that no pair of predictors has squared correlation above 0.98; to change this threshold, use --gene-prune <float>.
By default, LDAK will assign all predictors weighting one (equivalent to using --ignore-weights YES). You can provide your own weightings using --weights <weightsfile>, but we have not found this beneficial.
LDAK uses permutations to derive the null distribution of the likelihood ratio test statistic (copying the work of Listgarten el al.). By default, LDAK performs ten permutation for each predictor set; to change this, use --gene-permutations <integer>. By default, the permutations will vary across runs. However, you can prevent this using --random-seed <integer> (see Advanced Options for more details).
When analysing a binary phenotype, you can use --prevalence <float> to specify the population prevalence. LDAK will then additionally report heritability estimates on the liability scale.
Provisional association testing results will be saved in the file <folder>/remls.1 (or in the files <folder>/remls.1, <folder>/remls.2, ... if you have multiple partitions). Meanwhile, the estimated genetic contributions of each set will be saved in files with the prefix <folder>/prs.1 (or with prefixes <folder>/prs.1, <folder>/prs.2, ... if you have multiple partitions).
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Join results across sets:
The main argument is --join-genes-reml <folder>.
This requires no options.
The final association testing results will be saved in <folder>/remls.all. For each set, this file reports the estimated heritability and two p-values (for reasons explained above, we recommend that you use the second p-value, obtained using permutations). The estimated genetic contributions of each set will be saved in SP format in the files <folder>/prs.all.sp, <folder>/prs.all.bim and <folder>/prs.all.fam (see Genetic File Formats for details), while the corresponding p-values are saved in <folder>/prs.all.pvalues. Note that to ensure predictor names are unique, LDAK will prefix gene names by their index.
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
Example:
Here we use the binary PLINK files human.bed, human.bim and human.fam, the phenotype quant.pheno, and the annotations file anns.txt from the Test Datasets. The file anns.txt contains annotations for human genes on Chromosomes 21 & 22
head -n 5 anns.txt
MIR3648-2 21 9825831 9826011
MIR3687-2 21 9826202 9826263
TEKT4P2 21 9907188 9968594
TEKT4P2 21 9915249 9968594
TPTE 21 10906186 10990943
For example, we see that the Gene TEKT4P2 lies on Chromosome 21 from basepairs 9907189 to 9968594. Note that the same gene can have more than one entry (as is the case here). You may prefer to reduce to just one annotation per gene, but in our experience, it is fine to keep them all.
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
1 - Basic gene-based analysis
We can run a gene-based analysis using the commands
./ldak.out --cut-genes genes --bfile human --genefile anns.txt
./ldak.out --calc-genes-reml genes --pheno quant.pheno --bfile human --power -0.25
./ldak.out --join-genes-reml genes
The final association testing results are saved in genes/remls.all, while estimates of the genetic contributions of each gene are saved in the files genes/remls.all.sp, genes/remls.all.bim and genes/remls.all.fam (the latter files can be used for Clumping). We can find the most significant gene by sorting based on the likelihood ratio test statistic (Column 10)
head -n 1 genes/remls.all
Gene_Name Gene_Chr Gene_Start Gene_End Length Heritability SD Null_Likelihood Alt_Likelihood LRT_Stat LRT_P_Raw LRT_P_Perm Her_Liability SD_Liability
sort -n -k 10 genes/remls.all |tail -n 1
ABCC13 21 15646120 15673692 20 0.096276 0.052390 -603.7353 -590.3308 26.8090 1.1230e-07 1.3031e-07 NA NA
We see that the most significant gene is ABCC13. Its estimated heritability is 0.10 (SD 0.05), while its test statistic is 27. There are two p-value, 1.1e-7 and 1.3e-7. We prefer the second p-value (Column 12), which was obtained via permutations (note that this p-value is stochastic, so will vary between runs).
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
2 - Gene-based analysis including a genome-wide kinship matrix
We recommend computing the kinship matrix assuming a thinned version of the Uniform Model, for which we use the scripts
./ldak.out --thin le --bfile human --window-prune .05 --window-cm 1
./ldak.out --calc-kins-direct le --bfile human --power -1 --extract le.in
We then repeat the analysis above, except adding this kinship matrix when testing each set of predictors.
./ldak.out --cut-genes genes2 --bfile human --genefile anns.txt
./ldak.out --calc-genes-reml genes2 --pheno quant.pheno --bfile human --power -.25 --grm le
./ldak.out --join-genes-reml genes2
The final association testing results are saved in genes2/remls.all, while estimates of the genetic contributions of each gene are saved in the files genes2/remls.all.sp, genes2/remls.all.bim and genes2/remls.all.fam (the latter files can be used for Clumping). Again we sort genes based on the likelihood ratio test statistic (Column 10)
head -n 1 genes2/remls.all
Gene_Name Gene_Chr Gene_Start Gene_End Length Heritability SD Null_Likelihood Alt_Likelihood LRT_Stat LRT_P_Raw LRT_P_Perm Her_Liability SD_Liability
sort -n -k 10 genes2/remls.all |tail -n 1
ABCC13 21 15646120 15673692 20 0.087546 0.050512 -599.7448 -588.0336 23.4224 6.5027e-07 1.4419e-04 NA NA
We see that the most significant gene remains ABCC13. However, now its estimated heritability is 0.09 (SD 0.05), its test statistic is 23 and its (permuted) p-value is 1.3e-4.
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
3 - Basic chunk-based analysis
We can run a chunk-based analysis using the commands
./ldak.out --cut-genes chunks --bfile human --chunks-bp 100000
./ldak.out --calc-genes-reml chunks --pheno quant.pheno --bfile human --power -.25
./ldak.out --join-genes-reml chunks
Here we have used 100kb, overlapping chunks. The final association testing results are saved in chunks/remls.all, while estimates of the genetic contributions of each chunk are saved in the files chunks/remls.all.sp, chunks/remls.all.bim and chunks/remls.all.fam (the latter files can be used for Clumping). Again we sort genes based on the likelihood ratio test statistic (Column 10)
head -n 1 chunks/remls.all
Gene_Name Gene_Chr Gene_Start Gene_End Length Heritability SD Null_Likelihood Alt_Likelihood LRT_Stat LRT_P_Raw LRT_P_Perm Her_Liability SD_Liability
sort -n -k 10 chunks/remls.all |tail -n 1
Chunk_19 21 15550000 15650000 25 0.192006 0.086320 -603.7353 -580.9097 45.6511 7.0653e-12 3.1972e-13 NA NA
We see that the most significant chunk is on Chromosome 21 from basepairs 15,550,000 to 15,650,000. Its estimated heritability is 0.19 (SD 0.09), while its test statistic is 46 and its (permuted) p-value is 3.2e-12.
