vQTL Analysis

This feature was developed in collaboration with Elham Assary.

A vQTL analysis tests whether predictors are associated with changes in the variance of a phenotype (this contrasts with a regular GWAS, which tests whether predictors are associated with changes in the phenotypic mean). One of the main motivations for performing a vQTL is to find evidence for gene-by-environment interactions (if a genetic variant's impact on a phenotype depends on environmental factors, then this can result in an association between the variant and the phenotypic variance).

Here we explain how to perform a vQTL study using the Deviation Regression Model (DRM; Marderstein et al., 2020)  We selected this method based on the comparisons of different vQTL methods by Zhang et al., 2024 Suppose the vector Y contains the phenotype and the vector G contains the genotypes of the SNP being tested. DRM first constructs an adjusted phenotype Y', whose ith element is Y'[i]=|Y[i]-M[i]|, where M[i] is the median phenotype across all individuals whose genotype equals G[i] (e.g., G[i]=1, then M[i] is the median of Y when restricted to individuals with genotype one). DRM then regresses the adjusted phenotype on the SNP genotypes.

When analysing (pseudo) autosomal predictors, LDAK tests all individuals jointly, but when analysing Chromosome X predictors, LDAK tests males and females separately. Please note that for Chromosome X predictors, males should be coded 0 and 2 (males with genotype 1 will be ignored).

Always read the screen output, which suggests arguments and estimates memory usage.
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

By default, LDAK will perform classical DRM, which is designed for use with homogeneous datasets. If your dataset contains related individuals, you should use --relatives <pairsfile> to specify pairs of related individuals. <pairsfile> should have five columns: Columns 1 & 2 should provide the two IDs for the first individual in each pair, Columns 3 & 4 should provide the two IDs for the second individual in each pair, while Column 5 should specify the relatedness between the pair. See Relatives File for more details and scripts for constructing this file.

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).

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>).

You can add --permute YES in order to shuffle the predictor values. This is useful if wishing to perform permutation analysis to see the distribution of p-values or test statistics when there is no true signal.

The main output file is <outfile>.DRM.all, which contains results from regressing Y' on G using a linear model, followed by results from regressing Y' on G using a general model (the latter includes estimates of the additive and dominant effects, as well as a p-value from testing whether either effect is non-zero).
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

This requires the options

--bfile <datastem> - to specify the genetic data files in PLINK format (see File Formats). Note that DRM requires hard genotypes (each predictor value is 0, 1 or 2), and therefore it is not possible to provide genetic data in other formats.

--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).

--split-par <string> - to define the pseudo-autosomal regions (PAR) predictors. Use --split-par hg19 to define them based on the hg19 genome assembly (PAR predictors are those with basepairs below 2,699,520 or above 154,931,044) or --split-par hg38 to define them based on the hg38 genomic assembly (PAR predictors are those with basepairs below 2,781,479 or above 155,701,383). Note that --split-par ALL (--split-par NONE) instructs LDAK to treat all (none) of the predictors as PAR.

--sexfile <sexfile> - to specify the sexes of individuals (in PLINK format). Males should be coded 1, while females should be coded 2. Note that if <sexfile> has more than three columns (e.g., if you are using a covariate file that includes sex), then it must have a header row and a column labelled "Sex".

The main output files are <outfile>.DRM.all, <outfile>.DRM.males and <outfile>.DRM.females, which contain results from regressing Y' on G using a linear model, followed by results from regressing Y' on G using a general model (the first file contains results for the PAR predictors, while the last two files contain results for the non-PAR predictors).
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

Example:

Here we generate and analyse Simulated Data. Please note that because we are using randomly-sampled genotypes and phenotypes, your DRM results will be slightly different to those shown below. In addition, for the Chromosome X analyses, LDAK will warn that some males have heterogeneous genotypes (this reflects that when simulating genotypes, LDAK does not force males to be homogeneous for Chromosome X SNPs).
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

1 - Simulate genotypes and phenotypes.

We generate genotypes for 3000 individuals and 23000 SNPs. Note that in the command below, the option --num-chr 23 ensures that LDAK generates Chromosome X SNPs, while the options --family-size 3 and --relatedness 1 tell LDAK that each set of three individuals should be 100% related (i.e., we end up with genotypes for 1000 identical triplets).

./ldak.out --make-snps fake --num-samples 3000 --num-snps 23000 --family-size 3 --relatedness 1 --num-chr 23

The simulated genotypes are stored in binary PLINK format in the files fake.bed, fake.bim and fake.fam, with simulated covariates (sex and age) in fake.covar. The file fake.pairs lists the related pairs; for example, the top of the file shows that the three individuals with IDs "FAMILY1 IND1", "FAMILY1 IND2" and "FAMILY1 IND3" are 100% related to each other.

head -n 3 fake.pairs
FAMILY1 IND1 FAMILY1 IND2 1.000000
FAMILY1 IND1 FAMILY1 IND3 1.000000
FAMILY1 IND2 FAMILY1 IND3 1.000000

We generate a quantitative phenotype with heritability 0.5 and 100 causal SNPs, with effect sizes sampled according to the Human Default Heritability Model.

./ldak.out --make-phenos fake --bfile fake --num-phenos 1 --her 0.5 --num-causals 100 --power -0.25

The simulated phenotypes are saved in fake.pheno.
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

2 - Classical DRM for autosomal SNPs.

We perform classical DRM using the command

./ldak.out --linear drm --bfile fake --pheno fake.pheno --DRM AUTOSOMES

The main output file is drm.DRM.all. This contains results for only 22000 SNPs, because the Chromosome X SNPs were excluded. p-values from the linear test are stored in Column 9; if we sort on this column, we see that the smallest p-value is 5e-9 (i.e., highly significant, even after correction for multiple testing, and a consequence of the high levels of relatedness in the dataset).

sort -g -k 9 drm.DRM.all | cut -f 1-9 | head -n 2
Chromosome Predictor Basepair A1 A2 Effect_LIN SE_LIN Z_LIN P_LIN
22 SNP21135 1350000 A C -9.2853e-02 1.5881e-02 -5.8469 5.0092e-09

p-values from the general test are stored in Column 19; if we sort on this column, we see that again the smallest p-value is 7e-9 (so again, highly significant).

sort -g -k 19 drm.DRM.all | cut -f 1-5,18,19 | head -n 2
Chromosome Predictor Basepair A1 A2 CHISQ_GEN P_GEN
6 SNP5429 4290000 A C 33.4225 7.4163e-09
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

3 - Mixed-Model DRM for autosomal SNPs.

We switch from classical to mixed-model DRM by providing LDAK with a list of relatives

./ldak.out --linear drm2 --bfile fake --pheno fake.pheno --DRM AUTOSOMES --relatives fake.pairs

The main output file is drm2.DRM.all. There are no longer highly significant SNPs from the linear test

sort -g -k 9 drm2.DRM.all | cut -f 1-9 | head -n 2
Chromosome Predictor Basepair A1 A2 Effect_LIN SE_LIN Z_LIN P_LIN
22 SNP21135 1350000 A C -9.2853e-02 1.9962e-02 -4.6515 3.2947e-06

Nor are there highly significant SNPs from the general test

sort -g -k 19 drm2.DRM.all | cut -f 1-5,18,19 | head -n 2
Chromosome Predictor Basepair A1 A2 CHISQ_GEN P_GEN
6 SNP5429 4290000 A C 20.1350 7.2163e-06

We can measure the impact of allowing for relatedness by comparing the mean chi-squared test statistics between the classical and mixed-model analyses

awk < drm.DRM.all '(NR>1){c++}{a+=$8^2;b+=$18;c++}END{print "Classical DRM Mean Statistics:", a/c, "(Linear) and ", b/c, "(General)"}'
Classical DRM Mean Statistics: 0.791831 (Linear) and 0.916983 (General)

awk < drm2.DRM.all '(NR>1){c++}{a+=$8^2;b+=$18;c++}END{print "Mixed-Model DRM Mean Statistics:", a/c, "(Linear) and ", b/c, "(General)"}'
Mixed-Model DRM Mean Statistics: 0.506079 (Linear) and 0.500975 (General)

4 - Classical DRM for Chromosome X SNPs.

We perform classical DRM using the command

./ldak.out --linear drm3 --bfile fake --pheno fake.pheno --DRM CHRX --split-par hg19 --sexfile fake.covar

Note that this command works because fake.covar has a column named "Sex". The screen output tells us that 269 of the 1000 Chromosome X SNPs are within the pseudo-autosomal regions, and that LDAK has performed three analysis: a joint analysis of the 269 PAR predictors (results saved in drm3.DRM.all), a male-only analysis of the 731 non-PAR predictors (results saved in drm3.DRM.males), and a female-only analysis of the 731 non-PAR predictors (results saved in drm3.DRM.females). As noted above, the second analysis produces a warning, because there are heterogeneous non-PAR predictors.
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

5 - Mixed-Model DRM for Chromosome X SNPs.

We switch from classical to mixed-model DRM by providing LDAK with a list of relatives

./ldak.out --linear drm4 --bfile fake --pheno fake.pheno --DRM CHRX --split-par hg19 --sexfile fake.covar --relatives fake.pairs

The main output files are drm4.DRM.all, drm4.DRM.males and drm4.DRM.females.