Genotype PLINK File Quality Control

Runs kinship estimation, variant- and sample-level filtering, and LD pruning over a merged PLINK genotype set, producing the QC-ed genotypes and the pruned variant list used for PCA.

Overview

This module is the standard quality-control pass over a merged PLINK genotype set. It estimates kinship to identify related individuals, filters variants and samples on allele frequency, missingness and Hardy-Weinberg equilibrium, and LD-prunes what survives into the variant set used for PCA. Input may be PLINK1 (bed/bim/fam) or PLINK2 (pgen/pvar/psam); VCF has to be converted first by genotype_formatting.

Four workflows are available, and how they are combined depends on whether the cohort contains relatives.

• king estimates kinship with KING, flags pairs above --kinship, and splits the cohort into an *.unrelated and a *.related PLINK set. When no pair passes the threshold it writes an empty related-ID list and stops without splitting.

• qc applies the variant- and sample-level filters and then LD-prunes, writing both the filtered set and the .prune.in variant list.

• qc_no_prune is that same filtering step without the pruning stage, which is what to run when the pruned variant list already exists and only has to be extracted.

• genotype_phenotype_sample_overlap intersects the genotype sample list with the samples in a molecular phenotype file. A two-column --sample-participant-lookup (genotype_id, sample_id) translates between naming schemes; without one, the names are assumed to match already.

The usual sequence is qc_no_prune for basic filtering, genotype_phenotype_sample_overlap to restrict to samples that have phenotype data, king to detect relatives, and finally qc on the unrelated set to produce the pruned variants for PCA. If KING finds no relatives there is nothing to split, and qc runs on the full filtered set with --keep-samples instead.

The defaults follow Table 1 of this GWAS quality-control tutorial and are deliberately permissive, because the appropriate thresholds depend on the analysis that follows:

• kinship coefficient 0.0625, which is third-degree relatives and closer

• MAF and MAC both 0, so common and rare variants are kept

• variant-level and sample-level missingness both 0.1

• HWE 1e-15, which is very lenient

• LD pruning over a 50-variant window, shifting 10 variants at a time, at r2 0.1

For PCA a MAF cutoff of 0.01 is the usual choice, since common variants are the appropriate basis; for single-variant association a MAC floor of 5 is typical.

One constraint is structural: because both sample- and variant-level QC are applied together, all samples and chromosomes have to be merged into a single file first. That has been run at the scale of 200K exomes and 15 million variants on one merged PLINK set.

When to run it. After genotype_formatting has merged the per-chromosome genotypes, and before PCA and any association scan. The sample-overlap workflow additionally needs the molecular phenotype file to exist.

Input

• --genoFile – the merged genotype set, PLINK1 bed/bim/fam or PLINK2 pgen/pvar/psam, e.g. output/genotype_formatting/plink/protocol_example.genotype.merged.bed. king takes the QC-ed .bed, while genotype_phenotype_sample_overlap takes the .fam instead. Samples are identified by the FID and IID columns:

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

• --phenoFile – the molecular phenotype whose samples are matched against the genotype, bed.gz or tsv, e.g. input/rnaseq/protocol_example.rnaseq.bed.gz. Sample IDs are the column names after the four positional columns:

  #chr	start	end	ID	SAMPLE_001	SAMPLE_002	SAMPLE_003	SAMPLE_004	SAMPLE_005	SAMPLE_006	SAMP
  

• --name – string identifying the run; it is inserted into every output file name.

• --cwd – output directory, default output.

Restricting samples and variants:

• --keep-samples / --remove-samples – FID/IID lists limiting or excluding samples; --keep-samples is how the QC step is confined to, say, one ancestry group, or to the samples that have phenotype data.

• --keep-variants / --exclude-variants – variant ID lists. Supplying --keep-variants also inserts .extracted into the output file name.

• --sample-participant-lookup – two-column table (genotype_id, sample_id) used when genotype and phenotype name their samples differently.

Filter thresholds:

• --kinship – kinship coefficient above which a pair counts as related, default 0.0625.

• --kin-maf – MAF floor applied before the KING estimate, default 0.01.

• --maf-filter / --maf-max-filter and --mac-filter / --mac-max-filter – allele frequency and count bounds; 0 disables a bound.

• --geno-filter – maximum missingness per variant, default 0.1.

• --mind-filter – maximum missingness per sample, default 0.1.

• --hwe-filter – Hardy-Weinberg p-value cutoff, default 1e-15.

• --rm-dups – drop duplicate variants.

• --treat-dosage-missing – handle dosage data.

• --meta-only – write only the SNP and sample lists rather than a PLINK binary set.

• --other-args – extra PLINK arguments passed through, such as snps_only or write-samples.

LD pruning:

• --window, --shift, --r2 – pruning window in variants (default 50), how far it moves each time (10), and the r2 ceiling (0.1). Setting --r2 0 skips pruning altogether.

• --bad-ld – pass PLINK --bad-ld, which suppresses the pruning diagnostics that otherwise abort on regions of extreme LD.

Runtime:

• --numThreads, --job-size, --walltime, --mem – threads (default 20) and cluster resources.

• --modular-script-dir – location of the shell and R drivers, default code/script.

Output

• <genotype>.<name>.plink_qc.{bed,bim,fam} – the filtered genotype set, written by qc_no_prune and by the filtering step of qc. .extracted is inserted when --keep-variants was supplied, and --meta-only produces a .snplist instead of a binary set.

• <genotype>.<name>.plink_qc.prune.{bed,bim,fam} and <genotype>.<name>.plink_qc.prune.in – the LD-pruned genotype set and the variant list that defines it. The .prune.in list is what PCA consumes:

  chr1:820919_CACTACCTGCTTGTCCAGCAGGTCCACCCTGTCTACACTACCTGCCTGCAAAGCAGATCCACCCTGTCTACACTACCTGGCTGG
  chr1:820929_TTGTCCAGCAGGTCCACCCTGTCTACACTACCTGCCTGCAAAGCAGATCCACCCTGTCTACACTACCTGGCTGGCCAGTAGATC
  chr1:820929_TTGTCCAGCAGGTCCACCCTGTCTACACTACCTGCCTGCAAAGCAGATCCACCCTGTCTACACTACCTGGCTGGCCAGTAGATC
  

• <genotype>.<name>.kin0 – the KING kinship table, one row per pair above the threshold, with the kinship coefficient in the last column:

  #FID1	IID1	FID2	IID2	NSNP	HETHET	IBS0	KINSHIP
  0	SAMPLE_060	0	SAMPLE_059	59818	0.108178	0	0.248649
  

• <genotype>.<name>.related_id – the individuals dropped to leave a maximal unrelated set; the file is created empty when no pair passes --kinship:

  0 SAMPLE_060
  

• <genotype>.<name>.unrelated.{bed,bim,fam} and <genotype>.<name>.related.{bed,bim,fam} – the split cohort, written only when relatives were found.

• <phenotype>.sample_genotypes.txt – FID/IID of the samples present in both genotype and phenotype, in the form PLINK --keep expects:

  0	SAMPLE_001
  0	SAMPLE_002
  0	SAMPLE_003
  

• <phenotype>.sample_overlap.txt – the same overlap as a genotype_id / sample_id table:

  genotype_id	sample_id
  SAMPLE_001	SAMPLE_001
  SAMPLE_002	SAMPLE_002
  

Each step also writes PLINK .log files and .stdout / .stderr beside its output. Example results for the toy data are under output/gwas_qc/.

Minimal Working Example

Genotype QC of the protocol example data

The chr1_chr6 set was merged from the chr1 and chr6 data with the merge_plink command in genotype formatting. The steps below run in order, each consuming what the earlier ones produced.

Step 1. Basic QC (rare and common variants)

Apply the variant- and sample-level filters: missingness, HWE and MAC.

Timing: <1 min (on toy dataset)

sos run pipeline/GWAS_QC.ipynb qc_no_prune \
    --cwd output/gwas_qc/plink \
    --genoFile output/genotype_formatting/plink/protocol_example.genotype.merged.bed \
    --name protocol_example \
    --geno-filter 0.1 \
    --mind-filter 0.1 \
    --hwe-filter 1e-08 \
    --mac-filter 0

Step 2. Sample match with phenotype

Find the samples shared between genotype and phenotype and write the overlapping sample lists.

Timing: <1 min (on toy dataset)

sos run pipeline/GWAS_QC.ipynb genotype_phenotype_sample_overlap \
    --cwd output/gwas_qc/genotype \
    --genoFile output/gwas_qc/plink/protocol_example.genotype.merged.plink_qc.fam \
    --phenoFile input/rnaseq/protocol_example.rnaseq.bed.gz \
    --name protocol_example

Step 3. Kinship QC

Estimate kinship with KING and split the samples into related and unrelated sets. In the toy data SAMPLE_059 and SAMPLE_060 are a parent-offspring pair, so KING does find relatives and writes both *.related.bed and *.unrelated.bed.

Timing: <2 min (on toy dataset)

sos run pipeline/GWAS_QC.ipynb king \
    --cwd output/gwas_qc/kinship \
    --genoFile output/gwas_qc/plink/protocol_example.genotype.merged.plink_qc.bed \
    --name protocol_example.king \
    --keep-samples output/gwas_qc/genotype/protocol_example.rnaseq.bed.sample_genotypes.txt

Step 4. Prepare unrelated individuals for PCA

Because Step 3 found related individuals, run qc on the KING *.unrelated.bed to produce the LD-pruned, unrelated genotype set used downstream for PCA.

Timing: <1 min (on toy dataset)

sos run pipeline/GWAS_QC.ipynb qc \
    --cwd output/gwas_qc/genotype \
    --genoFile output/gwas_qc/kinship/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.bed \
    --name protocol_example \
    --mac-filter 5

Alternative to Step 4: no related individuals

If KING reports no related individuals and writes no *.unrelated.bed, run qc on the QC-ed genotype with --keep-samples instead. Relatives are present in this toy data, so Step 4 is the path that applies here and the command below is shown only for reference.

Timing: <1 min (on toy dataset)

sos run pipeline/GWAS_QC.ipynb qc \
    --cwd output/gwas_qc/genotype \
    --genoFile output/gwas_qc/plink/protocol_example.genotype.merged.plink_qc.bed \
    --name protocol_example \
    --mac-filter 5

Step 5. Extract pruned variants for PCA

Extract the LD-pruned variants from Step 4 out of the full genotype set, applying only sample-level missingness, in preparation for PCA.

Timing: TBD (on toy dataset)

sos run pipeline/GWAS_QC.ipynb qc_no_prune \
    --cwd output/gwas_qc/cache \
    --genoFile output/genotype_formatting/plink/protocol_example.genotype.merged.bed \
    --geno-filter 0 \
    --mind-filter 0.1 \
    --maf-filter 0 \
    --keep-variants output/gwas_qc/genotype/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.plink_qc.prune.in \
    --name for_pca

Command Interface

sos run pipeline/GWAS_QC.ipynb -h

usage: sos run pipeline/GWAS_QC.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:
  king
  qc_no_prune
  qc
  genotype_phenotype_sample_overlap

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --cwd output (as path)
                        the output directory for generated files
  --name VAL (as str, required)
                        A string to identify your analysis run
  --genoFile  paths

                        PLINK binary files (either BED/BIM/FAM or PGEN/PVAR/PSAM
                        format)
  --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)
  --keep-variants . (as path)
                        The path to the file that contains the list of variants
                        to keep
  --exclude-variants . (as path)
                        The path to the file that contains the list of variants
                        to exclude
  --kinship 0.0625 (as float)
                        Kinship coefficient threshold for related individuals
                        (e.g first degree above 0.25, second degree above 0.125,
                        third degree above 0.0625)
  --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

Sections
  king_1:               Inference of relationships in the sample to identify
                        closely related individuals
    Workflow Options:
      --kin-maf 0.01 (as float)
                        PLINK binary file
  king_2:               Select a list of unrelated individual with an attempt to
                        maximize the unrelated individuals selected from the
                        data
  king_3:               Split genotype data into related and unrelated samples,
                        if related individuals are detected
  qc_no_prune, qc_1:    Filter SNPs and select individuals
    Workflow Options:
      --maf-filter 0.0 (as float)
                        minimum MAF filter to use. 0 means do not apply this
                        filter.
      --maf-max-filter 0.0 (as float)
                        maximum MAF filter to use. 0 means do not apply this
                        filter.
      --mac-filter 0.0 (as float)
                        minimum MAC filter to use. 0 means do not apply this
                        filter.
      --mac-max-filter 0.0 (as float)
                        maximum MAC filter to use. 0 means do not apply this
                        filter.
      --geno-filter 0.1 (as float)
                        Maximum missingess per-variant
      --mind-filter 0.1 (as float)
                        Maximum missingness per-sample
      --hwe-filter 1e-15 (as float)
                        HWE filter -- a very lenient one
      --other-args  (as list)
                        Other PLINK arguments e.g snps_only, write-samples, etc
      --[no-]meta-only (default to False)
                        Only output SNP and sample list, rather than the PLINK
                        binary format of subset data
      --[no-]rm-dups (default to False)
                        Remove duplicate variants
      --[no-]treat-dosage-missing (default to False)
                        Add option to process dosage
  qc_2:                 LD prunning and remove related individuals (both ind of
                        a pair) Plink2 has multi-threaded calculation for LD
                        prunning
    Workflow Options:
      --window 50 (as int)
                        Window size
      --shift 10 (as int)
                        Shift window every 10 snps
      --r2 0.1 (as float)
      --mac-filter 0.0 (as float)
      --other-args  (as list)
      --[no-]bad-ld (default to False)
                        Use PLINK --bad-ld flag (skip LD pruning diagnostics for
                        regions with extreme LD)
  genotype_phenotype_sample_overlap: This workflow extracts overlapping samples
                        for genotype data with phenotype data, and output the
                        filtered sample genotype list as well as sample
                        phenotype list
    Workflow Options:
      --phenoFile VAL (as path, required)
                        A phenotype file, can be bed.gz or tsv
      --sample-participant-lookup . (as path)
                        If this file is provided, a genotype/phenotype sample
                        name match will be performed It must contain two column
                        names: genotype_id, sample_id

Workflow implementation

The SoS workflow definitions below are unchanged from the original protocol.

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
# the output directory for generated files
parameter: cwd = path("output")
# A string to identify your analysis run
parameter: name = str
# PLINK binary files (either BED/BIM/FAM or PGEN/PVAR/PSAM format)
parameter: genoFile = paths
# The path to the file that contains the list of samples to remove (format FID, IID)
parameter: remove_samples = path('.')
# The path to the file that contains the list of samples to keep (format FID, IID)
parameter: keep_samples = path('.')
# The path to the file that contains the list of variants to keep
parameter: keep_variants = path('.')
# The path to the file that contains the list of variants to exclude
parameter: exclude_variants = path('.')
# Kinship coefficient threshold for related individuals
# (e.g first degree above 0.25, second degree above 0.125, third degree above 0.0625)
parameter: kinship = 0.0625
# 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
# use this function to edit memory string for PLINK input
from sos.utils import expand_size
cwd = path(f"{cwd:a}")

# 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 based on the input file format
def get_plink_command_prefix(file_path):
    format_type = determine_plink_format(file_path)
    if format_type == 'plink1':
        return "--bfile"
    else:  # plink2
        return "--pfile"
        
# Generate the appropriate file extension based on the requested format
def get_output_extension(output_format, is_prune=False):
    if output_format == 'plink1':
        return '.bed' if not is_prune else '.prune.bed'
    else:  # plink2
        return '.pgen' if not is_prune else '.prune.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'

def get_other_args_flags(other_args):
    import re
    if not other_args:
        return ""
    args = [other_args] if isinstance(other_args, str) else list(other_args)
    flags = []
    for arg in args:
        arg = str(arg)
        if not re.fullmatch(r"[A-Za-z0-9][A-Za-z0-9_.:-]*", arg):
            raise ValueError(f"Cannot safely pass PLINK other_args value: {arg}")
        flags.append(f"--other-arg {arg}")
    return " ".join(flags)

Kinship estimation and cohort split

king_1 runs KING to produce the kinship table, king_2 chooses which individuals to drop, and king_3 writes the related and unrelated PLINK sets.

# Inference of relationships in the sample to identify closely related individuals
[king_1]
# PLINK binary file
parameter: kin_maf = 0.01
input: genoFile
output: f'{cwd}/{_input:bn}{("."+name) if name else ""}.kin0'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
plink_command = get_plink_command_prefix(genoFile)
bash: expand= "${ }", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout'
    bash ${modular_script_dir}/data_preprocessing/genotype/GWAS_QC.sh king \
        --cwd "${cwd}" \
        --genoFile "${genoFile}" \
        --plink-command "${plink_command}" \
        --out-prefix "${_output:n}" \
        --keep-samples "${keep_samples}" \
        --remove-samples "${remove_samples}" \
        --name "${name}" \
        --kinship ${kinship} \
        --kin-maf ${kin_maf} \
        --numThreads ${numThreads}

# Select a list of unrelated individual with an attempt to maximize the unrelated individuals selected from the data 
[king_2: shared = "related_id" ]
related_id = [x.strip() for x in open(_input).readlines() if not x.startswith("#")]
output: f'{_input:n}.related_id'
with open(_output, 'a'):
    pass
done_if(len(related_id) == 0, msg = f"No related individuals detected from {_input}.")
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}/data_preprocessing/genotype/GWAS_QC.R \
        --step king_2 \
        --input "${_input}" \
        --output "${_output}" \
        --kinship ${kinship}

# Split genotype data into related and unrelated samples, if related individuals are detected
[king_3]
depends: sos_variable("related_id")
input: output_from(2), genoFile
output_format = determine_plink_format(_input[1])
output_ext = get_output_extension(output_format)
output: unrelated_bed = f'{cwd}/{_input[0]:bn}.unrelated{output_ext}',
        related_bed = f'{cwd}/{_input[0]:bn}.related{output_ext}'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    bash ${modular_script_dir}/data_preprocessing/genotype/GWAS_QC.sh king_split \
        --cwd "${cwd}" \
        --genoFile "${_input[1]}" \
        --remove-samples "${_input[0]}" \
        ${('--keep-samples "%s"' % keep_samples) if keep_samples.is_file() else ""} \
        --out-unrelated "${_output[0]:n}" \
        --out-related "${_output[1]:n}" \
        --related-output "${_output[1]}" \
        --numThreads ${numThreads}

Variant and sample filtering

qc_no_prune, which doubles as the first step of qc, applies the frequency, missingness and HWE filters; qc_2 performs the LD pruning.

# Filter SNPs and select individuals 
[qc_no_prune, qc_1 (basic QC filters)]
# minimum MAF filter to use. 0 means do not apply this filter.
parameter: maf_filter = 0.0
# maximum MAF filter to use. 0 means do not apply this filter.
parameter: maf_max_filter = 0.0
# minimum MAC filter to use. 0 means do not apply this filter.
parameter: mac_filter = 0.0
# maximum MAC filter to use. 0 means do not apply this filter.
parameter: mac_max_filter = 0.0 
# Maximum missingess per-variant
parameter: geno_filter = 0.1
# Maximum missingness per-sample
parameter: mind_filter = 0.1
# HWE filter -- a very lenient one
parameter: hwe_filter = 1e-15
# Other PLINK arguments e.g snps_only, write-samples, etc
parameter: other_args = []
# Only output SNP and sample list, rather than the PLINK binary format of subset data
parameter: meta_only = False
# Remove duplicate variants
parameter: rm_dups = False
# Add option to process dosage
parameter: treat_dosage_missing = False

fail_if(not (keep_samples.is_file() or keep_samples == path('.')), msg = f'Cannot find ``{keep_samples}``')
fail_if(not (keep_variants.is_file() or keep_variants == path('.')), msg = f'Cannot find ``{keep_variants}``')
fail_if(not (remove_samples.is_file() or remove_samples == path('.')), msg = f'Cannot find ``{remove_samples}``')

input: genoFile, group_by=1
plink_command = get_plink_command_prefix(_input)
output_format = determine_plink_format(_input)
make_command = get_make_command(output_format) if not meta_only else "--write-snplist --write-samples"
output_ext = get_output_extension(output_format) if not meta_only else ".snplist"
other_args_flags = get_other_args_flags(other_args)
output: f'{cwd}/{_input:bn}{("." + name) if name else ""}.plink_qc{".extracted" if keep_variants.is_file() else ""}{output_ext}'
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/GWAS_QC.sh qc_no_prune \
        --cwd "${cwd}" \
        --genoFile "${_input}" \
        --plink-command "${plink_command}" \
        --output-format "${output_format}" \
        --make-command "${make_command}" \
        --out-prefix "${_output:n}" \
        --name "${name}" \
        --mac-filter ${mac_filter} \
        --maf-filter ${maf_filter} \
        --maf-max-filter ${maf_max_filter} \
        --mac-max-filter ${mac_max_filter} \
        --geno-filter ${geno_filter} \
        --mind-filter ${mind_filter} \
        --hwe-filter ${hwe_filter} \
        ${"--keep-samples " + str(keep_samples) if keep_samples.is_file() else ""} \
        ${"--remove-samples " + str(remove_samples) if remove_samples.is_file() else ""} \
        ${"--exclude-variants " + str(exclude_variants) if exclude_variants.is_file() else ""} \
        ${"--keep-variants " + str(keep_variants) if keep_variants.is_file() else ""} \
        ${"--meta-only" if meta_only else ""} \
        ${"--rm-dups" if rm_dups else ""} \
        ${"--treat-dosage-missing" if treat_dosage_missing else ""} \
        ${other_args_flags} \
        --numThreads ${numThreads}

# LD prunning and remove related individuals (both ind of a pair)
# Plink2 has multi-threaded calculation for LD prunning
[qc_2 (LD pruning)]
# Window size
parameter: window = 50
# Shift window every 10 snps
parameter: shift = 10
parameter: r2 = 0.1
parameter: mac_filter = 0.0
parameter: other_args = []
# Use PLINK --bad-ld flag (skip LD pruning diagnostics for regions with extreme LD)
parameter: bad_ld = False
stop_if(r2==0)
plink_command = get_plink_command_prefix(_input)
output_format = determine_plink_format(_input)
make_command = get_make_command(output_format)
output_ext = get_output_extension(output_format, is_prune=True)
other_args_flags = get_other_args_flags(other_args)
output: bed=f'{cwd}/{_input:bn}{output_ext}', prune=f'{cwd}/{_input:bn}.prune.in'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    bash ${modular_script_dir}/data_preprocessing/genotype/GWAS_QC.sh qc \
        --cwd "${cwd}" \
        --genoFile "${_input}" \
        --plink-command "${plink_command}" \
        --output-format "${output_format}" \
        --make-command "${make_command}" \
        --prune-prefix "${_output['prune']:nn}" \
        --out-prefix "${_output['bed']:n}" \
        --mac-filter ${mac_filter} \
        --window ${window} \
        --shift ${shift} \
        --r2 ${r2} \
        ${"--bad-ld" if bad_ld else ""} \
        ${other_args_flags} \
        --numThreads ${numThreads}

Genotype-phenotype sample matching

An auxiliary step matching genotype to phenotype samples through the optional look-up table; if no table is given, or the file is not found, the names are assumed to match already.

# This workflow extracts overlapping samples for genotype data with phenotype data, and output the filtered sample genotype list as well as sample phenotype list
[genotype_phenotype_sample_overlap]
# A genotype fam file
parameter: genoFile = path
# A phenotype file, can be bed.gz or tsv
parameter: phenoFile = path
# If this file is provided, a genotype/phenotype sample name match will be performed
# It must contain two column names: genotype_id, sample_id
parameter: sample_participant_lookup = path(".")
depends: executable('tabix'), executable('bgzip')
input: genoFile, phenoFile
output: f'{cwd:a}/{path(_input[1]):bn}.sample_overlap.txt', f'{cwd:a}/{path(_input[1]):bn}.sample_genotypes.txt'
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[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    bash ${modular_script_dir}/data_preprocessing/genotype/GWAS_QC.sh genotype_phenotype_sample_overlap \
        --cwd "${cwd}" \
        --genoFile "${genoFile}" \
        --phenoFile "${phenoFile}" \
        --name "${name}" \
        --numThreads ${numThreads}