A Worked Sample-and-Variant QC Pipeline

Why Work Through a QC Pipeline?

Genotyping quality control is not one calculation or one universal cutoff. It is a sequence of checks and decisions applied to data whose structure and intended use must first be understood. A filter can remove unreliable samples or variants, but applying it without recording its denominator, threshold, and order makes the result difficult to reproduce and can conceal why data were excluded.

This tutorial follows a small, hypothetical diploid autosomal dataset from initial file checks through sample and variant call-rate filters. The example uses simple missingness thresholds so every decision can be checked by hand. In practice, a pipeline may also include genotype-level quality filters and other diagnostics; their results should be reviewed and documented rather than silently folded into a call-rate calculation.

Start With File and Dataset Checks

Before calculating call rates, confirm that the files describe the same data and that the matrix can be interpreted as intended. Checks should include whether sample identifiers are unique and align across genotype, phenotype, and sample-information files; whether variant identifiers and genomic positions are plausible and consistently represented; and whether allele labels and genotype encodings are valid for the stated genome build and ploidy. Confirm that missing genotypes use a recognized representation, not a value that might be mistaken for a valid genotype.

When data have been converted between formats, compare basic counts and sample and variant identifiers before and after conversion. A successful conversion command does not by itself establish that allele coding, missing values, or sample order were preserved correctly. The earlier material on PLINK files and conversion checks is relevant here: document any conversion and verify the resulting data before filtering.

For the worked pipeline, assume a hypothetical dataset of six samples and five autosomal variants. Each genotype is either a diploid call with alleles A and G or missing. There are no duplicate sample or variant identifiers, no unexpected genotype codes, and the samples and variants agree across the input files. Those checks are assumed to have passed; the table gives the genotype matrix entering the call-rate filters.

For a dataset with \(n\) samples and \(m\) variants, let \(G_{ij}\) denote the genotype for sample \(i\) at variant \(j\). Define \(I(G_{ij}\text{ is called})\) to be 1 for a nonmissing genotype and 0 for a missing genotype. Then the sample and variant call rates are

These are descriptive proportions of observed calls in the eligible matrix. They do not establish why a genotype is missing, and they do not by themselves show that a called genotype is correct.

Apply the Sample Filter First

For this worked pipeline, retain samples with a call rate of at least 60%. This is a hypothetical teaching threshold, not a universal recommendation. The threshold and rationale for a real analysis should be selected for the data and study, as discussed in the earlier material on choosing QC thresholds. Here, the initial eligible variants are v1 through v5, so each sample has five opportunities for a call.

Record the actual decision, rather than only noting that a sample filter was run. In this example the sample rule is “retain call rate at least 60%,” S6 is the only failed sample, and S1–S5 remain. The full genotype matrix is then restricted to those five samples before variant call rates are evaluated.

Recalculate Variant Missingness on Retained Samples

After sample filtering, use the samples that remain as the denominator for each variant. For this example, retain variants with a call rate of at least 80%, equivalent to missingness of at most 20%. Again, this is a hypothetical threshold used to make the decisions transparent, not a default for every study.

The order of these steps matters because call rates depend on who or what is eligible for the denominator. For example, v4 had calls in S1–S4 and a missing genotype in S5, as well as a missing genotype in S6. Before sample filtering, its missingness was 2/6, or 33.3%, which would fail a 20% maximum. After S6 is removed, v4 has one missing genotype among five samples, or 20%, and passes. The final decision follows the specified pipeline: remove failed samples first, then calculate variant missingness among the retained samples.

Why the Denominator and Order Must Be Documented

The call rate is a ratio, so changing the eligible set can change the ratio even when the calls for a particular retained sample or variant have not changed. This follows directly from its definition. If \(C_j\) is the number of called genotypes at variant \(j\), and \(N\) is the number of eligible samples, then its call rate is \(C_j/N\). Removing a sample with a missing genotype at that variant leaves \(C_j\) unchanged but reduces \(N\); the new rate is \(C_j/(N-1)\), which is larger whenever \(C_j>0\). If the removed sample had a called genotype, both the numerator and denominator decrease, so the rate may increase or decrease. Therefore, the denominator cannot be inferred just from the variant name or final percentage.

Likewise, applying a variant filter first can change sample call rates, because removing variants changes the number of opportunities for a call in each sample. A pipeline therefore needs a declared order and a clear statement about whether rates are recalculated after each filtering step. A different defensible order may produce different retained sets; what matters is that the procedure is appropriate, consistently implemented, and recorded.

Compare the Starting and Retained Data

The final comparison should distinguish the number of samples and variants from the number of genotype cells and observed calls. Before filtering there are 6 × 5 = 30 genotype cells. The matrix has 9 missing genotypes, so there are 21 called genotypes and an overall call rate of 21/30 = 70.0%. This overall rate is simply a descriptive summary; it does not replace the separate sample and variant checks.

For a real pipeline, preserve enough information to reproduce this summary: the input files and any conversion checks, the order of operations, the eligible denominators, thresholds, counts excluded at each step, and the identifiers of excluded samples and variants. If another filter is added later, state whether its calls are converted to missing before call rates are calculated or whether it is applied at a separate stage. Those choices affect the matrix being evaluated.

Common Mistakes and Their Corrections

• Using the wrong denominator. A variant’s missingness after sample filtering is calculated among retained samples, not the original sample count. State which set is eligible at every step.
• Filtering at the boundary incorrectly. Under a rule requiring a call rate of at least 60%, a sample at exactly 60% passes. Write the inequality explicitly before classifying records.
• Assuming a cutoff is universal. The thresholds in the example are illustrative. Select thresholds for the dataset and analysis, and document the justification instead of treating a familiar number as automatically correct.
• Treating a higher retained call rate as proof of accurate genotypes. A call-rate filter addresses missingness, not all possible genotype errors. Other appropriate QC checks still matter.
• Reporting only a final percentage. A percentage without retained counts, denominators, and filtering order cannot show how the final matrix was produced.