Genotype Data Formatting

Genotype Data Formatting#

This module implements a collection of workflows used to format genotype data.

Overview#

Nothing here changes the genotypes; it changes how they are packaged. Tools in the protocol disagree about format - some want PLINK, some want VCF - and the association and fine-mapping steps run per region or per chromosome, so the genotype data has to be split the same way to be processed in parallel. These workflows do those conversions and splits, plus the LD matrix computation per region that summary-statistic fine-mapping needs.

When to run it. Run the relevant workflow after genotype quality control and before association scanning, LD calculation, or fine-mapping that requires the corresponding genotype layout.

This module provides workflows for the formatting of genotype files. This includes the conversion between VCF and PLINK formats, the splitting of data (by specified input, chromosomes or genes) and the merging of data (by specified input, or by chromosomes).

When to run it. After VCF_QC, before PCA and association testing. Which workflows you need depends on the tools you plan to run, not on your data.

Input#

  • --genoFile: the genotype data to convert. Its meaning follows the workflow – one or more PLINK prefixes for plink_to_vcf and the merge_* workflows, one or more VCFs for vcf_to_plink.

  • --name (required): output prefix for this run, e.g. protocol_example.genotype.merged.

  • --cwd (default output): work directory for the analysis.

  • --modular-script-dir (default code/script): directory holding the modular wrapper scripts.

vcf_to_plink – VCF to PLINK:

  • --add-chr / --no-add-chr (default on): prefix contig names with chr.

  • --remove-duplicates / --no-remove-duplicates (default off): drop duplicated variant IDs.

  • --keep-dosage / --no-keep-dosage (default off): retain dosages rather than hard calls.

  • --remove-samples / --keep-samples (both default .): optional sample lists.

genotype_by_region – split PLINK genotypes into regions:

  • --region-list (required): regions to cut the genotypes into.

  • --window (default 0): flanking bp added either side of each region.

ld_by_region_plink – LD matrices per region:

  • --region-list (required): regions to compute LD for.

  • --float-type (default 16): float width of the stored matrix.

genotype_by_chrom – split by chromosome:

  • --chrom (required): the chromosomes to emit, e.g. --chrom 21 22.

write_data_list – write a manifest of produced files:

  • --out (required) and --ext (required): manifest path and the file extension to collect.

  • --data-files: the files to list, when they are not discovered from the previous step.

merge_plink – merge PLINK filesets:

  • --keep-samples (default .): optional sample list.

  • --extra-plink-opts: extra options passed straight to PLINK.

Cluster resources: --numThreads (default 20), --job-size (default 1), --walltime (default 5h), --mem (default 16G).

Output#

  • <name>.vcf.gz (+ .tbi) – from plink_to_vcf: the genotypes as a bgzipped, tabix-indexed VCF.

  • <name>.bed / .bim / .fam, or .pgen / .pvar / .psam – from vcf_to_plink, and from merge_plink for the merged fileset: a PLINK1 or PLINK2 fileset.

  • <name>.<region>.bed / .bim / .fam – from genotype_by_region: one fileset per region in --region-list, widened by --window.

  • <name>.<region>.ld.rds – from ld_by_region_plink: the LD matrix for a region, stored at the width --float-type sets.

  • <name>.<chrom>.bed / .bim / .fam – from genotype_by_chrom: one fileset per --chrom.

  • the manifest named by --out, a plain .txt or .list file of the produced paths – from write_data_list. Downstream modules read this rather than globbing.

  • a .stdout / .stderr log and the generated script beside each target.

All paths are relative to --cwd, and <name> is whatever --name was set to.

A .pvar file carries meta_file headers ahead of the variant table, so it must be read with a parser that skips them rather than as a plain TSV.

The merged fileset from Route A, output/genotype_formatting/plink/protocol_example.genotype.merged, holds 371,802 variants across 60 samples. First rows of the .bim:

1	chr1:113886_CTCTT_C	0	113886	C	CTCTT
1	chr1:191223_C_*	0	191223	*	C

and of the .fam, which carries no phenotype (-9):

0 SAMPLE_001 0 0 0 -9
0 SAMPLE_002 0 0 0 -9

Route A step 3 writes one fileset per chromosome plus its manifest:

protocol_example.genotype.merged.10.bed
protocol_example.genotype.merged.10.bim
protocol_example.genotype.merged.10.fam
protocol_example.genotype_by_chrom_files.txt

Routes B and C have no committed example output. Route B (plink_to_vcf, merge_vcf) and the manifest step are runnable but have not been executed here; Route C additionally needs a region list, which the protocol does not yet ship.

Minimal Working Example#

Split by region and compute LD#

Both steps read the merged fileset from Route A step 2, so run that first.

Step 2. LD matrices per region#

ld_by_region_plink computes an LD matrix per region and writes a .ld.list index alongside. The protocol does not yet ship an example region list, so this command cannot be run against the toy data as it stands; input/genotype/protocol_example.genotype.region_list.txt is the path it expects.

Timing: TBD (on toy dataset)

sos run pipeline/genotype_formatting.ipynb ld_by_region_plink \
    --genoFile output/genotype_formatting/plink/protocol_example.genotype.merged \
    --region-list input/genotype/protocol_example.genotype.region_list.txt \
    --float-type 16 \
    --cwd output/genotype_formatting/ld \
    --name protocol_example.genotype

Building a file manifest by hand#

write_data_list collects paths into a manifest that downstream modules read instead of globbing. The _1 conversion steps emit their own *_files.txt automatically; this workflow is for building one by hand.

Timing: TBD (on toy dataset)

sos run pipeline/genotype_formatting.ipynb write_data_list \
    --data-files `ls output/genotype_formatting/plink/protocol_example.genotype.chr*.bed` \
    --out output/genotype_formatting/plink/protocol_example.genotype.plink_files.txt \
    --ext bed \
    --cwd output/genotype_formatting/plink \
    --name protocol_example.genotype

Command Interface#

sos run pipeline/genotype_formatting.ipynb -h

usage: sos run pipeline/genotype_formatting.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:
  plink_to_vcf
  vcf_to_plink
  genotype_by_region
  ld_by_region_plink
  genotype_by_chrom
  write_data_list
  merge_plink
  merge_vcf

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --cwd output (as path)
                        Work directory & output directory
  --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

                        the path to a bed file or VCF file, a vector of bed
                        files or VCF files, or a text file listing the bed files
                        or VCF files to process
  --name VAL (as str, required)

Sections
  plink_to_vcf_1:
  vcf_to_plink:
    Workflow Options:
      --[no-]remove-duplicates (default to False)
      --[no-]add-chr (default to True)
      --[no-]keep-dosage (default to False)
      --remove-samples . (as path)
                        The path to the file that contains the list of samples
                        to remove (format FID, IID)
      --keep-samples . (as path)
                        The path to the file that contains the list of samples
                        to keep (format FID, IID)
  genotype_by_region_1:
    Workflow Options:
      --window 0 (as int)
                        cis window size
      --region-list VAL (as path, required)
                        Region definition
  ld_by_region_plink_1:
    Workflow Options:
      --region-list VAL (as path, required)
                        Region definition
      --float-type 16 (as int)
  genotype_by_chrom_1:
    Workflow Options:
      --chrom VAL VAL ... (as type, required)
  genotype_by_chrom_2:
  plink_to_vcf_2:
  genotype_by_region_2:
  ld_by_region_*_2:
    Workflow Options:
      --region-list VAL (as path, required)
  write_data_list:
    Workflow Options:
      --out VAL (as path, required)
      --ext VAL (as str, required)
      --data-files  paths

  merge_plink:
    Workflow Options:
      --keep-samples . (as path)
                        The path to the file that contains the list of samples
                        to keep (format FID, IID)
      --extra-plink-opts  (as list)
  merge_vcf:

Workflow implementation#

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
# Work directory & output directory
parameter: cwd = path("output")
import re
# 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
# the path to a bed file or VCF file, a vector of bed files or VCF files, or a text file listing the bed files or VCF files to process
parameter: genoFile = paths
parameter: name = str
# use this function to edit memory string for PLINK input
from sos.utils import expand_size
cwd = f"{cwd:a}"

import os
def get_genotype_file(geno_file_paths):
    def valid_geno_file(x):
        suffixes = path(x).suffixes
        if suffixes[-1] in ['.bed', '.pgen']:
            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 paths(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

# Determine if the file is in PLINK1 (BED/BIM/FAM) or PLINK2 (PGEN/PVAR/PSAM) format
def determine_plink_format(file_path):
    """
    Determine the PLINK file format based on file extensions and companion files.
    
    Args:
        file_path (str or Path): Path to the input file
    
    Returns:
        str: 'plink1' or 'plink2'
    """
    # Convert to string if it's a Path object
    file_path = str(file_path)
    
    # Check direct file extensions
    if file_path.endswith('.bed'):
        return 'plink1'
    elif file_path.endswith('.pgen'):
        return 'plink2'
    
    # If the file doesn't have a standard extension, try to infer format
    try:
        # Remove the file extension if present
        base_path = file_path.rsplit('.', 1)[0] if '.' in file_path else file_path
        
        # Check for PLINK1 companion files
        plink1_companion_files = [
            f"{base_path}.bim",
            f"{base_path}.fam"
        ]
        
        # Check for PLINK2 companion files
        plink2_companion_files = [
            f"{base_path}.pvar",
            f"{base_path}.psam"
        ]
        
        # Check PLINK1 format
        if all(os.path.exists(f) for f in plink1_companion_files):
            return 'plink1'
        
        # Check PLINK2 format
        if all(os.path.exists(f) for f in plink2_companion_files):
            return 'plink2'
    
    except Exception as e:
        print(f"Error determining PLINK format: {e}")
    
    # Default to PLINK1 if can't determine
    return 'plink1'

# Get the appropriate PLINK command prefix
def get_plink_command_prefix(file_path):
    format_type = determine_plink_format(file_path)
    if format_type == 'plink1':
        return "--bfile"
    else:  # plink2
        return "--pfile"

# Get output file extension based on format
def get_output_extension(format_type):
    if format_type == 'plink1':
        return "bed"
    else:  # plink2
        return "pgen"       
     
# Choose the make-bed or make-pgen command based on desired output format
def get_make_command(output_format):
    if output_format == 'plink1':
        return '--make-bed'
    else:  # plink2
        return '--make-pgen'
genoFile = get_genotype_file(genoFile)

Compute LD matrices for given input region#

ldstore2 based implementation#

Not implemented. No ldstore2 section exists in this notebook and the tool is not in the container, so nothing here is runnable; ld_by_region_plink above is the only LD workflow. ldstore2 would suit larger cohorts such as UK Biobank, which is not a constraint in the FunGen-xQTL project.

For whoever picks this up, the draft lives outside this repository: LDstore workflow and genotype preparation. It needs a semicolon-separated master file with one dataset per line, Z-file IDs formatted as rsid:chrom:pos:a1:a2 with zero-padded chromosomes, plus a sample list.

Merge VCF files#

[merge_vcf]
skip_if(len(genoFile) == 1)
# File prefix for the analysis output
parameter: name = str
input: genoFile, group_by = 'all'
output:  f"{cwd}/{name}.vcf.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:n}.stderr', stdout = f'{_output:n}.stdout'
    bash ${modular_script_dir}/data_preprocessing/genotype/genotype_formatting.sh merge_vcf \
        --genoFile "${_input}" \
        --output "${_output}"

bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    Rscript ${modular_script_dir}/data_preprocessing/genotype/genotype_formatting.R \
        --step vcf_gz_summary \
        --data-files "${_output}"