GRM

The GRM module builds the genomic relationship matrices (GRMs) used as random-effect covariance matrices by UVLMM, REMMA, fastREML and the other mixed-model modules. It also post-processes existing GRMs. The files it writes use the GCTA format, so a GRM built by fastGMA can be used in GCTA, and a GRM built by GCTA can be used in fastGMA.

There are two entry points:

As elsewhere in fastGMA, a GRM is referred to by its prefix: --grm output/test means the files output/test.grm.bin, output/test.grm.id and so on.

1. Building a GRM (--make-grm) #

GRM types #

Choose the matrix with --grm-type. Let \(x\) be the genotype dosage (0/1/2 copies of the A1 allele), \(p\) the A1 frequency and \(q = 1 - p\).

Under Hardy–Weinberg equilibrium the dominance code is uncorrelated with the additive code, so fitting the add and dom GRMs together splits the genetic variance into the classical additive and dominance variances.

The scaling is set with --code-type:

Both additive scalings are special cases of --alpha, which assumes that the variance of a SNP effect is proportional to \((2pq)^\alpha\). With \(v_i = 2p_i q_i\), the additive GRM is

so that the expected diagonal is 1 (for --ffile input see “Inputs” below). --alpha -1 gives --code-type 2 and --alpha 0 gives --code-type 1, with identical output files. Under \(\alpha = -1\) every SNP explains the same share of the genetic variance; values above \(-1\) shift the variance toward common SNPs and values below \(-1\) toward rare SNPs. Speed et al. (2017, Nature Genetics 49:986–992) recommended \(\alpha = -0.25\) for human complex traits. To choose \(\alpha\) for your data, build GRMs over a grid of values (for example −1, −0.75, −0.5, −0.25, 0) and compare the REML log-likelihoods. --alpha cannot be combined with --code-type and applies only to add. With missing genotypes, each element is rescaled by \(M/N_{jk}\) as for the other scalings (see “Quality control and missing genotypes” below).

For the dominance GRM the two --code-type choices apply, using the variance of the dominance coding instead of \(2pq\). An epistatic GRM (for example additive-by-additive) is built from existing GRMs with --process-grm --make-epis (see Section 3).

Inputs #

Provide exactly one of:

With --ffile:

• <prefix>.fmat has one row per feature and one whitespace-separated column per sample, in .fam order. It has no header and no row names.
• <prefix>.bim needs one line per feature; only the line count is used.
• Each feature is centred by its mean over the observed values. With the default --alpha -1 (--code-type 2) it is also scaled by its standard deviation, so the GRM does not depend on the units of the features. Negative values are fine.
• --code-type 1 and --alpha use the observed variance of feature \(i\) as \(v_i\) in place of \(2p_i q_i\). Multiplying one feature by a constant \(c\) multiplies its contribution to the GRM by \(c^{2(1+\alpha)}\), so only \(\alpha = -1\) is independent of the units. Use other values only when all features are on a common scale (for example log-expression from one assay). \(\alpha\) then sets how strongly features with a large variance dominate the GRM; unlike for SNPs, it has no population-genetic interpretation. With complete data the mean diagonal is exactly 1 for every \(\alpha\).
• Missing values are written as any of the --missing-geno tokens (default NA NaN nan N/A ...).
• --maf does not apply, and only --grm-type add is supported.

Selecting individuals and SNPs #

A line with one column is matched to the IID; a line with two or more columns is matched to FID and IID, as in PLINK. A leading chr in chromosome names is ignored. The options can be combined: individuals are kept and then removed, and SNPs are extracted, excluded and then restricted by chromosome. The output keeps the .fam order of the individuals. Allele frequencies and missing rates are computed in the selected individuals, so the result is identical to running on a file subset with plink --keep-allele-order --make-bed. The SNP options require --bfile.

A leave-one-chromosome-out (LOCO) GRM for chromosome c is built with --not-chr c, and the GRM of one chromosome with --chr c.

Allele frequencies #

By default p is estimated from the selected individuals. --read-freq <file> uses given A1 frequencies instead, for example from a reference population, so that GRMs built on different subsets of individuals are on the same scale. The file needs a header line with the columns SNP, A1 and FREQ (other columns are ignored), which is the format written by --write-freq:

CHR	SNP	A1	A2	FREQ	NOBS
0	rs1155740	T	C	0.27383863	409

FREQ is the frequency of the allele in the A1 column. If that allele is the .bim A2, the frequency is converted to 1 − FREQ. SNPs that are absent from the file, or whose A1 matches neither .bim allele, are excluded and counted in the log. The frequencies are used for centring, scaling and the --maf filter; the missing-rate filter still uses the analysed individuals. With the same frequencies, the GRM is identical to GCTA --update-freq.

--write-freq writes <out>.freq with one line per selected SNP: the A1 frequency used (estimated or read) and NOBS, the number of individuals with an observed genotype. With --npart, only part 1 writes the file. Both options require --bfile.

Quality control and missing genotypes #

A SNP (or feature) is excluded when any of the following holds:

• its minor allele frequency \(\min(p, 1 - p)\) is below --maf (BED input only)
• its missing rate is above --missing-rate
• it has zero variance

Missing genotypes are handled as in GCTA. A missing call contributes nothing to the cross-product, and each pair of samples is scaled by N_jk, the number of SNPs observed in both samples. These counts are written to .grm.N.bin.

Note. Because each element has its own denominator, a GRM built from data with missing genotypes is not guaranteed to be positive semi-definite (GCTA has the same property). With complete data it always is. If you need to invert such a GRM, use a larger --val with --process-grm --make-inv.

Options #

The full genotype matrix is never held in memory: the SNPs are read in --npart-snp chunks. Memory is dominated by the GRM block (about \(8n^2\) bytes for \(n\) individuals) plus one chunk (about \(8nm/S\) bytes for \(m\) SNPs and \(S\) = --npart-snp chunks), so raise --npart-snp when \(m\) is large. For very large samples, split the lower triangle into row blocks with --npart and merge them afterwards:

for i in 1 2 3 4; do
  fastgma --make-grm --bfile test --npart 4 $i --out output/test
done
fastgma --process-grm --merge --grm output/test --npart 4

Each job writes output/test.4_<i>.grm.* and its own log output/test.4_<i>.log. The merge step concatenates the .grm.id, .grm.bin and .grm.N.bin blocks into output/test.grm.* and copies .grm.N from the first block. The result is identical to a single full run.

2. Output files #

These are exactly the files GCTA writes (plus the text .grm.N).

Every module that reads a GRM accepts both value widths of .grm.bin (float32 as written now, or float64 from older fastGMA versions; the width is detected from the file size) and a .grm.id with one column (IID) or two columns (FID, IID). Samples are matched by IID (the second column, or the only column).

Each run also writes a log. It is appended to <out>.log, so several steps that share one --out prefix collect in one file. A partitioned run (--npart n i) logs to <out>.<n>_<i>.log. Set the environment variable FASTGMA_LOG_LEVEL=debug for more detail.

Reading a GRM in Python:

import numpy as np

ids = [line.split()[-1] for line in open("output/test.grm.id") if line.strip()]
n = len(ids)
tri = np.fromfile("output/test.grm.bin", dtype=np.float32)
assert tri.size == n * (n + 1) // 2

A = np.zeros((n, n))
A[np.tril_indices(n)] = tri
A = A + A.T - np.diag(np.diag(A))

3. Processing an existing GRM (--process-grm) #

Pass --process-grm together with exactly one operation:

Output formats for --out-fmt:

Notes:

• --merge. The blocks are read from <grm>.<n>_<i>.grm.* for i = 1..n and the merged GRM is written to <grm>.grm.*. .grm.N.bin is merged when the blocks have it. --out only names the log file (without it, the log goes to the console only).
• --pca. The files follow the GCTA layout. <out>.eigenval has one eigenvalue per line, largest first. <out>.eigenvec has one line per sample in .grm.id order with the columns FID IID PC1 PC2 ... (tab-separated, no header line). With a one-column .grm.id, the IID is written in both ID columns. --num-pc 20 (default) computes only the 20 leading eigenpairs, which is faster than a full decomposition; --num-pc 0 writes all of them.
• --make-epis. The two GRMs must have the same .grm.id order. The output is float32 <out>.grm.bin plus <out>.grm.id; .grm.N is the concatenation of the two input .grm.N files when both exist. No .grm.N.bin is written.
• --group. The output goes next to the input GRM, not to --out. Besides the group tables, it writes a sparse GRM (.grm.sp.bin) that keeps only within-group pairs.
• --make-inv. The inverse is computed with a LAPACK Cholesky factorization (dpotrf, dpotri). A GRM centred by sample allele frequencies has rank at most \(n - 1\), so it cannot be inverted without a diagonal value. That is why --val defaults to 0.001 here, while --make-grm itself adds nothing. With missing genotypes the GRM can have slightly negative eigenvalues (Section 1). If the matrix is not positive definite, --make-inv stops, reports the smallest eigenvalue and asks to rerun with --val above a given value. With --out-fmt 3 the inverse is stored as float32, like every .grm.bin.

4. Compatibility with GCTA #

With the same SNPs, the values agree with GCTA to float32 precision (largest difference about 1e-7), with or without missing genotypes.

• GCTA GRMs in fastGMA. Pass the GCTA prefix directly, for example --grm gcta_out.
• fastGMA GRMs in GCTA. --reml, --mlma, --pca, --grm-cutoff, --make-bK-sparse and --mgrm all read the files directly.
• Files from older fastGMA versions (float64 .grm.bin, one-column .grm.id) are still read by every module (see Section 2).

Differences to keep in mind when comparing the two tools:

To reproduce a GCTA GRM exactly, use --maf 0 (or the same --maf in both tools), --missing-rate 1, and only autosomal SNPs (--chr 1-22 for human data).

LDAK. With the same SNPs and complete genotypes, --alpha \(\alpha\) gives the kinship matrix of LDAK --calc-kins-direct --ignore-weights YES --power \(\alpha\) up to a constant factor: LDAK divides by the observed mean diagonal, so its matrix has mean diagonal exactly 1, whereas fastGMA (like GCTA) divides by the expected value \(\sum_i v_i^{1+\alpha}\). After dividing the fastGMA matrix by its mean diagonal, the two agree to float32 precision. With missing genotypes they also differ in the treatment of missing calls (LDAK uses one denominator for all pairs instead of N_jk).

5. Examples #

Additive GRM for uvLMM #

fastgma --make-grm --bfile test --out output/test
fastgma --uvlmm --grm output/test --bfile test \
  --data test_simu_pheno.txt --trait pheno --out output/test_uvlmm

Leave-one-chromosome-out GRMs #

for c in $(seq 1 22); do
  fastgma --make-grm --bfile test --not-chr $c --out output/test_loco$c
done

output/test_loco<c> is used when testing the SNPs of chromosome c.

Choosing α #

Build a GRM for each value on a grid and compare the REML log-likelihoods of the null model. Without --bfile, --uvlmm fits only the null model:

for a in -1 -0.75 -0.5 -0.25 0; do
  fastgma --make-grm --bfile test --alpha $a --out output/alpha_$a
  fastgma --uvlmm --grm output/alpha_$a --data test_simu_pheno.txt --trait pheno --out output/alpha_$a
done
grep "logL (REML" output/alpha_*.log

All models have the same fixed effects and the same number of variance components, so the log-likelihoods can be compared directly; the highest one gives the best-supported \(\alpha\). For the example data:

With 427 individuals the log-likelihood changes by only 0.04 over the grid, so these data cannot distinguish the values of \(\alpha\); estimating \(\alpha\) needs large samples.

Additive, dominance and epistatic GRMs for a multi-GRM model #

fastgma --make-grm --grm-type add --bfile test --out output/add
fastgma --make-grm --grm-type dom --bfile test --out output/dom
fastgma --process-grm --make-epis --prefix1 output/add --prefix2 output/add --out output/axa

fastgma --uvlmm --model multi-grm --mgrm output/add output/dom output/axa \
  --bfile test --data test_simu_pheno.txt --trait pheno --out output/test_mgrm

GRM from gene expression #

# expr.fmat: genes x samples; expr.fam: samples; expr.bim: one line per gene
fastgma --make-grm --ffile expr --missing-rate 0.2 --out output/expr

PCA, related-sample pruning and inverse #

fastgma --process-grm --pca --grm output/test --out output/test             # top 20 PCs
fastgma --process-grm --pca --num-pc 0 --grm output/test --out output/test  # all PCs
fastgma --process-grm --out-close --grm output/test --cut-value 0.125 --out output/test
fastgma --process-grm --make-inv --grm output/test --out-fmt 3 --out output/test.inv

If the genotypes contain missing calls, --make-inv may stop with “not positive definite”. The message gives the smallest eigenvalue and the --val needed; rerun with a larger value, for example --val 0.02.

Use the GRM in GCTA #

gcta64 --grm output/test --pheno test.phen --reml --out output/test_gcta
gcta64 --grm output/test --make-bK-sparse 0.05 --out output/test_sp