- **PLINK 1.9** (`plink`) and **PLINK 2** (`plink2`) — install and add to your `PATH` ([PLINK 1.9 & 2 installation](#plink-192-installation)).
- **Tutorial genotype data (PLINK binary)** — ~1M variants, 504 EAS samples (1000 Genomes Phase 3); download from [01_Dataset](../01_Dataset/README.md) with `download_sampledata.sh` ([Download genotype data](#download-genotype-data)).

**Performance & File Formats:**
- PLINK 2 is significantly faster and uses less memory, especially for large datasets
- PLINK 2 supports new file formats (`.pgen`, `.pvar`, `.psam`) in addition to traditional `.bed/.bim/.fam` formats
- PLINK 2 can handle multiallelic variants natively

**Allele Concepts:**
- PLINK 1.9 uses **Major/Minor** alleles (A1/A2) - determined by frequency in your dataset
- PLINK 2 uses **Reference/Alternative** alleles (REF/ALT) - determined by the reference genome
- This means MAF in PLINK 1.9 ranges [0, 0.5], while ALT frequency in PLINK 2 ranges [0, 1]

**Command Syntax:**
- Most basic commands are similar, but PLINK 2 has updated syntax for some operations
- PLINK 2 uses `--bfile` for binary files (same as PLINK 1.9) but also supports `--pfile` for pgen format

**Features:**
- PLINK 2 includes KING-robust kinship estimator (more robust for population structure)
- Some advanced features may differ between versions
- PLINK 1.9 is still widely used and well-documented for standard QC workflows

**When to use which:**
- Use PLINK 1.9 for compatibility with existing workflows and extensive documentation
- Use PLINK 2 for faster processing of large datasets and when working with multiallelic variants

```bash
export PATH=$PATH:~/tools/bin/
```
This command will add the path to your current shell. 

If you restart the terminal, it will be lost. So you may need to add it to the Bash  configuration file. Then run 

```
echo "export PATH=$PATH:~/tools/bin/" >> ~/.bashrc
```

This will add a new line at the end of `.bashrc`, which will be run every time you open a new bash shell.

```bash
plink
PLINK v1.90b7.2 64-bit (11 Dec 2023)           www.cog-genomics.org/plink/1.9/
(C) 2005-2023 Shaun Purcell, Christopher Chang   GNU General Public License v3

plink <input flag(s)...> [command flag(s)...] [other flag(s)...]
plink --help [flag name(s)...]

Commands include --make-bed, --recode, --flip-scan, --merge-list,
--write-snplist, --list-duplicate-vars, --freqx, --missing, --test-mishap,
--hardy, --mendel, --ibc, --impute-sex, --indep-pairphase, --r2, --show-tags,
--blocks, --distance, --genome, --homozyg, --make-rel, --make-grm-gz,
--rel-cutoff, --cluster, --pca, --neighbour, --ibs-test, --regress-distance,
--model, --bd, --gxe, --logistic, --dosage, --lasso, --test-missing,
--make-perm-pheno, --tdt, --qfam, --annotate, --clump, --gene-report,
--meta-analysis, --epistasis, --fast-epistasis, and --score.

"plink --help | more" describes all functions (warning: long).
```

```bash
plink2
PLINK v2.00a5.9LM AVX2 AMD (12 Dec 2023)       www.cog-genomics.org/plink/2.0/
(C) 2005-2023 Shaun Purcell, Christopher Chang   GNU General Public License v3

plink2 <input flag(s)...> [command flag(s)...] [other flag(s)...]
plink2 --help [flag name(s)...]

Commands include --rm-dup list, --make-bpgen, --export, --freq, --geno-counts,
--sample-counts, --missing, --hardy, --het, --fst, --indep-pairwise, --ld,
--sample-diff, --make-king, --king-cutoff, --pmerge, --pgen-diff,
--write-samples, --write-snplist, --make-grm-list, --pca, --glm, --adjust-file,
--gwas-ssf, --clump, --score, --variant-score, --genotyping-rate, --pgen-info,
--validate, and --zst-decompress.

"plink2 --help | more" describes all functions.
```

```bash
cd ../01_Dataset
./download_sampledata.sh
```

And you will get the following three PLINK files:

```
-rw-r--r-- 1 yunye yunye 149M Dec 26 13:25 1KG.EAS.auto.snp.norm.nodup.split.rare002.common015.missing.bed
-rw-r--r-- 1 yunye yunye  40M Dec 26 13:25 1KG.EAS.auto.snp.norm.nodup.split.rare002.common015.missing.bim
-rw-r--r-- 1 yunye yunye  13K Dec 26 13:25 1KG.EAS.auto.snp.norm.nodup.split.rare002.common015.missing.fam
```

Check the bim file:

```bash
head 1KG.EAS.auto.snp.norm.nodup.split.rare002.common015.missing.bim
1       1:14930:A:G     0       14930   G       A
1       1:15774:G:A     0       15774   A       G
1       1:15777:A:G     0       15777   G       A
1       1:57292:C:T     0       57292   T       C
1       1:77874:G:A     0       77874   A       G
1       1:87360:C:T     0       87360   T       C
1       1:92917:T:A     0       92917   A       T
1       1:104186:T:C    0       104186  T       C
1       1:125271:C:T    0       125271  C       T
1       1:232449:G:A    0       232449  A       G
```

Check the fam file:
```bash
head 1KG.EAS.auto.snp.norm.nodup.split.rare002.common015.missing.fam
HG00403 HG00403 0 0 0 -9
HG00404 HG00404 0 0 0 -9
HG00406 HG00406 0 0 0 -9
HG00407 HG00407 0 0 0 -9
HG00409 HG00409 0 0 0 -9
HG00410 HG00410 0 0 0 -9
HG00419 HG00419 0 0 0 -9
HG00421 HG00421 0 0 0 -9
HG00422 HG00422 0 0 0 -9
HG00428 HG00428 0 0 0 -9
```

Suppose we have N samples and M SNPs for each sample.

For sample $j$ :

$$Sample\ Missing\ Rate_{j} = {{N_{missing\ SNPs\ for\ j}}\over{M}} = 1 - Call\ Rate_{sample, j}$$

For SNP $i$ :

$$SNP\ Missing\ Rate_{i} = {{N_{missing\ samples\ at\ i}}\over{N}} = 1 - Call\ Rate_{SNP, i}$$

```bash
cd ../04_Data_QC
genotypeFile="../01_Dataset/1KG.EAS.auto.snp.norm.nodup.split.rare002.common015.missing" #!!! Please add your own path here.  "1KG.EAS.auto.snp.norm.nodup.split.rare002.common015.missing" is the prefix of PLINK bed file. 

plink \
	--bfile ${genotypeFile} \
	--missing \
	--out plink_results
```
Remember to set the value for `${genotypeFile}`.

Note: The missing values were simulated based on normal distributions for each individual.  

Interactive plots (pan/zoom, hover counts): sample missing rate (top) and SNP missing rate (bottom). Regenerate with `python3 plot_missing_rate_interactive.py` in this folder (after `plink --missing`; the same script also builds the [F_het figure](#inbreeding-f-coefficient) and the [LD r² heatmap](#ld-calculation) when the corresponding PLINK outputs exist).

<iframe src="../assets/plots/04_data_qc/missing_rate_distributions.html" width="100%" height="740" style="border:0;" title="Missing rate distributions"></iframe>

Suppose the reference allele(REF) is A and the alternative allele(ALT) is B for a certain SNP. The possible genotypes are AA, AB and BB. In a population of N samples (2N alleles), $N = N_{AA} + N_{AB} + N_{BB}$:

- the number of A alleles:  $N_A = 2 \times N_{AA} + N_{AB}$
- the number of B alleles:  $N_B = 2 \times N_{BB} + N_{AB}$

So we can calculate the allele frequency:

- Reference Allele Frequency : $AF_{REF}= {{N_A}\over{N_A + N_B}}$
- Alternative Allele Frequency : $AF_{ALT}= {{N_B}\over{N_A + N_B}}$

The MAF for this SNP in this specific population is defined as:

$MAF = min( AF_{REF}, AF_{ALT} )$

```bash
# results from plink1.9
head plink_results.frq
CHR              SNP   A1   A2          MAF  NCHROBS
1      1:14930:A:G    G    A       0.4133     1004
1      1:15774:G:A    A    G      0.02794     1002
1      1:15777:A:G    G    A      0.07385     1002
1      1:57292:C:T    T    C       0.1054      996
1      1:77874:G:A    A    G      0.01996     1002
1      1:87360:C:T    T    C      0.02286     1006
1      1:92917:T:A    A    T     0.003018      994
1     1:104186:T:C    T    C        0.499     1002
1     1:125271:C:T    C    T      0.03088     1004
```

```bash
# results from plink2
head plink_results.afreq
#CHROM  ID      REF     ALT     PROVISIONAL_REF?        ALT_FREQS       OBS_CT
1       1:14930:A:G     A       G       Y       0.413347        1004
1       1:15774:G:A     G       A       Y       0.0279441       1002
1       1:15777:A:G     A       G       Y       0.0738523       1002
1       1:57292:C:T     C       T       Y       0.105422        996
1       1:77874:G:A     G       A       Y       0.0199601       1002
1       1:87360:C:T     C       T       Y       0.0228628       1006
1       1:92917:T:A     T       A       Y       0.00301811      994
1       1:104186:T:C    T       C       Y       0.500998        1002
1       1:125271:C:T    C       T       Y       0.969124        1004
```

$$P(N_{AB} = n_{AB} | N, n_A) = {{2^{n_{AB}}}N!\over{n_{AA}!n_{AB}!n_{BB}!}} \times {{n_A!n_B!}\over{(2N)!}} $$

To compute the Hardy-Weinberg equilibrium exact test statistics, we will sum up the probabilities of all configurations with probability equal to or less than the observed configuration :

$$P_{HWE} = \sum_{n^{*}_AB} I[P(N_{AB} = n_{AB} | N, n_A) \geqq P(N_{AB} = n^{*}_{AB} | N, n_A)] \times P(N_{AB} = n^{*}_{AB} | N, n_A)$$

$I(x)$ is the indicator function. If x is true, $I(x) = 1$; otherwise, $I(x) = 0$.

Reference : Wigginton, J. E., Cutler, D. J., & Abecasis, G. R. (2005). A note on exact tests of Hardy-Weinberg equilibrium. The American Journal of Human Genetics, 76(5), 887-893. [Link](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1199378/)

This code is converted from [here](https://github.com/jeremymcrae/snphwe/blob/master/src/snp_hwe.cpp) (Jeremy McRae) to python. Orginal citation: Wigginton, JE, Cutler, DJ, and Abecasis, GR (2005) A Note on Exact Tests of Hardy-Weinberg Equilibrium. AJHG 76: 887-893
```
def snphwe(obs_hets: int, obs_hom1: int, obs_hom2: int) -> float:
"""Calculate Hardy-Weinberg equilibrium exact test for a variant

Args:
    obs_hets (int): observed heterozygous number
    obs_hom1 (int): first observed homozygous number
    obs_hom2 (int): second observed homozygous number

Returns:
    float: Hardy-Weinberg equilibrium exact test p-value.
    
    - P > 0.05: No significant deviation from Hardy-Weinberg equilibrium (HWE).
    - P ≤ 0.05: Significant deviation from HWE, which could be due to factors such as population stratification, inbreeding, genotyping errors, or natural selection.
"""
obs_minor_homs = min(obs_hom1, obs_hom2)
obs_major_homs = max(obs_hom1, obs_hom2)

rare = 2 * obs_minor_homs + obs_hets
genotypes = obs_hets + obs_major_homs + obs_minor_homs

probs = [0.0 for i in range(rare +1)]

mid = rare * (2 * genotypes - rare) // (2 * genotypes)
if mid % 2 != rare%2:
    mid += 1

probs[mid] = 1.0
sum_p = 1 #probs[mid]

curr_homr = (rare - mid) // 2
curr_homc = genotypes - mid - curr_homr

for curr_hets in range(mid, 1, -2):
    probs[curr_hets - 2] = probs[curr_hets] * curr_hets * (curr_hets - 1.0)/ (4.0 * (curr_homr + 1.0) * (curr_homc + 1.0))
    sum_p+= probs[curr_hets - 2]
    curr_homr += 1
    curr_homc += 1

curr_homr = (rare - mid) // 2
curr_homc = genotypes - mid - curr_homr

for curr_hets in range(mid, rare-1, 2):
    probs[curr_hets + 2] = probs[curr_hets] * 4.0 * curr_homr * curr_homc/ ((curr_hets + 2.0) * (curr_hets + 1.0))
    sum_p += probs[curr_hets + 2]
    curr_homr -= 1
    curr_homc -= 1

target = probs[obs_hets]
p_hwe = 0.0
for p in probs:
    if p <= target :
        p_hwe += p / sum_p  

return min(p_hwe,1)

if __name__ == '__main__':
# For an example, we had 502 samples. 
# At position 1:14930:A:G, we count:
#     + 407 heterozygous genotypes (signed 0|1 or 1|0).
#     + 4 homozygous genotypes AA (signed 0|0)
#     + 91 homozygous genotypes GG (signed 1|1)
print(snphwe(407,4,91))
```

- consider a window of `x` SNPs
- calculate LD between each pair of SNPs in the window and remove one of a pair of SNPs if the LD is greater than `z`
- shift the window `y` SNPs forward and repeat the procedure.

Please check https://www.cog-genomics.org/plink/1.9/ld#indep for details.

You can check the PLINK log for how many variants were removed based on the filters you applied:
```
Total genotyping rate in remaining samples is 0.993916.
108837 variants removed due to missing genotype data (--geno).
--hwe: 9754 variants removed due to Hardy-Weinberg exact test.
87149 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
1029376 variants and 501 people pass filters and QC.
```

Let's take a look at the LD-pruned SNP file. Basically, it just contains one SNP id per line.

```bash
head plink_results.prune.in
1:15774:G:A
1:15777:A:G
1:77874:G:A
1:87360:C:T
1:125271:C:T
1:232449:G:A
1:533113:A:G
1:565697:A:G
1:566933:A:G
1:567092:T:C
```

$$F = {{O(HOM) - E(HOM)}\over{ M - E(HOM)}}$$

- $E(HOM)$ :Expected Homozygous Genotype Count 
- $O(HOM)$ :Observed Homozygous Genotype Count 
- M : Number of SNPs

High F may indicate a relatively high level of inbreeding. 

Low F may suggest the sample DNA was contaminated.

```bash
plink \
	--bfile ${genotypeFile} \
    --extract plink_results.prune.in \
	--het \
	--out plink_results
```

Check the output:

```bash
head plink_results.het
    FID       IID       O(HOM)       E(HOM)        N(NM)            F
HG00403   HG00403       180222    1.796e+05       217363      0.01698
HG00404   HG00404       180127    1.797e+05       217553      0.01023
HG00406   HG00406       178891    1.789e+05       216533   -0.0001138
HG00407   HG00407       178992     1.79e+05       216677   -0.0008034
HG00409   HG00409       179918    1.801e+05       218045    -0.006049
HG00410   HG00410       179782    1.801e+05       218028    -0.009268
HG00419   HG00419       178362    1.783e+05       215849     0.001315
HG00421   HG00421       178222    1.785e+05       216110    -0.008288
HG00422   HG00422       178316    1.784e+05       215938      -0.0022
```

Interactive histogram (pan/zoom, hover counts). Dashed lines: ±0.1 threshold used below. Regenerate with `python3 plot_missing_rate_interactive.py` in this folder (after `plink --het`).

<iframe src="../assets/plots/04_data_qc/f_het_distribution.html" width="100%" height="540" style="border:0;" title="F heterozygosity distribution"></iframe>

*Identity by state* (IBS) means two alleles *look the same* at a marker (same nucleotide), which can happen either because they were inherited from a common ancestor or simply because that allele is frequent in the population.     *Identity by descent* (IBD) means the alleles are *the same physical copy* from a shared ancestor—genealogical relatedness.

These differ in **who** is involved and **what** is randomized:

| Quantity | Symbol | Who | What it measures |
|----------|--------|-----|------------------|
| **Coefficient of kinship** (Malécot) | φ | **Two** people *A*, *B* | Probability that one random allele from *A* and one random allele from *B* (same autosomal locus) are **IBD**. The **KING kinship** column in the next table uses this scale (duplicates **≈ 0.5**). |
| **Coefficient of relationship** | *r* | **Two** people | **Fraction** of segregating genetic material in one person expected to be IBD with the other. For **outbred, non-inbred** pairs, **r = 2φ** (e.g. PO / full sibs **r ≈ 0.5**, **φ ≈ 0.25**; first cousins **r ≈ 0.125**, **φ ≈ 0.0625**). |
| **PI_HAT** (PLINK 1 `--genome`) | *π̂* | **Two** people | **P(Z1)/2 + P(Z2)** from `--genome`; **PI_HAT ≈ *r*** on ideal outbred pedigrees—same as the **PI_HAT** column below. |
| **Coefficient of inbreeding** | *f*, *F* | **One** person | Probability that its **two alleles at a locus** are IBD **with each other** (homozygosity by descent through the pedigree). **f = 0** if parents are unrelated and non-inbred; consanguinity raises *f* (e.g. child of first cousins often **f ≈ 1/16**). Alters *r* / φ for some pairs; many GWAS pipelines assume **f ≈ 0**. |

Ideal **outbred** pairs; real GWAS estimates scatter around these values.

| Relationship | Degree | PI_HAT | KING kinship |
|--------------|--------|--------|--------------|
| Unrelated | — | 0 | 0 |
| Parent–offspring | 1st | 0.5 | 0.25 |
| Full siblings | 1st | 0.5 | 0.25 |
| Half-siblings | 2nd | 0.25 | 0.125 |
| Grandparent–grandchild | 2nd | 0.25 | 0.125 |
| Uncle/aunt–niece/nephew | 2nd | 0.25 | 0.125 |
| First cousins | 3rd | 0.125 | 0.0625 |
| Identical (MZ) twins (or duplicate DNA sample) | — | 1 | 0.5 |

The prior probability of IBS sharing can be modeled as: 

$$P(I=i) = \sum^{z=i}_{z=0}P(I=i|Z=z)P(Z=z)$$

- I: IBS state (I = 0, 1, or 2)
- Z: IBD state (Z = 0, 1, or 2)
- $P(I|Z)$ is a function of allele frequency. PLINK will average over all SNPs to obtain the expected value for $P(I|Z)$.

So the proportion of alleles shared IBD ($\hat{\pi}$) can be estimated by:

$$\hat{\pi} = {{P(Z=1)}\over{2}} + P(Z=2)$$

For full column definitions, see the PLINK IBD documentation: https://www.cog-genomics.org/plink/1.9/ibd
```bash
head plink_results.genome
    FID1     IID1     FID2     IID2 RT    EZ      Z0      Z1      Z2  PI_HAT PHE       DST     PPC   RATIO
HG00403  HG00403  HG00404  HG00404 UN    NA  1.0000  0.0000  0.0000  0.0000  -1  0.858562  0.3679  1.9774
HG00403  HG00403  HG00406  HG00406 UN    NA  0.9805  0.0044  0.0151  0.0173  -1  0.858324  0.8183  2.0625
HG00403  HG00403  HG00407  HG00407 UN    NA  0.9790  0.0000  0.0210  0.0210  -1  0.857794  0.8034  2.0587
HG00403  HG00403  HG00409  HG00409 UN    NA  0.9912  0.0000  0.0088  0.0088  -1  0.857024  0.2637  1.9578
HG00403  HG00403  HG00410  HG00410 UN    NA  0.9699  0.0235  0.0066  0.0184  -1  0.858194  0.6889  2.0335
HG00403  HG00403  HG00419  HG00419 UN    NA  1.0000  0.0000  0.0000  0.0000  -1  0.857643  0.8597  2.0745
HG00403  HG00403  HG00421  HG00421 UN    NA  0.9773  0.0218  0.0010  0.0118  -1  0.857276  0.2186  1.9484
HG00403  HG00403  HG00422  HG00422 UN    NA  0.9880  0.0000  0.0120  0.0120  -1  0.857224  0.8277  2.0652
HG00403  HG00403  HG00428  HG00428 UN    NA  0.9801  0.0069  0.0130  0.0164  -1  0.858162  0.9812  2.1471
```

Manichaikul, A., Mychaleckyj, J. C., Rich, S. S., Daly, K., Sale, M., & Chen, W. M. (2010). Robust relationship inference in genome-wide association studies. Bioinformatics, 26(22), 2867-2873.

Pick a maximum allowed kinship, then **prune** so **no retained pair** is above it (**`--king-cutoff`** drops one sample per offending pair). Each value below is ≈ √(product of two adjacent expected kinships); **lower** *t* = **stricter**. Noise and population structure still matter.

| Prune if KING > … | √(product) | Typical meaning: exclude pairs at least this closely related | ≈ PI_HAT* |
|-------------------|------------|----------------------------------------------------------------|-----------|
| ~0.354 | 0.5 × 0.25 | Duplicate / MZ | ~0.708 |
| ~0.177 | 0.25 × 0.125 | Within 1st degree (and closer) | ~0.354 |
| ~0.0884 | 0.125 × 0.0625 | Within 2nd degree (and closer) | ~0.177 |
| ~0.0442 | 0.0625 × 0.03125 | Within 3rd degree (and closer, e.g. first cousins) | ~0.088 |

\*Illustrative **2 × KING** on pedigree expectations only; `--genome` vs KING differ in practice.

```bash
plink \
	    --bfile ${genotypeFile} \
        --chr 22 \
        --r2 \
        --out plink_results
```

```bash
head plink_results.ld
 CHR_A         BP_A             SNP_A  CHR_B         BP_B             SNP_B           R2
22     16069141   22:16069141:C:G     22     16071624   22:16071624:A:G     0.771226
22     16069784   22:16069784:A:T     22     16149743   22:16149743:T:A     0.217197
22     16069784   22:16069784:A:T     22     16150589   22:16150589:C:A     0.224992
22     16069784   22:16069784:A:T     22     16159060   22:16159060:G:A       0.2289
22     16149743   22:16149743:T:A     22     16150589   22:16150589:C:A     0.965109
22     16149743   22:16149743:T:A     22     16152606   22:16152606:T:C     0.692157
22     16149743   22:16149743:T:A     22     16159060   22:16159060:G:A     0.721796
22     16149743   22:16149743:T:A     22     16193549   22:16193549:C:T     0.336477
22     16149743   22:16149743:T:A     22     16212542   22:16212542:C:T     0.442424
```

From `plink_results.ld`: pairwise r² for the 10 variants with the smallest base positions among SNPs appearing in the file (ordered by chromosome, then BP). Cells with no value are pairs not reported in this `.ld` output for the chosen SNP set. Regenerate with `python3 plot_missing_rate_interactive.py` when `plink_results.ld` is present.

<iframe src="../assets/plots/04_data_qc/ld_r2_first10_heatmap.html" width="100%" height="660" style="border:0;" title="LD r² heatmap (first 10 SNPs)"></iframe>

```bash
#extract the 1000 samples with the pruned SNPs, and make a bed file.
plink \
	--bfile ${genotypeFile} \
	--extract plink_results.prune.in \
	--make-bed \
	--out plink_1000_pruned

#convert the bed/bim/fam to ped/map
plink \
        --bfile plink_1000_pruned \
        --recode \
        --out plink_1000_pruned
```

1224104 variants and 500 people pass filters and QC.

-rw-r--r--  1 yunye yunye 146M Dec 26 15:40 sample_data.clean.bed
-rw-r--r--  1 yunye yunye  39M Dec 26 15:40 sample_data.clean.bim
-rw-r--r--  1 yunye yunye  13K Dec 26 15:40 sample_data.clean.fam