Genotype Data Formatting

Genotype Data Formatting#

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

Description#

We include steps 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).

Input Files#

This minimal working example starts from the per-chromosome toy genotype VCFs under input/genotype/:

File

Description

input/genotype/protocol_example.genotype.chr{1..22}.vcf.gz

22 per-chromosome toy genotype VCFs (60 samples)

The pipeline also accepts a single whole-genome VCF or PLINK bed/bim/fam bundle, a list of VCF/PLINK files, or a two-column (chrom, file) list.

Output Files#

File

Description

output/genotype_formatting/plink/protocol_example.genotype.chr{1..22}.{bed,bim,fam}

Per-chromosome PLINK bundles (Step 1)

output/genotype_formatting/plink/protocol_example.genotype.merged.{bed,bim,fam}

Single merged genome-wide PLINK bundle (Step 2)

output/genotype_formatting/genotype_by_chrom/protocol_example.genotype.merged.{1..22}.{bed,bim,fam}

Merged bundle partitioned back by chromosome (Step 3)

The merged PLINK bundle is the input expected by downstream genotype workers (GWAS_QC, GRM, PCA).

Minimal Working Example#

The data and singularity container bioinfo.sif can be downloaded from Synapse.

Step 3. Partition Genotype Data by Chromosome#

Split the merged genome-wide bundle back into per-chromosome PLINK files. Timing: <1 min on the toy data.

sos run pipeline/genotype_formatting.ipynb genotype_by_chrom \
    --genoFile output/genotype_formatting/plink/protocol_example.genotype.merged.bed \
    --cwd output/genotype_formatting/genotype_by_chrom \
    --chrom `cut -f 1 output/genotype_formatting/plink/protocol_example.genotype.merged.bim | uniq | sed "s/chr//g"` \
    -j 4

Command Interface#

sos run genotype_formatting.ipynb -h

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
# Work directory & output directory
parameter: cwd = path("output")
# The filename name for containers
parameter: container = ''
import re
parameter: entrypoint= ""
# 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 = ""
# 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#

This is good for larger sample sizes such as data from UK Biobank although we are not facing this challenge in the FunGen-xQTL project.

FIXME: we need to build ldstore2 into a container image. According to Diana it should be

pip3 install https://files.pythonhosted.org/packages/a8/fd/f98ab7dea176f42cb61b80450b795ef19b329e8eb715b87b0d13c2a0854d/ldstore-0.1.9.tar.gz

FIXME: Diana, what’s the input for this workflow?

Create master file for ldstore2#

The master file is a semicolon-separated text file and contains no space. It contains the following mandatory column names and one dataset per line:

FIXME: Diana, this documentation is not clearly written. I cannot understand it. What are the mandatory column names? What does it mean by one data-set per line?

  • For the Z file, the format should be rsid:chrom:pos:a1:a2. Formatting for chromosome should be 01,02,03 etc

  • List of samples

The LDstore draft is currently availale here with the code to prepare for the genotypic input here. A minimal working example can be found [here]

Pls note that .pvar file have meta_file headers before the variant info dataframe#

sos run /home/rl3328/xqtl-protocol/code/data_preprocessing/genotype/genotype_formatting.ipynb genotype_by_chrom
–genoFile /home/rl3328/motor_data/geno_data/final/LD_pruning/chrAll_GRCh38_biallelic.motor_eQTL.plink_qc.pgen
–cwd /home/rl3328/motor_data/geno_data/final/genotype_by_chrom
–chrom grep -v '^#' /home/rl3328/motor_data/geno_data/final/LD_pruning/chrAll_GRCh38_biallelic.motor_eQTL.plink_qc.pvar | cut -f 1 | uniq | sed "s/chr//g"

[genotype_by_chrom_1]
stop_if(len(paths(genoFile))>1, msg = "This workflow expects one input genotype file.")
parameter: chrom = list
chrom = list(set(chrom))
input: genoFile, for_each = "chrom"
plink_command = get_plink_command_prefix(_input)
output_format = determine_plink_format(_input)
output_ext = get_output_extension(output_format)
make_command = get_make_command(output_format)
output: f'{cwd}/{_input:bn}.{_chrom}.{output_ext}'
# look up for genotype file
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', container = container, entrypoint = entrypoint
    bash ${modular_script_dir}/data_preprocessing/genotype/genotype_formatting.sh genotype_by_chrom \
        --cwd "${cwd}" \
        --genoFile "${_input}" \
        --name "${name}" \
        --output "${_output}" \
        --plink-command "${plink_command}" \
        --make-command "${make_command}" \
        --output-format "${output_format}" \
        --chrom ${_chrom} \
        --mem "${mem}" \
        --numThreads ${numThreads}

[genotype_by_chrom_2]
input: group_by = "all"
output_format = determine_plink_format(_input)
output_ext = get_output_extension(output_format)
output: f'{_input[0]:nn}.{step_name[:-2]}_files.txt'
sos_run("write_data_list", data_files = _input, out = _output, ext = f"{output_ext}")

[plink_to_vcf_2]
input: group_by = "all"
output: f'{_input[0]:nnn}.{step_name[:-2]}_files.txt'
sos_run("write_data_list", data_files = _input, out = _output, ext = "vcf.gz")

[genotype_by_region_2]
input: group_by = "all"
output: f'{_input[0]:nn}.{step_name[:-2]}_files.txt'
sos_run("write_data_list", data_files = _input, out = _output, ext = "bed")

[ld_by_region_*_2]
parameter: region_list = path
input: group_by = "all"
output: f'{cwd}/{region_list:bn}_LD/{genoFile:bn}.ld.list'
sos_run("write_data_list", data_files = _input, out = _output, ext = "npy.gz")

[write_data_list]
parameter: out = path
parameter: ext = str
parameter: data_files = paths
input: data_files
output: out
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', container = container, entrypoint = entrypoint
    python3 ${modular_script_dir}/data_preprocessing/genotype/genotype_formatting.py \
        --step write_data_list \
        --data-files "${"::".join([str(x) for x in _input])}" \
        --ext "${ext}" \
        --output "${_output}" \
        --numThreads ${numThreads}

Split VCF by Chromosome#

FIXME: add this as needed

Merge VCF files#

Anticipated Results#

The pipeline produces output files in the output/ subdirectory named after the workflow step. Verify success by checking that output files exist and are non-empty. See the Output section above for the expected file names and formats.

[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: container=container, expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', entrypoint=entrypoint
    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', container = container, entrypoint=entrypoint
    python3 ${modular_script_dir}/data_preprocessing/genotype/genotype_formatting.py \
        --step vcf_gz_summary \
        --data-files "${_output}"