---
title: "The simplePHENOTYPES v2 grammar"
subtitle: Samuel B Fernandes and Alexander E Lipka
author: "simplePHENOTYPES (`r Sys.Date()`)"
output:
  rmarkdown::html_document:
    toc: true
    toc_depth: 2
    number_sections: true
    theme: spacelab
    highlight: textmate

vignette: >
  %\VignetteIndexEntry{The simplePHENOTYPES v2 grammar}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#",
  fig.width = 7,
  fig.height = 4
)
set.seed(1)
```

```{r setup}
library(simplePHENOTYPES)
data("SNP55K_maize282_maf04")
geno <- SNP55K_maize282_maf04
```

This vignette walks through simulating phenotypes with the simplePHENOTYPES v2
grammar, using the bundled maize 282 panel throughout. The original
`create_phenotypes()` interface is covered at the end: it is still supported,
but it is superseded, and new work should use the grammar described here.

Two companion vignettes go deeper on the topics this one only introduces:

- `vignette("breeding-populations")` — crossing, selfing and doubled haploids:
  building a pedigree from real founders.
- `vignette("genetic-maps")` — what `cm` is for, and how to build or replace it.

# The example data

`SNP55K_maize282_maf04` is 10,650 markers scored on 280 inbred lines, coded
-1 / 0 / 1. The first five columns are marker metadata.

```{r data}
dim(geno)
geno[1:5, 1:8]
```

# The genetic map

`cm` holds a genetic map in centiMorgans. Earlier versions of the package left
this column empty, which meant nothing in the package could simulate meiosis.

```{r map}
aggregate(cm ~ chr, data = geno, FUN = max)
```

**This map is synthetic.** It is generated from the physical positions by
`synthetic_map()` and is not taken from any published maize linkage map, so its
cM values should not be quoted as measured recombination distances. It exists so
the crossing functions have realistic distances to work with. If you have a real
map for your population, use that instead — just put it in `cm`.

To build one for your own data:

```{r synthetic-map, eval = FALSE}
geno$cm <- synthetic_map(geno$chr, geno$pos)
```

The map only matters for the crossing functions; the grammar below does not use
it. See `vignette("genetic-maps")` for the recombination model, how to tune it,
and how to validate a real map.

# Reading genotype data

`as_numeric()` converts a genotype file into the numeric format used above,
handling HapMap, VCF, GDS and PLINK bed/ped input.

```{r as-numeric, eval = FALSE}
num <- as_numeric("my_genotypes.hmp.txt")
```

It handles the other formats the same way:

```{r as-numeric-vcf, eval = FALSE}
num <- as_numeric("my_genotypes.vcf", to_r = TRUE)
```

# Simulating phenotypes: the grammar

A simulation starts with `simulate_phenotype()`, which fixes the genetic
architecture, and then adds variance components as layers. Each layer's `prop`
is a requested marginal proportion of total phenotypic variance. Broad-sense
heritability is based on the additive, dominance, and epistatic genetic values;
the vQTL layer is residual heterogeneity and is not included.

```{r single-trait}
ph <- simulate_phenotype(geno, seed = 1) |>
  additive(prop = 0.5, n_qtn = 3)
ph
```

The pipe runs eagerly — there is no terminal `simulate()` call, and every object
already carries realized phenotypes.

```{r long-output}
head(phenotypes_long(ph))
```

## Dominance, epistasis and variance QTL

The other three layers work the same way. By default `dominance()` and `vqtl()`
reuse the additive QTNs; set `same_as_add = FALSE` to draw fresh ones.

```{r layers}
ph_ade <- simulate_phenotype(geno, seed = 2) |>
  additive(prop = 0.4, n_qtn = 5) |>
  dominance(prop = 0.1) |>
  epistasis(prop = 0.1, n_pairs = 2)
ph_ade
```

```{r vqtl}
ph_v <- simulate_phenotype(geno, seed = 3) |>
  additive(prop = 0.3, n_qtn = 4) |>
  vqtl(prop = 0.2)
ph_v
```

## Multiple traits and genetic correlation

`architecture = "pleiotropy"` puts the same QTNs behind several traits. With two
traits, `cor` sets the **genetic** correlation through the PleioArch algorithm,
which draws the pleiotropic QTN effects from a bivariate normal.

```{r pleiotropy}
ph_pleio <- simulate_phenotype(geno, architecture = "pleiotropy",
                               n_traits = 2, cor = 0.6, seed = 4) |>
  additive(prop = 0.5, n_qtn = 100)
ph_pleio

wide <- phenotypes_wide(ph_pleio)
round(cor(wide$Trait_1, wide$Trait_2), 2)
```

Two things are worth knowing before you compare that number with `cor`.

*It is measured on the wrong scale.* `cor` is the genetic correlation, but what
you observe above is the phenotypic one, and the residuals are drawn
independently. With uncorrelated residuals the two are related by
$r_P \approx r_G \sqrt{h^2_1 h^2_2}$, so a genetic correlation of 0.6 between
traits with $h^2 = 0.5$ shows up as a phenotypic correlation near 0.3.

*It is realized, not fixed.* `cor` is the target the effect draws aim at, not a
value hit exactly in any one simulation. The QTN effects are random draws, so the
realized correlation scatters around the target -- with few shared QTNs it is
also pulled slightly toward 0 on average -- and the scatter narrows as QTNs are
added (for loci in approximate linkage equilibrium), down to the sampling spread
of a correlation over the individuals you simulated. Measuring the genetic correlation directly, over 30 seeds:

```{r cor-scatter}
scatter <- data.frame(n_qtn = c(20, 100, 400), mean = NA_real_, sd = NA_real_)

for (i in seq_len(nrow(scatter))) {
  realized <- numeric(30)
  for (s in 1:30) {
    sim <- simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 2,
                              cor = 0.6, h2 = 0.5, seed = s) |>
      additive(n_qtn = scatter$n_qtn[i])
    g <- genetic_values(sim)
    realized[s] <- cor(g[, 1], g[, 2])
  }
  scatter$mean[i] <- round(mean(realized), 3)
  scatter$sd[i]   <- round(sd(realized), 3)
}

scatter
```

The mean stays near the target and the scatter narrows as QTNs are added, but on a
real panel not indefinitely, and not as $1/\sqrt{n_{qtn}}$: the correlation is
measured over a fixed number of individuals, whose sampling spread no number of
QTNs removes, linkage disequilibrium among the causal loci limits it further, and
with few shared QTNs
the realized value is also pulled slightly toward 0 on average. If you need a tight
correlation in one particular data set, use many QTNs — and check what you actually
got rather than assuming.

*Any number of traits.* The same engine covers three, five or twenty traits —
`cor` is either a scalar applied to every pair or a full `n_traits x n_traits`
matrix. The genetic covariance matrix must be positive semi-definite, which for
two traits is the familiar `cor^2 <= pi_1 * pi_2` and for more traits also rules
out mutually inconsistent requests; impossible targets are an error, not a
silent approximation.

```{r pleiotropy-multi}
R <- matrix(c(1, 0.8, -0.4,
              0.8, 1, -0.2,
             -0.4, -0.2, 1), 3, 3)

# Averaged over 20 seeds, since any single simulation scatters around the target
total <- matrix(0, 3, 3)
for (s in 1:20) {
  ph3 <- simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 3,
                            cor = R, h2 = 0.5, seed = s) |>
    additive(n_qtn = 300)
  total <- total + cor(genetic_values(ph3))
}

round(total / 20, 2)
```

*Major QTNs trade precision for realism.* `n_pleio_major` and `prop_var_major`
concentrate part of the pleiotropic variance in a few large-effect loci. Both
default to zero, because a locus carrying a large share of the variance makes
the realized correlation hinge on that one effect draw, and the convergence
above disappears — the standard deviation stops improving with `n_qtn`
altogether. Set them when a major pleiotropic gene is the object of study, not
otherwise.

## Linked, not pleiotropic: the "ld" architecture

`architecture = "ld"` is the spurious-pleiotropy scenario: two traits whose
**distinct** causal loci sit in linkage disequilibrium, so the traits covary
through linkage even though no locus is causal for both. It is a two-trait
construct (`n_traits = 2`) — one causal SNP per trait per locus-pair. `ld_type`
picks how the pair is linked: `"direct"` (the default) puts trait 1's and trait
2's causal SNPs directly in LD; `"indirect"` links each to a shared hidden,
non-causal "cause-of-LD" locus. `r2_min`/`r2_max` set the window the linked
pair's squared correlation must fall in. `qtn_table()` reports the linked pair
on each row: `QTN_t1` / `QTN_t2` (the two traits' causal SNPs) and `ld_r2` (the
r² between them).

```{r ld}
ph_ld <- simulate_phenotype(geno, architecture = "ld", n_traits = 2,
                            seed = 5) |>
  additive(prop = 0.4, n_qtn = 3)
ph_ld

# distinct causal loci per trait, with the r2 between the linked pair
qtn_table(ph_ld)[, c("trait", "snp", "QTN_t1", "QTN_t2", "ld_r2")]
```

## Combining architectures

`complex_phenotypes()` merges complete models under a common heritability,
weighting them by their genetic variances. This is how partial pleiotropy is
built: combine a pleiotropic model with an independent one.

```{r complex}
pleio <- simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 2,
                            cor = 0.5, seed = 6) |>
  additive(prop = 0.4, n_qtn = 5)
indep <- simulate_phenotype(geno, n_traits = 2, seed = 6) |>
  additive(prop = 0.3, n_qtn = 5)

both <- complex_phenotypes(pleio, indep, h2 = 0.5)
both
```

## One call instead of a pipe

Supplying `h2` together with `n_qtn` realizes a complete model in a single call.

```{r one-call}
simulate_phenotype(geno, h2 = 0.6, n_qtn = 5, model = "AD", seed = 7)
```

# Writing results

Phenotypes come back as an object; `write_phenotypes()` writes them to a file you
name, either as a long table (one row per individual and trait) or a wide one
(one column per trait). `sep` sets the field separator.

```{r exporters}
out_long <- file.path(tempdir(), "phenotypes_long.txt")
out_wide <- file.path(tempdir(), "phenotypes_wide.csv")
write_phenotypes(ph_pleio, file = out_long, format = "long")
write_phenotypes(ph_pleio, file = out_wide, format = "wide", sep = ",")
unlink(c(out_long, out_wide))
```

```{r wide}
head(phenotypes_wide(ph_pleio))
```

# Multi-generation crosses

The crossing functions simulate meiosis directly, so they can build a pedigree
from real founder genotypes. Phase the founders into a `Population`, then cross,
self or make doubled haploids:

```{r pedigree}
pop <- as_population(geno, individuals = c("33-16", "38-11"))
f1  <- cross(pop[1], pop[2], n = 1, seed = 11)
f2  <- selfcross(f1, n = 200, seed = 12)
f2
```

A `Population` is accepted anywhere `simulate_phenotype()` takes `geno`, so a
pedigree can be phenotyped without any conversion:

```{r phenotype-population}
ph_f2 <- simulate_phenotype(f2, seed = 14) |>
  additive(prop = 0.6, n_qtn = 5)
ph_f2
```

`vignette("breeding-populations")` covers this properly: backcrosses, recurrent
selfing to inbred lines, doubled haploids, and how to check that segregation and
recombination are behaving.

# Original implementation of create_phenotypes

`create_phenotypes()` is the original v1 function. It is **frozen**: it keeps
its original code paths and seed handling so that existing scripts keep
reproducing their published results, and it receives bug fixes only. New work
should use the grammar above, which is why it is documented as superseded.

```{r legacy, eval = FALSE}
create_phenotypes(
  geno_obj = geno,
  add_QTN_num = 3,
  add_effect = 0.2,
  rep = 10,
  h2 = 0.7,
  model = "A",
  to_r = TRUE,
  home_dir = tempdir()
)
```

# Session information

```{r session}
sessionInfo()
```