Genotype PLINK File Quality Control

Genotype PLINK File Quality Control#

This workflow implements some preliminary data QC steps for PLINK input files. 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 qc command to perform genotype data QC and LD pruning to generate a subset of variants in preparation for analysis such as PCA.
Run king first on either the original or a subset of common variants to identify unrelated individuals. The king pipeline will split samples to related and unrelated individuals. Then you perform qc on 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 whole genome PLINK bim/bed/fam bundle. 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#

sos run xqtl-protocol/pipeline/GWAS_QC.ipynb qc_no_prune \
   --cwd Genotype \
   --genoFile Genotype/ROSMAP_NIA_WGS.leftnorm.bcftools_qc.bed \
   --geno-filter 0.1 \
   --mind-filter 0.1 \
   --hwe-filter 1e-08   \
   --mac-filter 0 \
   --container /mnt/vast/hpc/csg/containers/bioinfo.sif \
   -J 1 -q csg -c csg.yml --mem 150G

v. Sample match with genotype#

Timing: <1 min

sos run pipeline/GWAS_QC.ipynb genotype_phenotype_sample_overlap \
        --cwd output/sample_meta \
        --genoFile input/protocol_example.genotype.chr21_22.fam  \
        --phenoFile input/protocol_example.protein.csv \
        --container containers/bioinfo.sif \
        --mem 5G

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:   /Users/alexmccreight/xqtl-protocol/output/sample_meta/protocol_example.protein.sample_overlap.txt /Users/alexmccreight/xqtl-protocol/output/sample_meta/protocol_example.protein.sample_genotypes.txt
INFO: Workflow genotype_phenotype_sample_overlap (ID=w71b4e35979654867) 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.txt to 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/kinship \
    --genoFile input/protocol_example.genotype.chr21_22.bed \
    --name pQTL \
    --keep-samples output/sample_meta/protocol_example.protein.sample_genotypes.txt \
    --container containers/bioinfo.sif \
    --no-maximize-unrelated \
    --mem 40G

INFO: Running king_1: Inference of relationships in the sample to identify closely related individuals
INFO: king_1 is completed.
INFO: king_1 output:   /Users/alexmccreight/xqtl-protocol/output/kinship/protocol_example.genotype.chr21_22.pQTL.kin0
INFO: Running king_2: Select a list of unrelated individual with an attempt to maximize the unrelated individuals selected from the data
INFO: king_2 is completed.
INFO: king_2 output:   /Users/alexmccreight/xqtl-protocol/output/kinship/protocol_example.genotype.chr21_22.pQTL.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:   /Users/alexmccreight/xqtl-protocol/output/kinship/protocol_example.genotype.chr21_22.pQTL.unrelated.bed /Users/alexmccreight/xqtl-protocol/output/kinship/protocol_example.genotype.chr21_22.pQTL.related.bed
INFO: Workflow king (ID=w7fad1d8b027ec781) is executed successfully with 3 completed steps.

vii. Prepare unrelated individuals data for PCA#

Here we write data to cache folder instead of output because this genotype data can be removed later after PCA. Also filter out minor allel accout < 5.

If your data has *.unrelated.bed generated, that means there are related individuals in your data. In cases, we will use output from the KING step for unrelated individuals.

Timing: <1 min

sos run pipeline/GWAS_QC.ipynb qc \
   --cwd output/cache \
   --genoFile output/kinship/protocol_example.genotype.chr21_22.pQTL.unrelated.bed \
   --mac-filter 5 \
   --container containers/bioinfo.sif \
   --mem 16G

INFO: Running basic QC filters: Filter SNPs and select individuals
INFO: basic QC filters is completed.
INFO: basic QC filters output:   /Users/alexmccreight/xqtl-protocol/output/cache/protocol_example.genotype.chr21_22.pQTL.unrelated.plink_qc.bed
INFO: Running LD pruning: LD prunning and remove related individuals (both ind of a pair) Plink2 has multi-threaded calculation for LD prunning
INFO: LD pruning is completed.
INFO: LD pruning output:   /Users/alexmccreight/xqtl-protocol/output/cache/protocol_example.genotype.chr21_22.pQTL.unrelated.plink_qc.prune.bed /Users/alexmccreight/xqtl-protocol/output/cache/protocol_example.genotype.chr21_22.pQTL.unrelated.plink_qc.prune.in
INFO: Workflow qc (ID=w3a34828bd2888342) is executed successfully with 2 completed steps.

In other cases eg ROSMAP proteomics data, message No related individuals detected from *.kin0 occured, there is no separate genotype data generated for unrelated individuals. In this case, we need to work from the original genotype data and must use --keep-samples to run qc to extract samples for PCA. For example:

Timing: <1 min

sos run pipeline/GWAS_QC.ipynb qc \
   --cwd output/cache \
   --genoFile input/protocol_example.genotype.chr21_22.bed \
   --keep-samples output/sample_meta/protocol_example.protein.sample_genotypes.txt \
   --name pQTL \
   --mac-filter 5 \
   --container containers/bioinfo.sif \
   --mem 40G

INFO: Running basic QC filters: Filter SNPs and select individuals
INFO: basic QC filters is completed.
INFO: basic QC filters output:   /Users/alexmccreight/xqtl-protocol/output/cache/protocol_example.genotype.chr21_22.pQTL.plink_qc.bed
INFO: Running LD pruning: LD prunning and remove related individuals (both ind of a pair) Plink2 has multi-threaded calculation for LD prunning
INFO: LD pruning is completed.
INFO: LD pruning output:   /Users/alexmccreight/xqtl-protocol/output/cache/protocol_example.genotype.chr21_22.pQTL.plink_qc.prune.bed /Users/alexmccreight/xqtl-protocol/output/cache/protocol_example.genotype.chr21_22.pQTL.plink_qc.prune.in
INFO: Workflow qc (ID=w3d8b01776519226f) is executed successfully with 2 completed steps.

[FIXME:] Extract previously selected variants from related individuals in preparation for PCA, only applying missingness filter at sample level,

sos run GWAS_QC.ipynb qc_no_prune \
    --cwd output/genotype \
    --genoFile output/genotype/chr1_chr6.20220110.related.bed \
    --keep-variants output/genotype/chr1_chr6.20220110.unrelated.for_pca.filtered.prune.in \
    --maf-filter 0 --geno-filter 0 --mind-filter 0.1 \
    --name for_pca \
    --container container/bioinfo.sif

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' + ' ' + container.split('/')[-1][:-4]) if container.endswith('.sif') 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
    Workflow Options:
      --[no-]maximize-unrelated (default to False)
                        If set to true, the unrelated individuals in a family
                        will be kept without being reported. Otherwise (use
                        `--no-maximize-unrelated`) the entire family will be
                        removed Note that attempting to maximize unrelated
                        individuals is computationally intensive on large 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]
# the output directory for generated files
parameter: cwd = path("output")
# A string to identify your analysis run
parameter: name = ""
# PLINK binary files
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 = ""
import re
parameter: entrypoint= ('micromamba run -a "" -n' + ' ' + re.sub(r'(_apptainer:latest|_docker:latest|\.sif)$', '', container.split('/')[-1])) if container else ""
# use this function to edit memory string for PLINK input
from sos.utils import expand_size
cwd = path(f"{cwd:a}")

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}'
bash: container=container, expand= "${ }", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', entrypoint=entrypoint
    plink2 \
      --bfile ${_input:n} \
      --make-king-table \
      --king-table-filter ${kinship} \
      ${('--keep %s' % keep_samples) if keep_samples.is_file() else ""} \
      ${('--remove %s' % remove_samples) if remove_samples.is_file() else ""} \
      --min-af ${kin_maf} \
      --max-af ${1-kin_maf} \
      --out ${_output:n} \
      --threads ${numThreads} \
      --memory ${int(expand_size(mem) * 0.9)/1e6} 
    
bash: expand= "${ }", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', container = container, entrypoint=entrypoint
    i="${_output}"
    output_size=$(ls -lh $i | cut -f 5 -d ' ')
    output_rows=$(cat $i | wc -l | cut -f 1 -d ' ')
    output_column=$(cat $i | head -1 | wc -w)
    output_preview=$(cat $i | grep -v "##" | head | cut -f 1,2,3,4,5,6)
    
    printf "output_info: %s\noutput_size: %s\noutput_rows: %s\noutput_column: %s\noutput_preview:\n%s\n" \
        "$i" "$output_size" "$output_rows" "$output_column" "$output_preview"

# Select a list of unrelated individual with an attempt to maximize the unrelated individuals selected from the data 
[king_2: shared = "related_id" ]
# If set to true, the unrelated individuals in a family will be kept without being reported. 
# Otherwise (use `--no-maximize-unrelated`) the entire family will be removed
# Note that attempting to maximize unrelated individuals is computationally intensive on large data.
parameter: maximize_unrelated = False
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   ## This should create an empty output such that both the empty relateness and unrelated samples bed will be created. So we dont need to subset the samples manually.
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}'
R:  container=container, expand= "${ }", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', entrypoint=entrypoint
    library(dplyr)
    library(igraph)
    # Remove related individuals while keeping maximum number of individuals
    # this function is simplified from: 
    # https://rdrr.io/cran/plinkQC/src/R/utils.R
    #' @param relatedness [data.frame] containing pair-wise relatedness estimates
    #' (in column [relatednessRelatedness]) for individual 1 (in column
    #' [relatednessIID1] and individual 2 (in column [relatednessIID1]). Columns
    #' relatednessIID1, relatednessIID2 and relatednessRelatedness have to present,
    #' while additional columns such as family IDs can be present. Default column
    #' names correspond to column names in output of plink --genome
    #' (\url{https://www.cog-genomics.org/plink/1.9/ibd}). All original
    #' columns for pair-wise highIBDTh fails will be returned in fail_IBD.
    #' @param relatednessTh [double] Threshold for filtering related individuals.
    #' Individuals, whose pair-wise relatedness estimates are greater than this
    #' threshold are considered related.
    relatednessFilter <- function(relatedness, 
                                  relatednessTh,
                                  relatednessIID1="IID1", 
                                  relatednessIID2="IID2",
                                  relatednessRelatedness="KINSHIP") {
        # format data
        if (!(relatednessIID1 %in% names(relatedness))) {
            stop(paste("Column", relatednessIID1, "for relatedness not found!"))
        }
        if (!(relatednessIID2 %in% names(relatedness))) {
            stop(paste("Column", relatednessIID1, "for relatedness not found!"))
        }
        if (!(relatednessRelatedness %in% names(relatedness))) {
            stop(paste("Column", relatednessRelatedness,
                       "for relatedness not found!"))
        }

        iid1_index <- which(colnames(relatedness) == relatednessIID1)
        iid2_index <- which(colnames(relatedness) == relatednessIID2)

        relatedness[,iid1_index] <- as.character(relatedness[,iid1_index])
        relatedness[,iid2_index] <- as.character(relatedness[,iid2_index])

        relatedness_names <- names(relatedness)
        names(relatedness)[iid1_index] <- "IID1"
        names(relatedness)[iid2_index] <- "IID2"
        names(relatedness)[names(relatedness) == relatednessRelatedness] <- "M"

        # Remove symmetric IID rows
        relatedness_original <- relatedness
        relatedness <- dplyr::select_(relatedness, ~IID1, ~IID2, ~M)

        sortedIDs <- data.frame(t(apply(relatedness, 1, function(pair) {
            c(sort(c(pair[1], pair[2])))
            })), stringsAsFactors=FALSE)
        keepIndex <- which(!duplicated(sortedIDs))

        relatedness_original <- relatedness_original[keepIndex,]
        relatedness <- relatedness[keepIndex,]

        # individuals with at least one pair-wise comparison > relatednessTh
        # return NULL to failIDs if no one fails the relatedness check
        highRelated <- dplyr::filter_(relatedness, ~M > relatednessTh)
        if (nrow(highRelated) == 0) {
            return(list(relatednessFails=NULL, failIDs=NULL))
        }

        # all samples with related individuals
        allRelated <- c(highRelated$IID1, highRelated$IID2)
        uniqueIIDs <- unique(allRelated)

        # Further selection of samples with relatives in cohort
        multipleRelative <- unique(allRelated[duplicated(allRelated)])
        singleRelative <- uniqueIIDs[!uniqueIIDs %in% multipleRelative]

        highRelatedMultiple <- highRelated[highRelated$IID1 %in% multipleRelative |
                                            highRelated$IID2 %in% multipleRelative,]
        highRelatedSingle <- highRelated[highRelated$IID1 %in% singleRelative &
                                           highRelated$IID2 %in% singleRelative,]

        # Only one related samples per individual
        if(length(singleRelative) != 0) {
          # randomly choose one to exclude
          failIDs_single <- highRelatedSingle[,1]
            
        } else {
          failIDs_single <- NULL
        }
  
        # An individual has multiple relatives
        if(length(multipleRelative) != 0) {
            relatedPerID <- lapply(multipleRelative, function(x) {
                tmp <- highRelatedMultiple[rowSums(
                    cbind(highRelatedMultiple$IID1 %in% x,
                          highRelatedMultiple$IID2 %in% x)) != 0,1:2]
                rel <- unique(unlist(tmp))
                return(rel)
            })
            names(relatedPerID) <- multipleRelative

            keepIDs_multiple <- lapply(relatedPerID, function(x) {
                pairwise <- t(combn(x, 2))
                index <- (highRelatedMultiple$IID1 %in% pairwise[,1] &
                              highRelatedMultiple$IID2 %in% pairwise[,2]) |
                    (highRelatedMultiple$IID1 %in% pairwise[,2] &
                         highRelatedMultiple$IID2 %in% pairwise[,1])
                combination <- highRelatedMultiple[index,]
                combination_graph <- igraph::graph_from_data_frame(combination,
                                                                   directed=FALSE)
                all_iv_set <- igraph::ivs(combination_graph)
                length_iv_set <- sapply(all_iv_set, function(x) length(x))

                if (all(length_iv_set == 1)) {
                    # check how often they occurr elsewhere
                    occurrence <- sapply(x, function(id) {
                        sum(sapply(relatedPerID, function(idlist) id %in% idlist))
                    })
                    # if occurrence the same everywhere, pick the first, else keep
                    # the one with minimum occurrence elsewhere
                    if (length(unique(occurrence)) == 1) {
                        nonRelated <- sort(x)[1]
                    } else {
                        nonRelated <- names(occurrence)[which.min(occurrence)]
                    }
                } else {
                    nonRelated <- all_iv_set[which.max(length_iv_set)]
                }
                return(nonRelated)
            })
            keepIDs_multiple <- unique(unlist(keepIDs_multiple))
            failIDs_multiple <- c(multipleRelative[!multipleRelative %in%
                                                       keepIDs_multiple])
        } else {
            failIDs_multiple <- NULL
        }
        allFailIIDs <- c(failIDs_single, failIDs_multiple)
        relatednessFails <- lapply(allFailIIDs, function(id) {
            fail_inorder <- relatedness_original$IID1 == id &
                relatedness_original$M > relatednessTh
            fail_inreverse <- relatedness_original$IID2 == id &
                relatedness_original$M > relatednessTh
            if (any(fail_inreverse)) {
                inreverse <- relatedness_original[fail_inreverse, ]
                id1 <- iid1_index
                id2 <- iid2_index
                inreverse[,c(id1, id2)] <- inreverse[,c(id2, id1)]
                names(inreverse) <- relatedness_names
            } else {
                inreverse <- NULL
            }
            inorder <- relatedness_original[fail_inorder, ]
            names(inorder) <- relatedness_names
            return(rbind(inorder, inreverse))
        })
        relatednessFails <- do.call(rbind, relatednessFails)
        if (nrow(relatednessFails) == 0) {
            relatednessFails <- NULL
            failIDs <- NULL
        } else {
            names(relatednessFails) <- relatedness_names
            rownames(relatednessFails) <- 1:nrow(relatednessFails)
            uniqueFails <- relatednessFails[!duplicated(relatednessFails[,iid1_index]),]
            failIDs <- uniqueFails[,iid1_index]
        }
        return(list(relatednessFails=relatednessFails, failIDs=failIDs))
    }
    
  
    # main code
    kin0 <- read.table(${_input:r}, header=F, stringsAsFactor=F)
    colnames(kin0) <- c("FID1","ID1","FID2","ID2","NSNP","HETHET","IBS0","KINSHIP")
    if (${"TRUE" if maximize_unrelated else "FALSE"}) {
        rel <- relatednessFilter(kin0, ${kinship}, "ID1", "ID2", "KINSHIP")$failIDs
        tmp1 <- kin0[,1:2]
        tmp2 <- kin0[,3:4]
        colnames(tmp1) = colnames(tmp2) = c("FID", "ID")
        # Get the family ID for these rels so there are two columns FID and IID in the output
        lookup <- dplyr::distinct(rbind(tmp1,tmp2))
        dat <- lookup[which(lookup[,2] %in% rel),]
    } else {
        rel <- kin0 %>% filter(KINSHIP >= ${kinship})
        dat = rbind(rel[,c("FID1","ID1")],setNames(rel[,c("FID2","ID2")],c("FID1","ID1")))
        dat = dat[!duplicated(dat),] ## This is to remove duplicated FID and IID caused by one sample being related to multiple samples
       }    

    cat("There are", nrow(dat),"related individuals using a kinship threshold of ${kinship}\n")
    write.table(dat,${_output:r}, quote=FALSE, row.names=FALSE, col.names=FALSE)
    
bash: expand= "${ }", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', container = container, entrypoint=entrypoint
    i="${_output}"
    output_size=$(ls -lh $i | cut -f 5 -d ' ')
    output_rows=$(cat $i | wc -l | cut -f 1 -d ' ')
    output_column=$(cat $i | head -1 | wc -w)
    output_preview=$(cat $i | grep -v "##" | head | cut -f 1,2,3,4,5,6)
    
    printf "output_info: %s\noutput_size: %s\noutput_rows: %s\noutput_column: %s\noutput_preview:\n%s\n" \
        "$i" "$output_size" "$output_rows" "$output_column" "$output_preview"

# 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: unrelated_bed = f'{cwd}/{_input[0]:bn}.unrelated.bed',
        related_bed = f'{cwd}/{_input[0]:bn}.related.bed'
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 \
      --bfile ${_input[1]:n} \
      --remove ${_input[0]} \
      ${('--keep %s' % keep_samples) if keep_samples.is_file() else ""} \
      --make-bed \
      --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 \
      --bfile ${_input[1]:n} \
      --keep ${_input[0]} \
      --make-bed \
      --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

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
output: f'{cwd}/{_input:bn}{("."+name) if name else ""}.plink_qc{".extracted" if keep_variants.is_file() else ""}{".bed" if not meta_only else ".snplist"}'
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
    plink2 \
      --bfile ${_input:n} \
      ${('--maf %s' % maf_filter) if maf_filter > 0 else ''} \
      ${('--max-maf %s' % maf_max_filter) if maf_max_filter > 0 else ''} \
      ${('--mac %s' % mac_filter) if mac_filter > 0 else ''} \
      ${('--max-mac %s' % mac_max_filter) if mac_max_filter > 0 else ''} \
      ${('--geno %s' % geno_filter) if geno_filter > 0 else ''} \
      ${('--hwe %s' % hwe_filter) if hwe_filter > 0 else ''} \
      ${('--mind %s' % mind_filter) if mind_filter > 0 else ''} \
      ${('--keep %s' % keep_samples) if keep_samples.is_file() else ""} \
      ${('--remove %s' % remove_samples) if remove_samples.is_file() else ""} \
      ${('--exclude %s' % exclude_variants) if exclude_variants.is_file() else ""} \
      ${('--extract %s' % keep_variants) if keep_variants.is_file() else ""} \
      ${('--make-bed') if not meta_only else "--write-snplist --write-samples"} \
      ${("") if not rm_dups else "--rm-dup force-first 'list'"} \
      ${paths(["--%s" % x for x in other_args]) if other_args else ""} \
      --out ${_output:n} \
      --threads ${numThreads} \
      --memory ${int(expand_size(mem) * 0.9)/1e6} --new-id-max-allele-len 1000 --set-all-var-ids chr@:#_\$r_\$a 
        
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint=entrypoint
    i="${_output}"
    output_size=$(ls -lh $i | cut -f 5 -d ' ')
    printf "output_info: %s\noutput_size: %s\n" "$i" "$output_size" >> ${_output:n}.stdout

# 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)
output: bed=f'{cwd}/{_input:bn}.prune.bed', 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: container=container, expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', entrypoint=entrypoint
    plink2 \
    --bfile ${_input:n} \
    --indep-pairwise ${window} ${shift} ${r2}  \
    --out ${_output["prune"]:nn} \
    --threads ${numThreads} \
    --memory ${int(expand_size(mem) * 0.9)/1e6}
   
    plink2 \
    --bfile ${_input:n} \
    --extract ${_output['prune']} \
    --make-bed \
    --out ${_output['bed']:n} \
    --threads ${numThreads} \
    --memory ${int(expand_size(mem) * 0.9)/1e6}
    
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint=entrypoint
    i="${_output[0]}"
    output_size=$(ls -lh $i | cut -f 5 -d ' ')
    printf "output_info: %s\noutput_size: %s\n" "$i" "$output_size" >> ${_output[0]:n}.stdout
    i="${_output[1]}"
    output_size=$(ls -lh $i | cut -f 5 -d ' ')
    output_rows=$(zcat $i | wc -l | cut -f 1 -d ' ')
    output_column=$(zcat $i | head -1 | wc -w)
    output_preview=$(cat $i | grep -v "##" | head | cut -f 1,2,3,4,5,6) 
    printf "output_info: %s\noutput_size: %s\noutput_rows: %s\noutput_column: %s\noutput_preview:\n%s\n" \
        "$i" "$output_size" "$output_rows" "$output_column" "$output_preview" >> ${_output[1]}.stdout

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}'
R: expand = "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint=entrypoint
    # Load required libraries
    library(dplyr)
    library(readr)
    library(data.table)

    # Read data files; let read_delim auto-determine the delimiter
    genoFam <- fread(${_input[0]:ar}, header=FALSE)
    phenoFile <- read_delim(${_input[1]:ar}, col_names=TRUE)
    if (${"TRUE" if sample_participant_lookup.is_file() else "FALSE"}) {
        sample_lookup <- fread(${sample_participant_lookup:ar}, header=TRUE)
        colnames(sample_lookup) <- c("sample_id", "genotype_id") # FIXME: This is for the old lookup table with columns c("sample_id", "participant_id")
        # rename phenotype file according to lookup file
        colnames(phenoFile)[-c(1:4)] <- phenoFile %>%
          colnames() %>%
          .[-c(1:4)] %>%
          match(sample_lookup$sample_id) %>%
          sample_lookup$genotype_id[.]
      output_file <- paste0(${phenoFile:nnr}, '.rename_sample.bed')
      # rename phenotype file
      phenoFile$start <- as.integer(phenoFile$start) 
      phenoFile$end <- as.integer(phenoFile$end) 
      fwrite(phenoFile, output_file, sep = '\t')
      # Compress the file using bgzip
      system(paste("bgzip -c", output_file, ">", paste0(output_file, ".gz")))
      # Create the index file using tabix
      system(paste("tabix -p bed", paste0(output_file, ".gz")))
      system(paste("rm", output_file))
    } else {
        sample_lookup <- cbind(genoFam[,2], genoFam[,2])
        colnames(sample_lookup) <- c("genotype_id", "sample_id")
    }
    sample_lookup <- sample_lookup %>%
    filter(
        genotype_id %in% genoFam$V2,
        sample_id %in% colnames(phenoFile)
    )
    
    genoFam %>%
    filter(
        V2 %in% sample_lookup$genotype_id,
    ) %>%
    select(V1, V2) %>%
    fwrite(${_output[1]:r}, col.names=FALSE, sep="\t")

    sample_lookup %>%
    fwrite(${_output[0]:r}, sep="\t")

Genotype PLINK File Quality Control

Contents

Genotype PLINK File Quality Control#

Description#

Methods#

Default Parameters: QC#

Input#

Minimal Working Example#

iii. Perform QC on both rare and common variants#

v. Sample match with genotype#

vi. Kinship QC#

vii. Prepare unrelated individuals data for PCA#

Command Interface#

Estimate kinship in the sample#

Genotype and sample QC#

Extract genotype based on overlap with phenotype#