Genotype PLINK File Quality Control#
This workflow implements some preliminary data QC steps for PLINK input files. It supports both PLINK1 binary format (BED/BIM/FAM) and PLINK2 format (PGEN/PVAR/PSAM). VCF format of inputs will be converted to PLINK before performing QC.
Description#
This notebook includes workflow for
Compute kinship matrix in sample and estimate related individuals
Genotype and sample QC: by MAF, missing data and HWE
LD pruning for follow up PCA analysis on genotype, as needed
A potential limitation is that the workflow requires all samples and chromosomes to be merged as one single file, in order to perform both sample and variant level QC. However, in our experience using this pipeline with 200K exomes with 15 million variants, this pipeline works on the single merged PLINK file.
Methods#
Depending on the context of your problem, the workflow can be executed in two ways:
Run
qccommand to perform genotype data QC and LD pruning to generate a subset of variants in preparation for analysis such as PCA.Run
kingfirst on either the original or a subset of common variants to identify unrelated individuals. Thekingpipeline will split samples to related and unrelated individuals. Then you performqcon these individuals only and finally extract the same set of QC-ed variants for related individuals.
Default Parameters: QC#
Kinship coefficient for related individuals: 0.0625
MAF and MAC default: 0
Above default includes both common and are variant
Recommand MAF for PCA: 0.01, we should stick to common variants
Recommand MAC for single variant analysis: 5.
Variant level missingness threshold: 0.1
Sample level missingness threshold: 0.1
LD pruning via PLINK for PCA analysis:
window 50
shift 10
r2 0.1
HWE default: 1E-15 which is very lenient
Input#
The genotype files in either PLINK1 binary format (BED/BIM/FAM) or PLINK2 format (PGEN/PVAR/PSAM). For input in VCF format and/or per-chromosome VCF or PLINK format, please use vcf_to_plink and merge_plink in genotype formatting pipeline to convert them to PLINK file bundle.
Minimal Working Example#
Minimal working example data-set as well as the singularity container bioinfo.sif can be downloaded from Synapse.
The chr1_chr6 data-set was merged from chr1 and chr6 data, using merge_plink command from genotype formatting pipeline.
iii. Perform QC on both rare and common variants#
Timing: <1 min
sos run pipeline/GWAS_QC.ipynb qc_no_prune \
--cwd output/plink \
--genoFile output/plink/wgs.merged.bed \
--geno-filter 0.1 \
--mind-filter 0.1 \
--hwe-filter 1e-08 \
--mac-filter 0
INFO: Running qc_no_prune: Filter SNPs and select individuals
INFO: qc_no_prune is completed.
INFO: qc_no_prune output: /restricted/projectnb/xqtl/xqtl_protocol/toy_xqtl_protocol/output/plink/wgs.merged.plink_qc.bed
INFO: Workflow qc_no_prune (ID=w6697f77cea0f6dc2) is executed successfully with 1 completed step.
v. Sample match with genotype#
Timing: <1 min
sos run pipeline/GWAS_QC.ipynb genotype_phenotype_sample_overlap \
--cwd output/genotype/ \
--genoFile output/plink/wgs.merged.plink_qc.fam \
--phenoFile output/rnaseq/bulk_rnaseq_tmp_matrix.low_expression_filtered.outlier_removed.tmm.expression.bed.bed.gz
INFO: Running 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
INFO: genotype_phenotype_sample_overlap is completed.
INFO: genotype_phenotype_sample_overlap output: /restricted/projectnb/xqtl/xqtl_protocol/toy_xqtl_protocol/output/genotype/bulk_rnaseq_tmp_matrix.low_expression_filtered.outlier_removed.tmm.expression.bed.bed.sample_overlap.txt /restricted/projectnb/xqtl/xqtl_protocol/toy_xqtl_protocol/output/genotype/bulk_rnaseq_tmp_matrix.low_expression_filtered.outlier_removed.tmm.expression.bed.bed.sample_genotypes.txt
INFO: Workflow genotype_phenotype_sample_overlap (ID=w3667bce205c8538b) is executed successfully with 1 completed step.
vi. Kinship QC#
To accuratly estimate the PCs for the genotype. We split participants based on their kinship coefficients, estimated by KING.
Variant level and sample level QC on unrelated individuals using missingness > 10%, and LD-prunning in preparation for PCA analysis.
There is no related samples in these ROSMAP samples, so there is an additional step to only keep those samples in
rosmap_pheno.sample_genotypes.txtto do PCA.
Be aware:
If the message from king_2 shown as No related individuals detected from *.kin0, this means no related individuals detected for the samples in --keep_samples. In this case, there will be no output for related individuals from this step.
Timing: <2 min
sos run pipeline/GWAS_QC.ipynb king \
--cwd output/genotype/kinship \
--genoFile output/plink/wgs.merged.plink_qc.bed \
--name wgs.merged.king \
--keep-samples output/genotype/bulk_rnaseq_tmp_matrix.low_expression_filtered.outlier_removed.tmm.expression.bed.bed.sample_genotypes.txt
INFO: Running king_1: Inference of relationships in the sample to identify closely related individuals
INFO: king_1 is completed.
INFO: king_1 output: /restricted/projectnb/xqtl/xqtl_protocol/toy_xqtl_protocol/output/genotype/kinship/wgs.merged.plink_qc.wgs.merged.king.kin0
INFO: Running king_2: Select a list of unrelated individual with an attempt to maximize the unrelated individuals selected from the data
INFO: No related individuals detected from /restricted/projectnb/xqtl/xqtl_protocol/toy_xqtl_protocol/output/genotype/kinship/wgs.merged.plink_qc.wgs.merged.king.kin0.
INFO: king_2 output: /restricted/projectnb/xqtl/xqtl_protocol/toy_xqtl_protocol/output/genotype/kinship/wgs.merged.plink_qc.wgs.merged.king.related_id
INFO: Running king_3: Split genotype data into related and unrelated samples, if related individuals are detected
INFO: king_3 is completed.
INFO: king_3 output: /restricted/projectnb/xqtl/xqtl_protocol/toy_xqtl_protocol/output/genotype/kinship/wgs.merged.plink_qc.wgs.merged.king.unrelated.bed /restricted/projectnb/xqtl/xqtl_protocol/toy_xqtl_protocol/output/genotype/kinship/wgs.merged.plink_qc.wgs.merged.king.related.bed
INFO: Workflow king (ID=w414e098f8f3bd993) is executed successfully with 2 completed steps and 1 ignored step.
Command Interface#
sos run GWAS_QC.ipynb -h
usage: sos run 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:
--cwd output (as path)
the output directory for generated files
--name ''
A string to identify your analysis run
--genoFile paths
PLINK binary files
--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
--container ''
Software container option
--entrypoint ('micromamba run -a "" -n' + ' ' + re.sub(r'(_apptainer:latest|_docker:latest|\.sif)$', '', container.split('/')[-1])) if container else ""
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
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)
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
[global]
parameter: renovated_code_dir = path('renovated_code/script') # override with --renovated-code-dir
# the output directory for generated files
parameter: cwd = path("output")
# A string to identify your analysis run
parameter: name = ""
# 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
# Software container option
parameter: container = ""
parameter: entrypoint= ""
# 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)
Estimate kinship in the sample#
The output is a list of related individuals, as well as the kinship matrix
# 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', container = container, entrypoint = entrypoint
bash ${renovated_code_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', container = container, entrypoint = entrypoint
Rscript ${renovated_code_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
plink_command = get_plink_command_prefix(_input[1])
output_format = determine_plink_format(_input[1])
make_command = get_make_command(output_format)
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', container = container, entrypoint=entrypoint
plink2 \
${plink_command} ${_input[1]:n} \
--remove ${_input[0]} \
${('--keep %s' % keep_samples) if keep_samples.is_file() else ""} \
${make_command} \
--out ${_output[0]:n} \
--threads ${numThreads} \
--memory ${int(expand_size(mem) * 0.9)/1e6} --new-id-max-allele-len 1000 --set-all-var-ids chr@:#_\$r_\$a
if [ ${len(related_id)} -ne 0 ] ; then
plink2 \
${plink_command} ${_input[1]:n} \
--keep ${_input[0]} \
${make_command} \
--out ${_output[1]:n} \
--threads ${numThreads} \
--memory ${int(expand_size(mem) * 0.9)/1e6} --new-id-max-allele-len 1000 --set-all-var-ids chr@:#_\$r_\$a
else
touch ${_output[1]}
fi
Genotype and sample QC#
QC the genetic data based on MAF, sample and variant missigness and Hardy-Weinberg Equilibrium (HWE).
In this step you may also provide a list of samples to keep, for example in the case when you would like to subset a sample based on their ancestries to perform independent analyses on each of these groups.
The default parameters are set to reflect some suggestions in Table 1 of this paper.
# 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', container = container, entrypoint = entrypoint
bash ${renovated_code_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
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)
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', container = container, entrypoint = entrypoint
bash ${renovated_code_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} \
--numThreads ${numThreads}
Extract genotype based on overlap with phenotype#
This is an auxiliary step to match genotype and phenotype based on the data and look-up table. The look up table should contain two columns: sample_id, genotype_id. If the look up table is not provided or look-up table file not found, then we will assume the names have already been matched.
# 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', container = container, entrypoint = entrypoint
bash ${renovated_code_dir}/data_preprocessing/genotype/GWAS_QC.sh genotype_phenotype_sample_overlap \
--cwd "${cwd}" \
--genoFile "${genoFile}" \
--phenoFile "${phenoFile}" \
--name "${name}" \
--numThreads ${numThreads}