bio-clinical-databases-polygenic-risk

## Version Compatibility

Reference examples tested with: LDpred2 1.14+, PRSice-2 2.3+, numpy 1.26+, scipy 1.12+

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show <package>` then `help(module.function)` to check signatures
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters
- CLI: `<tool> --version` then `<tool> --help` to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# Polygenic Risk Scores

**"Calculate polygenic risk scores for my cohort"** → Compute genome-wide risk scores from GWAS summary statistics and individual genotypes to predict disease susceptibility.
- CLI: `PRSice_linux --base gwas.txt --target genotypes --out prs_results`
- R: `bigsnpr::snp_ldpred2_auto()` for LDpred2 Bayesian PRS

## PRSice-2 Workflow

**Goal:** Calculate polygenic risk scores from GWAS summary statistics using clumping and thresholding.

**Approach:** Run PRSice-2 with GWAS summary stats and target genotypes, applying LD clumping and multiple p-value thresholds.

### Basic PRS Calculation

```bash
# PRSice-2 with clumping and thresholding
PRSice_linux \
    --base gwas_summary.txt \
    --target genotypes \
    --snp SNP \
    --chr CHR \
    --bp BP \
    --A1 A1 \
    --A2 A2 \
    --pvalue P \
    --beta BETA \
    --clump-kb 250 \
    --clump-r2 0.1 \
    --bar-levels 5e-8,1e-5,1e-3,0.01,0.05,0.1,0.5,1 \
    --fastscore \
    --all-score \
    --out prs_results
```

### PRSice-2 with Covariates

```bash
PRSice_linux \
    --base gwas_summary.txt \
    --target genotypes \
    --pheno phenotype.txt \
    --cov covariates.txt \
    --cov-col @PC[1-10],Age,Sex \
    --binary-target T \
    --clump-kb 250 \
    --clump-r2 0.1 \
    --out prs_with_cov
```

## GWAS Summary Statistics Format

```
SNP          CHR  BP        A1  A2  BETA    SE      P
rs12345      1    10000     A   G   0.05    0.01    1e-8
rs67890      1    20000     T   C   -0.03   0.02    0.001
```

## LDpred2 (R)

**Goal:** Compute Bayesian polygenic risk scores with automatic hyperparameter tuning via LDpred2-auto.

**Approach:** Load genotypes with bigsnpr, match GWAS variants, compute LD matrix, estimate heritability with LD score regression, then run LDpred2-auto.

### Setup and Run

```r
library(bigsnpr)
library(data.table)

# Load genotype data (plink bed/bim/fam)
obj.bigsnp <- snp_attach('genotypes.rds')
G <- obj.bigsnp$genotypes
map <- obj.bigsnp$map

# Load and format GWAS summary stats
sumstats <- fread('gwas_summary.txt')

# Match variants
df_beta <- snp_match(sumstats, map, strand_flip = TRUE)

# Compute LD matrix (correlation)
# Uses reference panel or in-sample LD
corr <- snp_cor(G, ind.col = df_beta$`_NUM_ID_`)

# LDpred2-auto (recommended - automatic hyperparameter tuning)
ldsc <- snp_ldsc2(corr, df_beta)
h2_est <- ldsc[['h2']]

multi_auto <- snp_ldpred2_auto(
    corr,
    df_beta,
    h2_init = h2_est,
    vec_p_init = seq_log(1e-4, 0.2, 30),
    ncores = 4
)

# Extract posterior effect sizes
beta_auto <- sapply(multi_auto, function(x) x$beta_est)
pred_auto <- big_prodMat(G, beta_auto)
```

### LDpred2 Grid Model

```r
# Grid of hyperparameters
h2_seq <- round(h2_est * c(0.7, 1, 1.4), 4)
p_seq <- signif(seq_log(1e-5, 1, 21), 2)
params <- expand.grid(p = p_seq, h2 = h2_seq, sparse = c(FALSE, TRUE))

# Run LDpred2-grid
beta_grid <- snp_ldpred2_grid(corr, df_beta, params, ncores = 4)
pred_grid <- big_prodMat(G, beta_grid)

# Select best parameters by validation R2
auc_grid <- apply(pred_grid, 2, function(x) {
    AUC(x, obj.bigsnp$fam$affection - 1)
})
best_params <- params[which.max(auc_grid), ]
```

## PRS-CS

**Goal:** Compute PRS using continuous shrinkage priors with an external LD reference panel.

**Approach:** Run PRS-CS to estimate posterior effect sizes, then score with plink.

```bash
# PRS-CS with external LD reference
python PRScs.py \
    --ref_dir=ldblk_1kg_eur \
    --bim_prefix=target \
    --sst_file=gwas_summary.txt \
    --n_gwas=100000 \
    --out_dir=prscs_output

# Score with plink
plink --bfile target \
    --score prscs_output_pst_eff_a1_b0.5_phi1e-02.txt 2 4 6 \
    --out prs_scores
```

## Score Normalization

**Goal:** Normalize raw PRS values to Z-scores and population percentiles for interpretable reporting.

**Approach:** Z-score normalize against a reference distribution, then convert to percentiles via the normal CDF.

```python
import numpy as np
from scipy import stats

def normalize_prs(scores, reference_scores=None):
    '''Z-score normalize PRS

    Args:
        scores: Array of PRS values
        reference_scores: Population reference (if None, use scores)

    Returns:
        Z-scored PRS values
    '''
    if reference_scores is None:
        reference_scores = scores
    mean = np.mean(reference_scores)
    std = np.std(reference_scores)
    return (scores - mean) / std

def prs_to_percentile(z_score):
    '''Convert Z-scored PRS to population percentile'''
    return stats.norm.cdf(z_score) * 100

# Example
prs_raw = np.array([0.5, 1.2, -0.3, 2.1, 0.8])
prs_z = normalize_prs(prs_raw)
percentiles = prs_to_percentile(prs_z)
```

## Risk Stratification

**Goal:** Categorize individuals into clinical risk groups based on their Z-scored PRS.

**Approach:** Apply population-distribution-based thresholds to assign Low/Average/High/Very High risk tiers.

```python
def stratify_risk(prs_z, thresholds=None):
    '''Categorize PRS into risk groups

    Default thresholds based on population distribution:
    - Low: < -1 SD (bottom 16%)
    - Average: -1 to 1 SD (middle 68%)
    - High: > 1 SD (top 16%)
    - Very high: > 2 SD (top 2.5%)
    '''
    if thresholds is None:
        thresholds = {'low': -1, 'high': 1, 'very_high': 2}

    if prs_z > thresholds['very_high']:
        return 'Very High Risk'
    elif prs_z > thresholds['high']:
        return 'High Risk'
    elif prs_z < threshol