Genotype preprocessing

Convert raw genotype VCFs into quality-controlled, sample-matched PLINK datasets and genetic principal components for downstream xQTL analysis.

Miniprotocol Timing¶

This is the total duration for all phases; module-specific timings appear on their respective pages.

Timing: ~3-5 min (on the toy dataset)

Overview¶

This mini-protocol walks through how raw genotype VCFs are normalized, converted to PLINK, quality-controlled, matched to molecular-phenotype samples, and prepared for population-structure adjustment. Each step calls a workflow from VCF_QC.ipynb, genotype_formatting.ipynb, GWAS_QC.ipynb, or PCA.ipynb.

The commands are organized as selectable routes rather than one mandatory 12-step chain. The standard route estimates principal components in unrelated individuals. Use the conditional extension only when related individuals will remain in the analysis and must be projected into the same PCA space before the datasets are recombined.

Steps¶

Choose a route before running commands; the 12 commands are not one mandatory chain.

Analysis goal	Commands to run, in order	Inputs
Produce a QC-passed, chromosome-partitioned PLINK dataset	1 → 2 → 3 → 4	`tests/fixtures/vcf_qc/protocol_example.genotype.chr22.vcf.gz` `input/reference_data/00-All.add_chr.variants.gz` `input/reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta`
Match samples and estimate PCA in unrelated individuals	1 → 2 → 3 → 4 → 5 → 6 → 7 → 8	Basic-route inputs plus `tests/fixtures/gene_annotation/protocol_example.rnaseq.bed.gz`
Retain related individuals by PCA projection	1 → 2 → 3 → 4 → 5 → 6 → 7 → 8 → 9 → 10 → 11 → 12	Standard-PCA inputs plus `tests/fixtures/pca/protocol_example.pca_pheno.txt`

Run only the row matching the intended analysis goal, following its commands in numerical order.

1. Quality-control the input VCF¶

What it does: Normalize and filter the toy VCF against dbSNP and GRCh38.

sos run pipeline/VCF_QC.ipynb qc \
    --genoFile tests/fixtures/vcf_qc/protocol_example.genotype.chr22.vcf.gz \
    --dbsnp-variants <path/to/00-All.add_chr.variants.gz> \
    --reference-genome <path/to/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta> \
    --cwd output/vcf_qc \
    --skip-vcf-header-filtering True

2. Convert the QC-passed VCF to PLINK and merge chromosomes¶

What it does: Convert the output of step 1 to PLINK. The merge command generalizes to multiple chromosome-level files.

sos run pipeline/genotype_formatting.ipynb vcf_to_plink \
    --genoFile output/vcf_qc/protocol_example.genotype.chr22.leftnorm.vcf.gz \
    --cwd output/genotype_formatting/plink \
    --name protocol_example \
    -j 4

sos run pipeline/genotype_formatting.ipynb merge_plink \
    --genoFile `ls output/genotype_formatting/plink/protocol_example.genotype.chr*.bed` \
    --name protocol_example.genotype.merged \
    --cwd output/genotype_formatting/plink \
    -j 2

3. Apply PLINK-level quality control¶

What it does: Apply genotype-, sample- and Hardy-Weinberg-equilibrium filters to the merged PLINK dataset.

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

4. Partition the QC-passed genotype data by chromosome¶

What it does: Create chromosome-specific PLINK files required by chromosome-oriented downstream workflows.

sos run pipeline/genotype_formatting.ipynb genotype_by_chrom \
    --genoFile output/gwas_qc/plink/protocol_example.genotype.merged.plink_qc.bed \
    --cwd output/genotype_by_chrom \
    --chrom `cut -f 1 output/gwas_qc/plink/protocol_example.genotype.merged.plink_qc.bim | uniq | sed "s/chr//g"` \
    -j 4

5. Match genotype and molecular-phenotype samples¶

What it does: Retain the sample intersection between the QC-passed genotype data and molecular phenotype.

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

What it does: Use KING to identify related pairs and produce related and unrelated subsets.

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

7. Prepare the unrelated, LD-pruned PCA subset¶

What it does: Apply the minor-allele-count filter and LD pruning used to estimate ancestry axes.

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

8. Estimate principal components in unrelated individuals¶

What it does: Estimate the PCA model and scores in unrelated individuals.

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 it does: Restrict the related subset to the variants used by the unrelated-sample PCA model.

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

What it does: Project related individuals into the PCA space and identify ancestry-space outliers.

sos run pipeline/PCA.ipynb project_samples \
    --cwd output/pca_uf \
    --genoFile output/pca_related/protocol_example.genotype.merged.plink_qc.protocol_example.king.related.for_pca.plink_qc.extracted.bed \
    --phenoFile tests/fixtures/pca/protocol_example.pca_pheno.txt \
    --pca-model output/pca_uf/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.plink_qc.prune.protocol_example.pca.rds \
    --label-col race --pop-col race --name protocol_example --maha-k 2

11. Remove projected PCA outliers¶

What it does: Remove projected outliers before recombining samples.

sos run pipeline/GWAS_QC.ipynb qc_no_prune \
    --cwd output/pca_related \
    --genoFile output/pca_related/protocol_example.genotype.merged.plink_qc.protocol_example.king.related.for_pca.plink_qc.extracted.bed \
    --remove-samples output/pca_uf/protocol_example.pca_pheno.pca.projected.outliers \
    --name no_outlier

What it does: Merge the unrelated PCA subset with the retained projected related samples for downstream analysis.

sos run pipeline/genotype_formatting.ipynb merge_plink \
    --genoFile output/gwas_qc/genotype/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.plink_qc.prune.bed \
               output/pca_related/protocol_example.genotype.merged.plink_qc.protocol_example.king.related.for_pca.plink_qc.extracted.no_outlier.plink_qc.bed \
    --cwd output/genotype_final \
    --name protocol_example.qced

Output¶

Route or step	Products
Step 1: VCF quality control	`output/vcf_qc/*.leftnorm.vcf.gz` and its index
Step 2: PLINK conversion and merge	`output/genotype_formatting/plink/protocol_example.genotype.merged.{bed,bim,fam}`
Step 3: PLINK quality control	`output/gwas_qc/plink/*.plink_qc.{bed,bim,fam}`
Step 4: Chromosome partitioning	`output/genotype_by_chrom/*.{bed,bim,fam}`
Step 5: Sample matching	`output/gwas_qc/genotype/*.sample_genotypes.txt`
Step 6: Kinship	KING results and related/unrelated PLINK subsets under `output/gwas_qc/kinship/`
Steps 7–8: Unrelated-sample PCA	LD-pruned PLINK data under `output/gwas_qc/genotype/`, plus the PCA model, scores and plots under `output/pca_uf/`
Steps 9–11: Related-sample projection	Variant-matched related PLINK data under `output/pca_related/`, projected PCA scores and the outlier list under `output/pca_uf/`, and an outlier-filtered related subset
Step 12: Recombined analysis dataset	`output/genotype_final/protocol_example.qced.{bed,bim,fam}`

Anticipated Results¶

The QC route produces normalized genotype data in PLINK format, including chromosome-partitioned files for downstream xQTL workflows. The standard PCA route additionally produces sample-matched related and unrelated subsets, an LD-pruned unrelated dataset, and principal-component scores estimated without close relatives.

When related individuals are retained, steps 9–12 project them using the unrelated-sample PCA loadings, remove designated projection outliers, and recombine the retained samples into the final analysis dataset.

Proceed to phenotype and covariate preprocessing before association testing.

Command Interface¶

List the workflows and parameters available in each module used by this mini-protocol.

sos run pipeline/VCF_QC.ipynb -h
sos run pipeline/genotype_formatting.ipynb -h
sos run pipeline/GWAS_QC.ipynb -h
sos run pipeline/PCA.ipynb -h