Genomic Relationship Matrices

Genomic Relationship Matrices#

This workflow generates genomic relationship matrices (GRM) under the leave-one-chromosome-out (LOCO) theme.

Overview#

A genomic relationship matrix summarises how genetically similar every pair of samples is, and association models use it to absorb that relatedness rather than mistake it for signal. Computing it from all chromosomes at once creates a problem: the variant being tested contributes to the matrix that is supposed to control for background relatedness, which deflates the very association being measured.

The leave-one-chromosome-out construction avoids that. For each chromosome the GRM is built from every other chromosome, so the variants under test never appear in the matrix correcting them. GCTA does the computation, and a final step reformats the result for APEX.

Calling the grm workflow runs four numbered steps in sequence: write the leave-one-out file list per chromosome, compute the GRM with GCTA, format it for APEX, then collect the per-chromosome results into one list.

When to run it. Once per genotype set, before any association scan that takes a GRM.

Input#

--genoFile: a two-column, tab-separated list with one row per chromosome, giving the chromosome number and the path to that chromosome’s PLINK 1.0 .bed. Its .bim and .fam must sit beside each .bed. Common variants, MAF above 1%, are recommended. Example input/genotype/protocol_example.genotype.list, 3 rows:
```
20  output/plink/protocol_example.genotype.chr20.bed
21  output/plink/protocol_example.genotype.chr21.bed
22  output/plink/protocol_example.genotype.chr22.bed
```
--cwd: the directory outputs are written to.
--numThreads: threads handed to GCTA, 20 by default.
--job-size, --walltime and --mem: cluster job resources.

Output#

{cwd}/{name}.chr<N>.loco.txt - for chromosome N, the list of every other chromosome’s genotype file. This is what makes it leave-one-out. Example output/grm_uf/protocol_example.genotype.chr20.loco.txt, 2 lines:
```
output/plink/protocol_example.genotype.chr21
output/plink/protocol_example.genotype.chr22
```
{cwd}/{name}.chr<N>.loco.grm.gz - the GRM for chromosome N computed by GCTA from all other chromosomes, with .grm.id naming the samples and .log recording the run.
{cwd}/{name}.chr<N>.loco.apex.grm - the same matrix reformatted for APEX.

{cwd}/{name}.loco_grm_list.txt - the per-chromosome results collected into one table, #chr and the GRM directory. Example:

#chr	dir
chr20chr21chr22	output/grm_uf/protocol_example.genotype.chr20.apex.grm
chr20chr21chr22	output/grm_uf/protocol_example.genotype.chr21.apex.grm

Minimal Working Example#

Timing: TBD (on toy dataset)

sos run pipeline/GRM.ipynb grm \
    --cwd output/grm_uf \
    --genoFile input/genotype/protocol_example.genotype.list

Command Interface#

sos run GRM.ipynb -h

[91mERROR[0m: [91mNotebook JSON is invalid: %s[0m
usage: sos run code/SoS/data_preprocessing/genotype/GRM.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  grm

Global Workflow Options:
  --cwd output (as path)
                        Work directory & output directory
  --modular-script-dir code/script (as path)
                        The filename name for output data
  --job-size 1 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 5h
                        Wall clock time expected
  --mem 16G
                        Memory expected
  --numThreads 20 (as int)
                        Number of threads
  --genoFile  paths


Sections
  grm_1:                Generate LOCO file list
  grm_2:                Compute LOCO-GRM
  grm_3:                Format output for APEX
  grm_4:                Generate APEX input list

Workflow implementation#

[global]
# Work directory & output directory
parameter: cwd = path("output")
# The filename name for output data
parameter: modular_script_dir = path('code/script')
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "5h"
# Memory expected
parameter: mem = "16G"
# Number of threads
parameter: numThreads = 20
parameter: genoFile = paths
cwd = f"{cwd:a}"

## Code below are copied from genotype_formatting.ipynb. Change should be made to codes in that document and sync with this chunk.
from sos.utils import expand_size
import os

def get_genotype_file(geno_file_paths):
    #
    def valid_geno_file(x):
        suffixes = path(x).suffixes
        if suffixes[-1] == '.bed':
            return True
        if len(suffixes)>1 and ''.join(suffixes[-2:]) == ".vcf.gz":
            return True
        return False
    #
    def complete_geno_path(x, geno_file):
        if not valid_geno_file(x):
            raise ValueError(f"Genotype file {x} should be VCF (end with .vcf.gz) or PLINK bed file (end with .bed)")
        if not os.path.isfile(x):
            # relative path
            if not os.path.isfile(f'{geno_file:ad}/' + x):
                raise ValueError(f"Cannot find genotype file {x}")
            else:
                x = f'{geno_file:ad}/' + x
        return x
    # 
    def format_chrom(chrom):
        if chrom.startswith('chr'):
            chrom = chrom[3:]
        return chrom
    # Inputs are either VCF or bed, or a vector of them 
    if len(geno_file_paths) > 1:
        if all([valid_geno_file(x) for x in geno_file_paths]):
            return paths(geno_file_paths)
        else: 
            raise ValueError(f"Invalid input {geno_file_paths}")
    # Input is one genotype file or text list of genotype files
    geno_file = geno_file_paths[0]
    if valid_geno_file(geno_file):
        return path(geno_file)
    else: 
        units = [x.strip().split() for x in open(geno_file).readlines() if x.strip() and not x.strip().startswith('#')]
        if all([len(x) == 1 for x in units]):
            return paths([complete_geno_path(x[0], geno_file) for x in units])
        elif all([len(x) == 2 for x in units]):
            genos = dict([(format_chrom(x[0]), path(complete_geno_path(x[1], geno_file))) for x in units])
        else:
            raise ValueError(f"{geno_file} should contain one column of file names, or two columns of chrom number and corresponding file name")
        return genos
                        
genotypes = get_genotype_file(genoFile)

Generate LOCO file list#

# Generate LOCO file list
[grm_1]
chrom = list(genotypes.keys())
input: for_each = 'chrom'
output: f'{cwd}/{genotypes[_chrom]:bn}.loco.txt'
with open(_output, 'w') as f:
    f.write('\n'.join([str(genotypes[x].with_suffix('')) for x in genotypes if x != _chrom]))

Compute LOCO-GRM#

# Compute LOCO-GRM
[grm_2]
output: f'{_input:nn}.grm.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand= "$[ ]", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout'
   gcta \
   --mbfile $[_input] \
   --make-grm-gz \
   --out $[_output:nn]

Format output for APEX#

# Format output for APEX
[grm_3]
output: f'{_input:nn}.apex.grm'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand= "$[ ]", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout'
    Rscript $[modular_script_dir:a]/data_preprocessing/genotype/GRM.R --step format_apex \
        --grm $[_input] \
        --id $[_input:n].id \
        --output $[_output]

Generate APEX input list#

# Generate APEX input list
[grm_4]
input: group_by = "all"
output: f'{cwd}/{path(list(genotypes.values())[1]):bnn}.loco_grm_list.txt'
bash: expand= "$[ ]", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout'
    Rscript $[modular_script_dir:a]/data_preprocessing/genotype/GRM.R --step apex_list \
        --input $[_input:r,] \
        --chroms $[" ".join(genotypes.keys())] \
        --output $[_output]