Principal Component Analysis

Principal Component Analysis#

The intention of this notebook is to perform the PCA analysis on genotype data and generate plots.

Overview#

Population structure is the classic confounder in genetic association: if ancestry correlates with both genotype and phenotype, unadjusted tests return associations that are real but not causal. The remedy is to compute principal components of the genotype matrix and carry the leading ones as covariates.

Components are computed on unrelated individuals and the remaining related samples are projected back into that space, so relatives cannot distort the axes but every sample still gets coordinates. The sequence is: remove related individuals, LD-prune the variants, run PCA on the unrelated set, then exclude PCA-space outliers. Relatedness estimation and sample QC happen upstream in GWAS_QC.ipynb.

The workflows cover the whole job rather than just the decomposition: pca_plink and flashpca compute components, project_samples places new samples onto axes computed elsewhere, plot_pca produces the scree and scatter plots you inspect, and detect_outliers flags samples that sit far from any cluster.

Steps to generate a PCA include:

  • removing related individuals

  • pruning variants in linkage disequilibrium (LD)

  • performing PCA analysis on the genotypes of the unrelated individuals

  • excluding outlier samples in the PCA space for individuals of homogeneous self-reported ancestry. These outliers may suggest poor genotyping quality or distant relatedness.

Limitations:

  1. Some of the PCs may capture LD structure rather than population structure, which decreases power to detect associations in these regions of high LD.

  2. When projecting a new study dataset onto the PCA space computed from a reference dataset, the projected PCs are shrunk toward 0 in the new dataset.

  3. PC scores may capture outliers arising from family structure, population structure or other reasons. Detecting and removing these individuals can be worthwhile, either to maximise the population structure captured by the PCA when only a few outliers are involved, or to restrict analyses to genetically homogeneous samples.

When to run it. After genotype QC and LD pruning, before covariate preprocessing; the components become part of the covariate file used by the association scan.

Input#

  • --genoFile output/gwas_qc/genotype/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.plink_qc.prune.bed (the LD-pruned PLINK bundle for the unrelated individuals, as produced by GWAS QC. PCA is computed on these samples)

  • --phenoFile: a phenotype or label table carrying an IID column, optionally FID, used to label and colour the projected samples. A PLINK .fam may be used instead. Its extra columns are what --pop-col and --label-col select from. Example input/covariate/protocol_example.pca_pheno.txt, 60 samples and 7 columns:

    FID         IID         msex  age_death  pmi    study  race
SAMPLE_001  SAMPLE_001  1     90.97      10.57  1      1
SAMPLE_002  SAMPLE_002  1     80.24      7.87   0      1
SAMPLE_003  SAMPLE_003  1     83.9       9.93   0      1

  • --pca-model: the fitted PCA model to project onto, used by the projection step rather than the fitting step.

  • --name: the stem of the output files. Required.

  • --cwd: the directory outputs are written to.

  • --k: how many principal components to compute, 20 by default.

  • --maha-k: how many components the Mahalanobis outlier distance is computed over, 5 by default.

  • --pop-col and --pops: the phenotype column naming each sample’s population and the subset of populations to analyse, for admixed or multi-ancestry data.

  • --label-col: the phenotype column used to colour points in the PCA plots.

  • --keep-samples and --remove-samples: sample lists to restrict or exclude before fitting.

  • --keep-variants: a variant list to restrict to.

  • --maf-filter, --mac-filter, --geno-filter and --mind-filter: the PLINK frequency, count, variant-missingness and sample-missingness thresholds applied before PCA.

Output#

  • {cwd}/{name}.pca.rds and {cwd}/{name}.pca.txt - the fitted PCA model and the PC scores for the unrelated samples. Every path below sits under --cwd and shares the stem set by --name.

  • {name}.pca.scree.txt and {name}.pca.scree.png - variance explained per component, as a table and a plot. The table is what covariate_formatting reads to choose how many PCs to keep:

    PCs  PVE   PVE_cum
1    0.05  0.05
2    0.05  0.1

  • {name}.pca.mahalanobis.rds, .mahalanobis_hist.png and .mahalanobis_qq.png - the per-population Mahalanobis distances and their diagnostic plots.

  • {name}.pca.outliers - the samples flagged as PC outliers, an index and a sample ID per line:

    0  SAMPLE_001
0  SAMPLE_026
0  SAMPLE_050

  • {name}.pca.pc.png - the PCA plot.

  • {name}.projected.rds, .projected.txt and .projected.pc.png - the related individuals projected back onto the model, and the plot including them.

  • {name}.pca.analysis_summary.md - a written summary of the run.

The pca_plink sanity-check workflow declares {genoFile}.pca.eigenvec, but no .eigenvec file ships in output/pca_uf; only the flashpca products are present.

Minimal Working Example#

Note: parameters set for the MWE are meant to let the MWE work to show the workflow procedures. They may be unrealistic and should not be used in practice. The pipeline has reasonable default values for what is suggested in practice for most of the parameters.

Step 1. Estimate kinship in the sample (prerequisite)#

Before PCA it is necessary to know which individuals are related, because PCA must be computed on an unrelated subset and the related individuals are projected back afterwards. This step runs the king workflow from GWAS_QC to estimate pairwise kinship and split the samples into unrelated and related sets. It is an upstream prerequisite — its outputs (the king.unrelated and king.related PLINK bundles) feed the PCA steps below.

Timing: TBD (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 2. Sample selection and QC of the genotype data (prerequisite)#

QC the genotypes that will go into PCA: filter on MAF, sample/variant missingness and Hardy-Weinberg equilibrium, and LD-prune the variants. You can restrict to one population or to the related / unrelated split. Here the QC and LD-pruning are applied to the unrelated individuals to build the PCA reference, and separately extract the same variants for the related individuals so they can be projected later. These qc / qc_no_prune calls come from GWAS_QC and are upstream prerequisites for the PCA steps.

Variant level and sample level QC on unrelated individuals, in preparation for PCA analysis#

Timing: TBD (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 \
    --mac-filter 5

Extract the previously selected (LD-pruned) variants from the related individuals, applying only a sample-level missingness filter so the related samples carry exactly the same variant set as the unrelated reference.

Timing: TBD (on toy dataset)

sos run pipeline/GWAS_QC.ipynb qc_no_prune \
    --cwd output/pca_related \
    --genoFile output/gwas_qc/kinship/protocol_example.genotype.merged.plink_qc.protocol_example.king.related.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

Step 3. Run PCA on the unrelated samples (flashpca)#

This is the core PCA step. It runs flashpca on the LD-pruned unrelated PLINK bundle to compute the principal components, then automatically computes Mahalanobis distances and flags PC outliers. Outputs: the fitted PCA model (.pca.rds), the PC score table (.pca.txt), a scree plot of variance explained, and the outlier list. Timing: <2 min.

The --name argument tags all output files; point --genoFile at the LD-pruned unrelated bed produced upstream.

Timing: TBD (on toy dataset)

sos run pipeline/PCA.ipynb flashpca \
    --cwd output/pca_uf \
    --genoFile output/gwas_qc/genotype/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.plink_qc.prune.bed \
    --name protocol_example

What you should see: the workflow runs flashpca_1 (PCA), then detect_outliers (Mahalanobis distance per population) and plot_pca, finishing with “executed successfully with 5 completed steps”. The PCA model and PC scores are written to the --cwd directory. The scree plot below shows the proportion of variance explained by each PC — use it to decide how many PCs to keep.

%preview output/pca_uf/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.plink_qc.prune.protocol_example.pca.scree.png

Analysis by population (admixed / multi-ancestry data)#

The Steps 1–4 above produce a single PCA for the whole cohort, which is appropriate for a homogeneous population. If your cohort contains multiple ancestry groups, run the per-population workflow below instead: split the samples by population, then repeat the same QC → PCA → projection steps separately within each population.

pheno <- read.table("input/covariate/protocol_example.pca_pheno.txt", header = TRUE, stringsAsFactors = FALSE)

if (all(pheno$FID == pheno$IID)) pheno$FID <- "0"

for (i in 1:3) {
    race    <- subset(pheno, race == i)
    race_id <- cbind(race$FID, race$IID)
    write.table(race_id,
                paste0("output/pca_uf/protocol_example.ID.", i, ".txt"),
                quote = FALSE, sep = "\t", col.names = FALSE, row.names = FALSE)
}

Step 1. QC on unrelated samples, per population#

Same as Step 2 previously (QC), but looped over each population. Below only Populations 1 and 3 are shown as an example.

for i in 3; do
    sos run pipeline/GWAS_QC.ipynb qc \
        --cwd output/pca_uf \
        --genoFile output/gwas_qc/genotype/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.plink_qc.bed \
        --keep-samples output/pca_uf/protocol_example.ID.$i.txt \
        --mac-filter 5 \
        --bad-ld True \
        --name pop_$i
done

Step 3. Run PCA (flashpca) on unrelated samples, per population#

Same as Step 3 above, looped per population. Here the number of PCs is set to 5 (--k 5) because each population’s sample size is small.

for i in 3; do
    sos run pipeline/PCA.ipynb flashpca \
        --name pop_$i \
        --cwd output/pca_uf \
        --genoFile output/pca_uf/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.plink_qc.pop_$i.plink_qc.prune.bed \
        --phenoFile input/covariate/protocol_example.pca_pheno.txt \
        --label-col race \
        --pop-col race \
        --maha-k 2 \
        --k 5
done

Command Interface#

sos run pipeline/PCA.ipynb -h

usage: sos run pipeline/PCA.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:
  pca_plink
  flashpca
  project_samples
  plot_pca
  detect_outliers

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
  --pop-col ''
                        Name of the population column in the phenoFile
  --pops  (as list)
                        Name of the populations (from the population column) you
                        would like to plot and show on the PCA plot
  --label-col ''
                        Name of the color label column in the phenoFile; can be
                        the same as population column. Can also be a separate
                        column eg a "super population" column as a way to enable
                        you to combine selected populations based on another
                        column.
  --k 20 (as int)
                        Number of Principal Components to output,must be
                        consistant between flashpca run and project samples run
                        (flashpca partial PCA method).
  --maha-k 5 (as int)
                        Number of Principal Components based on which outliers
                        should be evaluated. Default is 5 but this should be
                        based on examine the scree plot
  --[no-]homogeneous (default to False)
                        Homogeneity of populations. Set to --homogeneous when
                        true and --no-homogeneous when false
  --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 10 (as int)
                        Number of threads

Sections
  pca_plink:            PCA command with PLINK, as a sanity check
    Workflow Options:
      --genoFile VAL (as path, required)
                        PLINK binary file
  flashpca_1:           Run PCA analysis using flashpca
    Workflow Options:
      --genoFile VAL (as path, required)
                        Plink binary file
      --phenoFile  path(f'{genoFile}'.replace(".bed", ".fam").replace(".pgen", ".psam"))

                        The phenotypic file
      --min-pop-size 2 (as int)
                        minimum population size to consider in the analysis
      --stand binom2
                        How to standardize X before PCA
  project_samples_1:    Project back to PCA model additional samples
    Workflow Options:
      --genoFile VAL (as path, required)
                        Plink binary file
      --phenoFile  path(f'{genoFile}'.replace(".bed", ".fam").replace(".pgen", ".psam"))

                        The phenotypic file
      --stand binom2
                        How to standardize X before PCA
      --pca-model  f'{cwd}/{phenoFile:bn}{("."+name) if name else ""}.{(suffix+".") if suffix != "" else ""}pca.rds'

  plot_pca:             Plot PCA results. Can be used independently as
                        "plot_pca" or combined with other workflow as eg
                        "flashpca+plot_pca"
    Workflow Options:
      --outlier-file . (as path)
      --plot-data VAL (as path, required)
      --min-axis ''
      --max-axis ''
  detect_outliers:      Calculate Mahalanobis distance per population and report
                        outliers
    Workflow Options:
      --prob 0.997 (as float)
                        Set the probability to remove outliers eg 0.95 or 0.997
      --pval 0.05 (as float)
                        Mahalanobis distance p-value cutoff
      --[no-]robust (default to True)
                        Robust Mahalanobis to outliers
      --pca-result VAL (as path, required)
  flashpca_2, project_samples_2:
    Workflow Options:
      --prob 0.997 (as float)
                        Set the probability to remove outliers eg 0.95 or 0.997
      --[no-]robust (default to True)
                        Robust Mahalanobis to outliers
  flashpca_3, project_samples_3:

Workflow implementation#

[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
# Name of the population column in the phenoFile
parameter: pop_col = ""
# Name of the populations (from the population column) you would like to plot and show on the PCA plot
parameter: pops = []
# Name of the color label column in the phenoFile; can be the same as population column. Can also be a separate column eg a "super population" column as a way to enable you to combine selected populations based on another column.
parameter: label_col = ""
# Number of Principal Components to output,must be consistant between flashpca run and project samples run (flashpca partial PCA method).
parameter: k = 20
# Number of Principal Components based on which outliers should be evaluated. Default is 5 but this should be based on examine the scree plot
parameter: maha_k = 5
# Homogeneity of populations. Set to --homogeneous when true and --no-homogeneous when false
parameter: homogeneous = False
import re
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "5h"
# Memory expected
parameter: mem = "16G"
# Number of threads
parameter: numThreads = 10
suffix = '_'.join(pops)
cwd = path(f"{cwd:a}")
if not pop_col:
    homogeneous = True

# 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"

def normalize_phenoFile(pheno_path, out_dir):
    """Return a phenoFile path where FID=0 if the file has FID==IID in every data row.
    Writes a normalized copy to out_dir when needed; otherwise returns the original path."""
    import os, csv
    pheno_path = str(pheno_path)
    with open(pheno_path) as fh:
        rows = [r for r in csv.reader(fh, delimiter='\t')]
    if not rows:
        return path(pheno_path)
    # Detect header: first row is a header when its first cell is not a sample value
    has_header = rows[0][0].upper() in ('FID', '#FID', 'IID', 'SAMPLE', 'SAMPLE_ID')
    data_rows = rows[1:] if has_header else rows
    if not data_rows or not all(len(r) >= 2 and r[0] == r[1] for r in data_rows):
        return path(pheno_path)
    # FID==IID everywhere — write a normalized copy to out_dir with FID=0
    out_path = os.path.join(str(out_dir), os.path.basename(pheno_path))
    with open(out_path, 'w', newline='') as fh:
        w = csv.writer(fh, delimiter='\t')
        if has_header:
            w.writerow(rows[0])
        for r in data_rows:
            r[0] = '0'
            w.writerow(r)
    return path(out_path)

# PCA command with PLINK, as a sanity check
[pca_plink]
# PLINK binary file
parameter: genoFile = path
input: genoFile
plink_command = get_plink_command_prefix(_input)
output: f'{cwd}/{genoFile:bn}.pca.eigenvec'
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'
    plink2 ${plink_command} ${_input:n} --out ${_output:n} --pca ${k}

PCA analysis#

# Run PCA analysis using flashpca 
[flashpca_1]
# Plink binary file
parameter: genoFile = path
# The phenotypic file
parameter: phenoFile = path(f'{genoFile}'.replace(".bed", ".fam").replace(".pgen", ".psam"))
# minimum population size to consider in the analysis
parameter: min_pop_size = 2
# How to standardize X before PCA
parameter: stand = "binom2"
## Input genoFile here is for unrelated samples
phenoFile = normalize_phenoFile(phenoFile, cwd)
input: genoFile, phenoFile
output: f'{cwd}/{phenoFile:bn}{("."+name) if name else ""}.{(suffix+".") if suffix != "" else ""}pca.rds'
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:n}.stderr', stdout = f'{_output:n}.stdout'
    Rscript ${modular_script_dir}/data_preprocessing/genotype/PCA.R \
        --step flashpca \
        --cwd "${cwd}" \
        --genoFile "${_input[0]}" \
        --phenoFile "${_input[1]}" \
        --output "${_output}" \
        --stand "${stand}" \
        --min-pop-size ${min_pop_size} \
        --homogeneous "${homogeneous}" \
        --pop-col "${pop_col}" \
        --label-col "${label_col}" \
        --pops "${",".join([str(x) for x in pops])}" \
        --k ${k} \
        --maha-k ${maha_k} \
        --numThreads ${numThreads}

Plot PCA results#

# Plot PCA results. 
# Can be used independently as "plot_pca" or combined with other workflow as eg "flashpca+plot_pca"
[plot_pca]
parameter: outlier_file = path()
parameter: plot_data = path
parameter: min_axis = ""
parameter: max_axis = ""
input: plot_data
output: f'{cwd}/{_input:bn}.pc.png',
        f'{cwd}/{_input:bn}.scree.png',
        f'{cwd}/{_input:bn}.scree.txt'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = 1, tags = f'{step_name}_{_output[0]:bn}'
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    Rscript ${modular_script_dir}/data_preprocessing/genotype/PCA.R \
        --step plot_pca \
        --cwd "${cwd}" \
        --plot-data "${_input}" \
        --outlier-file "${outlier_file}" \
        --min-axis "${min_axis}" \
        --max-axis "${max_axis}" \
        --pop-col "${pop_col}" \
        --label-col "${label_col}" \
        --pops "${",".join([str(x) for x in pops])}" \
        --k ${k} \
        --output-pc-plot "${_output[0]}" \
        --output-scree-plot "${_output[1]}" \
        --output-scree-text "${_output[2]}" \
        --numThreads ${numThreads}

Detect outliers#

# Calculate Mahalanobis distance per population and report outliers
[detect_outliers]
# Set the probability to remove outliers eg 0.95 or 0.997
parameter: prob = 0.997
# Mahalanobis distance p-value cutoff
parameter: pval = 0.05
# Robust Mahalanobis to outliers
parameter: robust = True
parameter: pca_result = path
input: pca_result
output: distance=f'{_input:n}.mahalanobis.rds',
        identified_outliers=f'{_input:n}.outliers',
        analysis_summary=f'{_input:n}.analysis_summary.md',
        qqplot_mahalanobis=f'{_input:n}.mahalanobis_qq.png',
        hist_mahalanobis=f'{_input:n}.mahalanobis_hist.png'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = 1, tags = f'{step_name}_{_output[0]:bn}'
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    Rscript ${modular_script_dir}/data_preprocessing/genotype/PCA.R \
        --step detect_outliers \
        --cwd "${cwd}" \
        --pca-result "${_input}" \
        --prob ${prob} \
        --pval ${pval} \
        --robust ${"TRUE" if robust else "FALSE"} \
        --pop-col "${pop_col}" \
        --k ${k} \
        --distance-output "${_output['distance']}" \
        --identified-outliers-output "${_output['identified_outliers']}" \
        --analysis-summary-output "${_output['analysis_summary']}" \
        --qqplot-output "${_output['qqplot_mahalanobis']}" \
        --hist-output "${_output['hist_mahalanobis']}" \
        --numThreads ${numThreads}

Add plot and outlier detection to PCA steps#

[flashpca_2, project_samples_2]
# Set the probability to remove outliers eg 0.95 or 0.997
parameter: prob = 0.997
# Robust Mahalanobis to outliers
parameter: robust = True
output: distance=f'{_input:n}.mahalanobis.rds',
        identified_outliers=f'{_input:n}.outliers',
        analysis_summary=f'{_input:n}.analysis_summary.md',
        qqplot_mahalanobis=f'{_input:n}.mahalanobis_qq.png',
        hist_mahalanobis=f'{_input:n}.mahalanobis_hist.png'
sos_run("detect_outliers", pca_result=_input, prob=prob, robust=robust)

[flashpca_3, project_samples_3]
input: output_from(1), output_from(2)['identified_outliers']
outliers = [x.strip() for x in open(_input[1]).readlines() if x.strip()]
output: f"{cwd}/{_input[0]:bn}.pc.png",
        f"{cwd}/{_input[0]:bn}.scree.png"
sos_run("plot_pca", plot_data = _input[0], outlier_file = _input[1] if len(outliers) else path())