GWAS MLMA Batch QC

This repository provides a server-ready workflow for batch quality control of GCTA/GEMMA-style GWAS .mlma summary result files across multiple populations and traits.

It is designed for use cases like:

• thousands of .mlma result files
• automated Manhattan plot, QQ plot, lambda GC, P-value sanity checks, and QC flagging

What The Pipeline Produces

For every GWAS result file, the pipeline writes:

• combined Manhattan + QQ plot
• optional compressed cleaned summary table
• per-trait QC metrics

For the whole batch, it writes:

• qc_summary.tsv: one row per GWAS result
• qc_fail.tsv: results requiring review
• qc_warn.tsv: borderline results
• qc_pass.tsv: results passing automated checks
• qc_recommendations.tsv: likely problem pattern and suggested adjustment

Recommended Directory Layout On Server

project/
  gwas-results/
    pop1/
    pop2/
    pop3/
  manifest.tsv
  gwas-mlma-batch-qc/

The manifest controls which files are analyzed.

Manifest Format

Create a tab-separated manifest.tsv with these columns:

population	trait	file
pop1	height	/path/to/gwas-results/pop1/height.mlma
pop2	height	/path/to/gwas-results/pop2/height.mlma
pop3	height	/path/to/gwas-results/pop3/height.mlma

Optional columns are allowed and preserved in the output, for example sample_size, batch, or note.

An example file is provided at examples/manifest.example.tsv.

Install Fast Python Dependencies

The recommended server workflow uses the Python fast script.

conda install -c conda-forge numpy scipy matplotlib polars

Optional GPU acceleration for QQ sorting requires a CUDA-matched CuPy build. Example for CUDA 12:

conda install -c conda-forge cupy cuda-version=12.0

If polars is unavailable, the script can still run with a slower Python csv fallback, but installing Polars is strongly recommended. For whitespace-delimited GCTA .mlma files, the fast script uses numpy.loadtxt(usecols=...) to read only the required columns. Use Python 3.8 or newer.

GPU note: .mlma text parsing is usually the main bottleneck, so GPU mainly helps the QQ sorting step for very large files. It does not remove the need to read the text files from disk.

Install R Dependencies

On most Linux servers:

Rscript -e 'install.packages(c("data.table", "ggplot2"), repos="https://cloud.r-project.org")'

No other R packages are required.

The R script is kept as a stable fallback. For large batches, prefer scripts/gwas_mlma_qc_fast.py.

Three-Step QC Workflow

For large batches, use the three-step workflow below. It is much faster than drawing Manhattan plots for every .mlma file.

Step 1: Metrics

This reads each .mlma file, calculates QC metrics, and does not draw plots.

Full run:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_step1_metrics \
  --threads 8 \
  --mode metrics

Partial run by population:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_step1_metrics_C4 \
  --threads 4 \
  --mode metrics \
  --only-pop C4

Partial run for a quick first-N test:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_step1_metrics_first4 \
  --threads 4 \
  --mode metrics \
  --limit 4

For pilot tests on sampled .mlma files, use --assessment-mode sample. This still calculates lambda GC, minimum P value, valid P count, invalid P count, and significant hit count, but it skips final thresholds that require the complete file, such as total SNP count and excess significant hits.

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest_sample100k.tsv \
  --outdir qc_sample100k_metrics \
  --threads 4 \
  --mode metrics \
  --assessment-mode sample

To calculate metrics only and avoid all PASS/WARN/FAIL classification:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest_sample100k.tsv \
  --outdir qc_sample100k_calculate \
  --threads 4 \
  --mode metrics \
  --assessment-mode calculate

Main outputs:

• qc_summary.tsv
• qc_pass.tsv
• qc_warn.tsv
• qc_fail.tsv
• qc_calculated.tsv
• qc_quality_solutions.tsv

Step 2: QQ Plots

This draws QQ plots for all manifest rows and classifies the QQ shape.

Full run:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_step2_qq \
  --threads 4 \
  --mode qq \
  --qq-max-points 100000 \
  --qq-confidence 0.95

GPU QQ sorting, if CuPy is installed:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_step2_qq_gpu \
  --threads 2 \
  --mode qq \
  --qq-engine gpu \
  --qq-max-points 100000

Use fewer threads with GPU mode so multiple processes do not fight for the same GPU memory.

Two-GPU QQ Workflow

When a server has 2 GPUs, split the manifest into two files and bind one Python process to each GPU. This is usually more stable than using many threads with one process.

Create split manifests:

mkdir -p QC_02_2026_results/gpu_split

head -n 1 QC_02_2026_results/manifest.tsv > QC_02_2026_results/gpu_split/manifest_gpu0.tsv
head -n 1 QC_02_2026_results/manifest.tsv > QC_02_2026_results/gpu_split/manifest_gpu1.tsv

awk 'NR==1{next} NR%2==0{print >> "QC_02_2026_results/gpu_split/manifest_gpu0.tsv"} NR%2==1{print >> "QC_02_2026_results/gpu_split/manifest_gpu1.tsv"}' \
  QC_02_2026_results/manifest.tsv

Test each GPU with a small subset first:

CUDA_VISIBLE_DEVICES=0 python GWAS_QC/scripts/gwas_mlma_qc_fast.py \
  --manifest QC_02_2026_results/gpu_split/manifest_gpu0.tsv \
  --outdir QC_02_2026_results/test_gpu0_qq \
  --threads 1 \
  --mode qq \
  --assessment-mode full \
  --qq-engine gpu \
  --qq-max-points 100000 \
  --limit 4

CUDA_VISIBLE_DEVICES=1 python GWAS_QC/scripts/gwas_mlma_qc_fast.py \
  --manifest QC_02_2026_results/gpu_split/manifest_gpu1.tsv \
  --outdir QC_02_2026_results/test_gpu1_qq \
  --threads 1 \
  --mode qq \
  --assessment-mode full \
  --qq-engine gpu \
  --qq-max-points 100000 \
  --limit 4

Run the full QQ workflow on both GPUs:

mkdir -p QC_02_2026_results/logs

CUDA_VISIBLE_DEVICES=0 python GWAS_QC/scripts/gwas_mlma_qc_fast.py \
  --manifest QC_02_2026_results/gpu_split/manifest_gpu0.tsv \
  --outdir QC_02_2026_results/step2_qq_gpu0 \
  --threads 1 \
  --mode qq \
  --assessment-mode full \
  --qq-engine gpu \
  --qq-max-points 100000 \
  > QC_02_2026_results/logs/step2_qq_gpu0.log 2>&1 &

CUDA_VISIBLE_DEVICES=1 python GWAS_QC/scripts/gwas_mlma_qc_fast.py \
  --manifest QC_02_2026_results/gpu_split/manifest_gpu1.tsv \
  --outdir QC_02_2026_results/step2_qq_gpu1 \
  --threads 1 \
  --mode qq \
  --assessment-mode full \
  --qq-engine gpu \
  --qq-max-points 100000 \
  > QC_02_2026_results/logs/step2_qq_gpu1.log 2>&1 &

Check progress:

jobs
tail -f QC_02_2026_results/logs/step2_qq_gpu0.log
tail -f QC_02_2026_results/logs/step2_qq_gpu1.log

Merge QQ summary tables after both jobs finish:

mkdir -p QC_02_2026_results/step2_qq_gpu_merged

head -n 1 QC_02_2026_results/step2_qq_gpu0/qq_shape_summary.tsv \
  > QC_02_2026_results/step2_qq_gpu_merged/qq_shape_summary.tsv

tail -n +2 QC_02_2026_results/step2_qq_gpu0/qq_shape_summary.tsv \
  >> QC_02_2026_results/step2_qq_gpu_merged/qq_shape_summary.tsv

tail -n +2 QC_02_2026_results/step2_qq_gpu1/qq_shape_summary.tsv \
  >> QC_02_2026_results/step2_qq_gpu_merged/qq_shape_summary.tsv

head -n 1 QC_02_2026_results/step2_qq_gpu0/qc_quality_solutions.tsv \
  > QC_02_2026_results/step2_qq_gpu_merged/qc_quality_solutions.tsv

tail -n +2 QC_02_2026_results/step2_qq_gpu0/qc_quality_solutions.tsv \
  >> QC_02_2026_results/step2_qq_gpu_merged/qc_quality_solutions.tsv

tail -n +2 QC_02_2026_results/step2_qq_gpu1/qc_quality_solutions.tsv \
  >> QC_02_2026_results/step2_qq_gpu_merged/qc_quality_solutions.tsv

Optionally collect QQ images into one directory:

mkdir -p QC_02_2026_results/step2_qq_gpu_merged/qq_plots

cp -r QC_02_2026_results/step2_qq_gpu0/qq_plots/* QC_02_2026_results/step2_qq_gpu_merged/qq_plots/
cp -r QC_02_2026_results/step2_qq_gpu1/qq_plots/* QC_02_2026_results/step2_qq_gpu_merged/qq_plots/

Recommended resource settings:

• metrics: CPU only, --threads 4
• QQ GPU: 2 shell jobs, each with CUDA_VISIBLE_DEVICES=<gpu_id> and --threads 1
• Manhattan: CPU only, --threads 2 to 4
• if GPU memory is limited, keep --qq-max-points 100000

Partial run by population:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_step2_qq_C4 \
  --threads 4 \
  --mode qq \
  --only-pop C4

Partial run for one trait:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_step2_qq_C4_FL1 \
  --threads 1 \
  --mode qq \
  --only-pop C4 \
  --only-trait FL1

Main outputs:

• qq_plots/
• qq_shape_summary.tsv
• qc_quality_solutions.tsv

QQ plot style:

• points are downsampled for drawing using --qq-max-points, but lambda GC and QC metrics use all valid P values
• each plot includes a confidence interval ribbon when SciPy is available
• each plot annotates lambdaGC and valid SNP count

QQ shape classes:

QQ class	Meaning
QQ_IDEAL	QQ body follows the diagonal; acceptable for downstream review
QQ_GOOD_WITH_TAIL_SIGNAL	QQ body is acceptable, with strong tail deviation compatible with true loci
QQ_INFLATED	global upward deviation; likely structure, relatedness, batch effect, or confounding
QQ_DEFLATED	global downward deviation; likely overcorrection or overfitted model
QQ_NOISY	irregular curve; check sample size, missingness, convergence, or low-count variants

Step 3: Manhattan Plots

After reviewing Step 1 and Step 2, draw Manhattan plots. This mode does not redraw QQ plots.

Full run:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_step3_manhattan \
  --threads 4 \
  --mode manhattan \
  --max-plot-points 200000 \
  --plot-keep-p 1e-4

Partial run by population:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_step3_manhattan_C6 \
  --threads 4 \
  --mode manhattan \
  --only-pop C6 \
  --max-plot-points 200000 \
  --plot-keep-p 1e-4

You can also draw one specific population-trait pair:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_step3_manhattan_C6_height \
  --threads 1 \
  --mode manhattan \
  --only-pop C6 \
  --only-trait height \
  --max-plot-points 300000

Main outputs:

• manhattan_plots/
• manhattan_plot_summary.tsv
• qc_summary.tsv
• qc_quality_solutions.tsv

For faster Manhattan plotting on very large .mlma files, the Python script keeps all points with P <= 1e-4 and samples the remaining points when there are too many:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_step3_manhattan_C6_fastplot \
  --threads 4 \
  --mode manhattan \
  --only-pop C6 \
  --max-plot-points 200000 \
  --plot-keep-p 1e-4

Increase --max-plot-points for publication-style figures. Decrease it for faster visual QC. Manhattan colors use a colorblind-safe blue/orange palette, with a dark red genome-wide significance line.

Legacy Combined Plot

The older combined Manhattan + QQ mode is still available if needed:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_legacy_full_C6 \
  --threads 4 \
  --mode full \
  --only-pop C6

R Fallback

The R script supports the same three modes:

Rscript scripts/gwas_mlma_qc.R \
  --manifest manifest.tsv \
  --outdir qc_step1_metrics_R \
  --threads 8 \
  --mode metrics

Run QC Legacy One-Step Example

Rscript scripts/gwas_mlma_qc.R \
  --manifest manifest.tsv \
  --outdir qc_output \
  --threads 8 \
  --mode full

For small test runs:

Rscript examples/create_smoke_data.R smoke_data
Rscript scripts/gwas_mlma_qc.R \
  --manifest smoke_data/manifest.tsv \
  --outdir qc_test \
  --threads 2 \
  --mode metrics

Fast Python smoke test:

python examples/create_smoke_data.py smoke_data
python scripts/gwas_mlma_qc_fast.py \
  --manifest smoke_data/manifest.tsv \
  --outdir qc_test_fast \
  --threads 2 \
  --mode qq

Input Column Detection

The script automatically recognizes common .mlma columns:

Meaning	Accepted names
chromosome	Chr, CHR, chrom, chromosome
position	bp, BP, pos, position
SNP ID	SNP, rs, ID, marker, variant
P value	p, P, pval, PVAL, p_value
effect	b, BETA, beta, effect
standard error	se, SE, stderr
allele frequency	Freq, freq, AF, maf, MAF

If your files use unusual column names, use the explicit options:

python scripts/gwas_mlma_qc_fast.py \
  --manifest manifest.tsv \
  --outdir qc_output \
  --chr-col Chr \
  --bp-col bp \
  --snp-col SNP \
  --p-col p

Default QC Rules

The defaults are intentionally conservative for sample sizes around 200-540.

QC status	Rule
PASS	no major automated problem detected
WARN	borderline lambda GC, few SNPs, many missing P values, or unusual significant count
FAIL	severe lambda GC deviation, invalid P values dominate, missing required columns, unreadable file

Default thresholds:

Parameter	Default
lambda GC fail high	1.20
lambda GC warn high	1.10
lambda GC warn low	0.95
lambda GC fail low	0.90
minimum SNP count fail	10000
minimum SNP count warn	50000
genome-wide significance	5e-8

You can override thresholds:

Rscript scripts/gwas_mlma_qc.R \
  --manifest manifest.tsv \
  --outdir qc_output \
  --mode metrics \
  --lambda-warn-high 1.10 \
  --lambda-fail-high 1.20 \
  --lambda-warn-low 0.95 \
  --lambda-fail-low 0.90 \
  --min-snps-warn 50000 \
  --min-snps-fail 10000

How To Interpret Problem Patterns

Pattern	Common cause	Recommended adjustment
QQ plot globally above diagonal, lambda GC high	population structure, relatedness, batch effect, phenotype confounding	add/top up PCA covariates; check GRM/kinship; include batch/sex/age; rerun by population then meta-analyze
QQ plot below diagonal, lambda GC low	overcorrection, too many covariates, phenotype residualization too strong	reduce redundant covariates; check kinship and residualization strategy
Many significant SNPs across many chromosomes	phenotype or batch artifact	test phenotype against batch/plate/family/sex; inspect outliers; winsorize or transform trait if justified
One chromosome or region has broad elevation	map/build mismatch, local genotyping artifact, structural region, low-quality imputation	check genome build, allele coding, INFO/MAF/missingness; rerun excluding problematic low-quality SNPs
Very few SNPs or many invalid P values	file truncation, ID mismatch, overly strict filtering, model convergence failures	check input genotype filters, sample ID matching, model logs, and per-file line counts
Only small sample population is noisy	low power and unstable estimates	use as sensitivity analysis; prioritize larger groups or meta-analysis

Final Quality Classes

The file qc_quality_solutions.tsv summarizes the combined quality decision and recommended action.

Quality class	Interpretation	Main action
PASS_IDEAL	Numeric QC passes; no major QQ problem detected	Keep for downstream analysis
PASS_WITH_TAIL_SIGNAL	QQ body is good with tail signal	Draw full Manhattan + QQ and inspect whether peaks are localized
WARN_EXCESS_SIGNAL	Too many significant SNPs	Check phenotype-batch associations and outliers
WARN_NOISY_QQ	QQ is irregular	Check sample size, missingness, low MAC variants, and convergence
FAIL_OR_WARN_INFLATION	Lambda/QQ inflation	Add/check PCs, kinship, batch covariates, and population-specific analysis
FAIL_OR_WARN_DEFLATION	Lambda/QQ deflation	Check overcorrection, residualization, and overfitted mixed model
WARN_DATA_COMPLETENESS	Fewer SNPs or more invalid P values than expected	Check filters and logs; may still be usable after review
FAIL_DATA_COMPLETENESS	Too few valid SNPs or many invalid P values	Check file completeness, SNP filters, imputation/MAF thresholds, and model logs
ERROR_INPUT	File or column problem	Check path, delimiter, and required MLMA columns
REVIEW	Borderline mixed signals	Review QQ plus full Manhattan plot

Suggested Review Workflow

1. Run batch QC for all files.
2. Open qc_summary.tsv and sort by status, lambda_gc, n_significant, and n_valid_p.
3. Review only FAIL and WARN plots first.
4. Fix upstream model/filtering problems.
5. Rerun the failed subset using a smaller manifest.
6. Preserve the first and rerun QC outputs for audit.

GitHub Usage

Initial upload from a server:

git init
git add README.md scripts config examples .gitignore
git commit -m "Add batch GWAS MLMA QC workflow"
git branch -M main
git remote add origin git@github.com:YOUR_ORG/gwas-mlma-batch-qc.git
git push -u origin main

After updates:

git add README.md scripts config examples .gitignore
git commit -m "Update GWAS QC workflow"
git push

Notes

• This pipeline performs summary-statistics QC. It does not replace genotype-level QC.
• For publication-grade interpretation, inspect cohort design, kinship, phenotype transformation, covariates, and genotype filtering logs.
• If ancestry or population structure differs across the 3 groups, analyze each group separately and combine with meta-analysis rather than forcing a single pooled model.