---
title: "R/qtl2pattern Vignette"
author: "Brian S. Yandell"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{demo qtl2pattern features}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 7, fig.height = 5)
```

**Look for bold stuff to remove or adjust code.**
**Note that changes to summary broke part of qtl2shiny in SNP/Gene Action**

This vignette continues an example developed in library [qtl2ggplot](https://CRAN.R-project.org/package=qtl2ggplot) using data
from [qtl2data](https://github.com/rqtl/qtl2data) with tools found in [qtl2pattern](https://github.com/byandell/qtl2pattern).
Calculations are repeated from that example, but only plots that help tell the story.
Focus is entirely on chromosome `2` of the [DOex](https://github.com/rqtl/qtl2data/tree/master/DOex) from 'qtl2data'. See vignette in 'qtl2ggplot' for complementary analysis.

Package 'qtl2pattern' does not depend on 'qtl2ggplot', and can be used independently.
However this vignette shows some new 'ggplot2' based routines that extend 'qtl2ggplot'
to the functionality of this package.
This vignette illustrates these 'ggplot2'-derived routines, which is why 'qtl2ggplot' is
suggested.

Package `qtl2pattern` has `summary` methods for `scan1` and `scan1coef`.
Not sure these are really needed. The `R/qtl2` routine `find_peaks` takes care of some but not all of this functionality.
Sections below are

- Download of Diversity Outbred Example Data
- Genome Scan Using 36 Allele Pairs
- Genome Scan Using 8 Alleles
- SNP Association Mapping near Peak
- Strain Distribution Pattern (SDP) Scan

```{r}
library(qtl2pattern)
```

## Download of Diversity Outbred Example Data

This vignette uses the example Diversity Outbred (DO) data.
See [Recla et al. (2014)](https://pubmed.ncbi.nlm.nih.gov/pubmed/24700285/).
These mice were derived from 8 founder strains, yielding up to 36 allele pairs at any marker.
The data can be found in <https://github.com/rqtl/qtl2data> as `DOex`.
While these data span three chromosomes (`2, 3, X`), we focus here on
chromosome `2`.

```{r}
dirpath <- "https://raw.githubusercontent.com/rqtl/qtl2data/master/DOex"
```

```{r}
DOex <-
  subset(
    qtl2::read_cross2(file.path(dirpath, "DOex.zip")),
    chr = "2")
```

Download 36 genotype probabilities for allele pairs by individual and marker.

```{r}
tmpfile <- tempfile()
download.file(file.path(dirpath, "DOex_genoprobs_2.rds"), tmpfile, quiet=TRUE)
pr <- readRDS(tmpfile)
unlink(tmpfile)
```

Or alternatively calculate those genotype probabilities

```{r eval = FALSE}
pr <- qtl2::calc_genoprob(DOex, error_prob=0.002)
```

## SNP Association Mapping near Peak

This section briefly examines SNP association mapping,
which involves collapsing the 36 allele pair genotype probabilities to
3 allele probabilities in the region of a QTL peak.

Download snp info from web and read as RDS.

```{r}
tmpfile <- tempfile()
download.file(file.path(dirpath, "c2_snpinfo.rds"), tmpfile, quiet=TRUE)
snpinfo <- readRDS(tmpfile)
unlink(tmpfile)
```

Rename the SNP position as `pos`,
and add SNP index information.

```{r}
snpinfo <- dplyr::rename(snpinfo, pos = pos_Mbp)
snpinfo <- qtl2::index_snps(DOex$pmap, snpinfo)
```

Convert genotype probabilities to SNP probabilities. Notice that starting with the allele probabilities only ends up with 2 SNP alleles per marker (`"A", "B"`), but we want all three (`"AA" "AB" "BB"`). Recall that the SNP alleles are in context of the strain distribution pattern (`sdp`) for that SNP.
The `sdp` for markers (`rsnnn`) are imputed in 'qtl2' from the adjacent SNPs.
That is, using allele probabilities,

```{r}
apr <- qtl2::genoprob_to_alleleprob(pr)
snpapr <- qtl2::genoprob_to_snpprob(apr, snpinfo)
dim(snpapr[["2"]])
dimnames(snpapr[["2"]])[[2]]
```

while using genotype probabilities,

```{r}
snppr <- qtl2::genoprob_to_snpprob(pr, snpinfo)
dim(snppr[["2"]])
dimnames(snppr[["2"]])[[2]]
```

```{r}
rm(snpapr)
```

Perform SNP association analysis (here, ignoring residual kinship).

```{r}
scan_snppr <- qtl2::scan1(snppr, DOex$pheno)
```

Package 'qtl2' provides a summary of the peak using `qtl2::find_peaks`.

```{r}
qtl2::find_peaks(scan_snppr, snpinfo)
```

## Strain Distribution Pattern (SDP) Scan

Strain distribution patterns (SDPs) separate out SNPs based on their SDP and plot the top patterns.
For instance `sdp = 52` corresponds to pattern `ABDGH:CEF`. That is, the SNP genotype `"AA"` resulting from `qtl2::genoprob_to_snpprob` applied to `pr` corresponds to any of the 36 allele pairs with the two alleles drawn from the reference (`ref`) set of `ABDGH` (15 pairs: `AA, AB, AD, AG, AH, BB, BD, BG, BH, DD, DG, DH, GG, GH, HH`), `"BB"` has two alleles from the alternate (`alt`) set `CEF` (6 pairs: `CC, CE, CF, EE, EF, FF`), and `"AB"` has one from each for the heterogeneous (`het`) set (15 pairs: `AC, AE, ..., HF`).
There are 255 possible SDPs, but only a few (4 in our example) that need be examined carefully. One can think of these as a subset of markers for
genome scan, where interest is only in those SNPS following a particular `sdp`; as with genome scans, we can fill in for missing data.
That is, only a few SNPs may show a particular pattern, but key differences might be seen nearby if we impute SNPs of the same pattern.

The `top_snps_pattern` routine is an extension of `qtl2::top_snps`, which provides more detail on SDPs.

```{r}
top_snps_tbl <- top_snps_pattern(scan_snppr, snpinfo)
```

This default `summary` is nearly identical to the `summary` of the SNP scan object above:

```{r}
(patterns <- summary(top_snps_tbl))
```

There may be multiple SNPs identified for an SDP with the same LOD, covering a range of positions.
If there is a range of `lod` values, there are additional SNPs beyond those reported in the summary with `lod` values below the maximum (see `ABDFGH:CE` and `ABCDFGH:E` patterns in plot above).
The following summary shows details for the first 10 SNPS with top `lod` per SDP:

```{r}
head(summary(top_snps_tbl, "best"))
```

## SDP Pattern Scan

A new routine `scan1pattern` scans the peak region for each of the 4 patterns provided.
That is, the SDP scan only considers SNPs that have a particular strain distribution pattern and uses the SNP probabilities. However, for each markers (or any genome position), one can impute a SNP with any particular SDP and compute its SDP probabilities. That is, we could consider the SDP probabilities at every marker position for the SDP 48 (pattern `ABCDGH:EF`) and do an SDP scan with these SDP probabilities. That is what the routine `scan1pattern` does.

```{r}
scan_pat <- scan1pattern(pr, DOex$pheno,
                         map = DOex$pmap,
                         patterns = patterns)
```

```{r}
summary(scan_pat, DOex$pmap)
```

Here is a scan around the region where we have SNP information. The SDPs 52 and 53 are both higher than the other two patterns, and sustain the SDP peak several cM to the right of the allele peak region. This difference in pattern, if real, might suggest some form of allele interaction is important.

```{r}
pat_names <- paste(c(52,43,16,20), colnames(scan_pat$scan), sep = "_")
```

```{r}
plot(scan_pat$scan, DOex$pmap, lodcolumn = 2, col = "red", xlim = c(90,110), type = "b")
abline(v = c(96.5, 98.5), col = "darkgray", lwd = 2, lty = 2)
plot(scan_pat$scan, DOex$pmap, lodcolumn = 1, add = TRUE, col = "blue", type = "b")
plot(scan_pat$scan, DOex$pmap, lodcolumn = 3, add = TRUE, col = "purple", type = "b")
plot(scan_pat$scan, DOex$pmap, lodcolumn = 4, add = TRUE, col = "green", type = "b")
title("Scans for SDP 52 (blue), 43 (red), 16 (purple), 20 (green)")
```

Here is a scan over the whole chromosome, imputing SDP. These correspond to reducing the 8 alleles to 2 corresponding to the `ref` (`ABDGH` for sdp 52) and `alt` (`CEF` for sdp 52) composite alleles for the blue curve.

```{r}
plot(scan_pat$scan, DOex$pmap, lodcolumn = 2, col = "red")
abline(v = c(96.5, 98.5), col = "darkgray", lwd = 2, lty = 2)
plot(scan_pat$scan, DOex$pmap, lodcolumn = 1, add = TRUE, col = "blue")
plot(scan_pat$scan, DOex$pmap, lodcolumn = 3, add = TRUE, col = "purple")
plot(scan_pat$scan, DOex$pmap, lodcolumn = 4, add = TRUE, col = "green")
title("Scans for SDP 52 (blue), 43 (red), 16 (purple), 20 (green)")
```

### SDP Effects Scan

In a similar fashion to a two-allele cross, we can examine the SDP effects. The top row of the plot below is for the patterns with the higher LOD score.

```{r}
oldpar <- par(mfrow = c(2,2))
cols <- c("blue","purple","red")
for(i in 1:4) {
  plot(scan_pat$coef[[i]], DOex$pmap, columns = 1:3, col = cols, xlim = c(90,110), type = "b")
  title(paste("Coefficients for", pat_names[i]))
  if(i == 3) {
    abline(v = c(96.5, 98.5), col = "darkgray", lwd = 2, lty = 2)
    legend(104, -3, legend = c("ref","het","alt"), 
       col = cols, lty = 1, lwd = 2)
  }
}
par(oldpar)
```

From these, the coefficients for pattern `ABDGH:CEF` (sdp 52) appear to be problematic, whereas pattern `ABCDGH:EF` show a clear pattern of dominance of the reference allele, consistent with the earlier plot of allele coefficients where `E = NZO` and `F = CAST`. The pattern `ABCFGH:CE` shows a similar dominance, but the allele `C = 129` does not stand out in the allele plot. This could be that `129` is important in combination with `NZO`. This is even more apparent when looking at coefficients over the entire chromosome.

```{r}
par(mfrow = c(2,2))
cols <- c("blue","purple","red")
for(i in 1:4) {
  plot(scan_pat$coef[[i]], DOex$pmap, columns = 1:3, col = cols)
  title(paste("Coefficients for", pat_names[i]))
  if(i == 3) {
    abline(v = c(96.5, 98.5), col = "darkgray", lwd = 2, lty = 2)
    legend(125, -3, legend = c("ref","het","alt"), 
       col = cols, lty = 1, lwd = 2)
  }
}
```

### 36 Allele Pair Scan

We are now going to look at some of the 36 coefficients. These are unwieldy, so we collapse the allele pairs for SDP 52, pattern `ABCDGH:EF` into `"ref", "het", "alt"` for summaries, looking individually only at pairs `"EE", "EF", "FF"`.

```{r}
coefs2 <- qtl2::scan1coef(pr, DOex$pheno[,"OF_immobile_pct"], zerosum = FALSE)
```

The following is messy code to just pull out weighted averages corresponding to reference, het, and alternative genotypes at the `pr` peak.

```{r}
x <- LETTERS[1:8]
ref <- outer(x,x,paste0)
ref <- ref[upper.tri(ref, diag=TRUE)]
ss <- 2 - sapply(stringr::str_split(ref, ""), function(x) x[1] == x[2])
alt <- ref[!stringr::str_detect(ref, c("A|B|C|D|G|H"))]
het <- ref[stringr::str_detect(ref, c("A|B|C|D|G|H")) & stringr::str_detect(ref, c("E|F"))]
altw <- ss[match(alt, ref, nomatch = 0)]
hetw <- ss[match(het, ref, nomatch = 0)]
refw <- ss[!stringr::str_detect(ref, c("E|F"))]
ref <- ref[!stringr::str_detect(ref, c("E|F"))]
```

```{r}
coefsum <- coefs2
coefsum[,"AA"] <- apply(coefsum[,ref], 1, weighted.mean, w = refw)
coefsum[,"AB"] <- apply(coefsum[,het], 1, weighted.mean, w = hetw)
coefsum[,"BB"] <- apply(coefsum[,alt], 1, weighted.mean, w = altw)
colnames(coefsum)[1:3] <- c("ref","het","alt")
```

```{r}
plot(coefsum, DOex$pmap, c("ref", "het", "alt", alt), xlim = c(90,110), ylim = c(-500, 500),
     col = c(1,8,7,2:4), type = "b")
abline(h = mean(DOex$pheno[,"OF_immobile_pct"]), lwd = 2, col = "darkgrey", lty = 2)
abline(v = c(96.5, 98.5), col = "darkgray", lwd = 2, lty = 2)
title("Allele Coefficients for OF_immobile_pct")
legend(105, 80, legend = c("ref", "het", "alt", alt), 
       col = c(1,8,7,2:4), lty = 1, lwd = 2)
```

Estimates for allele pairs are imprecise, so let's drop the extreme `"EE"`. The `"ref"`, `"het"` and `"alt"` are weighted mean of 15, 6, 15 allele pairs, respectively, assuming Hardy-Weinberg equilibrium.
We see that the allele pairs `"EF", "FF"` lie near `"het"` and below `"ref"` in the peak region
in the graph below.

```{r}
plot(coefsum, DOex$pmap, c("ref", "het", alt[-1]), xlim = c(90,110), ylim = c(55,90),
     col = c(1,8,3:4), type = "b")
abline(h = mean(DOex$pheno[,"OF_immobile_pct"]), lwd = 2, col = "darkgrey", lty = 2)
abline(v = c(96.5, 98.5), col = "darkgray", lwd = 2, lty = 2)
title("Allele Pair Coefficients for OF_immobile_pct")
legend(107, 75, legend = c("ref", "het", alt[-1]), 
       col = c(1,8,3:4), lty = 1, lwd = 2)
```

### Effect estimates

**These look off somehow. Probably don't have best choice of position.**
**Also, allele1 is fragile if you don't provide stuff.**
**Also, want to add plot of data.**

```{r}
scan_apr <- qtl2::scan1(apr, DOex$pheno)
coefs <- qtl2::scan1coef(apr, DOex$pheno)
coefs2 <- qtl2::scan1coef(pr, DOex$pheno)
```

```{r}
alleles <- allele1(probD = pr,
                   scanH = scan_apr, coefH = coefs, coefD = coefs2, 
                   scan_pat = scan_pat, map = DOex$pmap, alt = "E")
```

```{r}
ggplot2::autoplot(alleles, scan_apr, DOex$pmap, frame = FALSE)
```

Look at subset with haplotype and the 3 well-behaved SDPs.

```{r}
aa <- subset(alleles, sources = levels(alleles$source)[c(1,4:6)])
ggplot2::autoplot(aa, scan_apr, DOex$pmap, frame = FALSE)
```

The following plot of data uses the closest marker to the SDP 48 max. Hopefully this marker has SDP 48. We could likely improve things with more work by using the SDP probabilities. Check sign on calculations.

```{r}
pos = patterns$max_pos[patterns$sdp == 48]
wh <- which.min(abs(pos - DOex$pmap[[1]]))
geno <- apply(snppr[[1]][,,wh], 1, function(x) which.max(x))
geno <- factor(c("ref","het","alt")[geno], c("ref","het","alt"))
dat <- data.frame(DOex$pheno, geno = geno)
dat$x <- jitter(rep(1, nrow(dat)))
```

```{r}
ggplot2::ggplot(dat) +
  ggplot2::aes(x = x, y = OF_immobile_pct, col = geno, label = geno) +
  ggplot2::geom_boxplot() +
  ggplot2::geom_point() +
  ggplot2::facet_wrap(~ geno)
```

These plots anticipate some of the possibilities with mediation.

## SNP features and Gene action

Download Gene info for DOex from web via RDS. For the example below,
genes are presented without exons. 

```{r}
tmpfile <- tempfile()
download.file(file.path(dirpath, "c2_genes.rds"), tmpfile, quiet=TRUE)
gene_tbl <- readRDS(tmpfile)
unlink(tmpfile)
class(gene_tbl) <- c("feature_tbl", class(gene_tbl))
```

Variants are usually SNPs, but may be of other types, such as deletions (`DEL`), insertions (`INS`), or insertions and deletions (`InDel`). The routine cannot yet handle major chromosomal rearrangements and translocations. In addition variants may be intronic, exonic, or intergenic, depending on the "consequence" as reported in the variant database.

```{r}
out <- merge_feature(top_snps_tbl, snpinfo, scan_snppr, exons = gene_tbl)
summary(out, "pattern")
```

Add bogus consequence for show.

```{r}
out$snp_type <- sample(c(rep("intron", 40), rep("exon", nrow(out) - 40)))
```

```{r}
ggplot2::autoplot(out, "OF_immobile_pct")
```

```{r}
ggplot2::autoplot(out, "OF_immobile_pct", "consequence")
```

It is also possible to consider different types of gene action for SDPs. With no restrictions, there are 2 degrees of freedom (`add+dom`). Various one df options include `additive`, `dominant`, `recessive` and `non-add`itive. It is also possible if sex is encoded to separate analysis by sex.

**This currently ignores snp_type**

## Gene and exon scan

**Add fake exons.**

```{r eval = FALSE}
# Not quite right.
gene <- "Lrrc4c"
geneinfo <- dplyr::filter(gene_tbl, Name == gene)
parts <- geneinfo$start + diff(unlist(geneinfo[,c("start","stop")])) * (0:5) / 5
geneinfo <- geneinfo[rep(1,3),]
geneinfo$type <- "exon"
geneinfo$start <- parts[c(1,3,5)]
geneinfo$stop <- parts[c(2,4,6)]
gene_tbl <- rbind(gene_tbl[1,], geneinfo, gene_tbl[-1,])
```


Plot Genes within some distance of high SNPs.

```{r}
qtl2::plot_genes(gene_tbl, xlim = c(96,99))
```

```{r}
ggplot2::autoplot(gene_tbl)
```

Adding the result of `top_snps_pattern` overlays significant SNPS on the plot of genes.

```{r}
ggplot2::autoplot(gene_tbl, top_snps_tbl = top_snps_tbl)
```
The routine `gene_exon` examines individual genes and their exons.
This example does not have exons, so there is just the figure of the gene. Exons would appear on different rows. Again, the vertical dashed lines correspond to SNPs with high LODs.

```{r}
# Get Gene exon information.
out <- gene_exon(top_snps_tbl, gene_tbl)
summary(out, gene = out$gene[1])
```

```{r}
ggplot2::autoplot(out, top_snps_tbl)
```

## Database Query and Large Crosses

The package [R/qtl2](https://kbroman.org/qtl2/) is designed for large crosses, in which phenotype, genotype and variant information may require substantial space and efficient access and calculations. Here is some additional information that is useful when working with large crosses. See the following information at :

- [R/qtl2 User Guide](https://kbroman.org/qtl2/assets/vignettes/user_guide.html)
  + [R/qtl2 Connecting to SNP and gene databases](https://kbroman.org/qtl2/assets/vignettes/user_guide.html#Connecting_to_SNP_and_gene_databases)
- [qtl2fst user guide](https://kbroman.org/qtl2/assets/vignettes/qtl2fst.html)

Large systems genetics studies have databases that need not be uploaded to R, but are better accessed
via query searchers. Package 'qtl2' has two functions to create query functions for SNP and other variants (`create_variante_query_func`) and genes (`create_gene_query_func`). Package 'qtl2pattern' adds a function for genotype probabilities (`create_probs_query_func`).

Here are query function creations using the 'qtl2' routines, assuming that variants and genes are stored in an `SQLite` file database in folder `qtl2shinyData/CCmouse`:

```{r}
query_variants <- 
  qtl2::create_variant_query_func(
    file.path("qtl2shinyData", "CCmouse", "cc_variants.sqlite"))
```

```{r}
query_genes <- 
  qtl2::create_gene_query_func(
    file.path("qtl2shinyData", "CCmouse", "mouse_genes.sqlite"))
```

Here are the objects in the local `query_variants()` environment:

```{r}
objects(envir = environment(query_variants))
```

Here is the value of `dbfile` in the local `query_variants()` environment:

```{r}
get("dbfile", envir = environment(query_variants))
```

While these can be used interactively, their real value comes in saving them for
use separately, such as with 'qtl2shiny'. This can be done using `saveRDS`, with later access using `readRDS`.

```{r eval = FALSE}
saveRDS(query_variants, "query_variants.rds")
#
# other code
#
query_variants <- readRDS("query_variants.rds")
```

The new routine in 'qtl2pattern' is for genotype probabilities, which can now be quite large when considering SNP probabilities. These are also specific to a cross, whereas the variants and genes are specific to a taxa. The strategy is to place genotype probabilities in a folder (default names is `genoprob`) either in RDS or FST format.

```{r eval = FALSE}
query_probs <- 
  create_probs_query_func(
    file.path("qtl2shinyData", "CCmouse", "Recla"))
```

The current preferred storage for large probability files uses package [fst](https://www.fstpackage.org/) via the package [qtl2fst](https://kbroman.org/qtl2/assets/vignettes/qtl2fst.html).
Here is some of the size information from the full Recla genotype data.

```
> dirpath <- "inst/qtl2shinyApp/qtl2shinyData/CCmouse/Recla/genoprob/"
> fpr <- readRDS(file.path(dirpath, "fst_probs.rds"))
> print(object.size(fpr), units = "Mb")
1.5 Mb
> dim(fpr)
      1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19   X
ind 261 261 261 261 261 261 261 261 261 261 261 261 261 261 261 261 261 261 261 261
gen  36  36  36  36  36  36  36  36  36  36  36  36  36  36  36  36  36  36  36  44
mar 536 542 469 452 438 444 423 373 386 430 393 355 362 369 303 299 270 260 196 439
>
> file.info(file.path(dirpath, "fst_probs.rds"))$size
[1] 70827
> file.info(file.path(dirpath, "probs_1.fst"))$size
[1] 40314202
```

Thus the master FST file `"fst_probs.rds"` external to 'R' is 71Kb, while the FST database for Chr 1 external to 'R' is 40Mb. There is a separate file for each chromosome. All combined, the external storage for the master FST and chromosome-specific data is over 800Mb; adding the allele probabilities, which are much smaller, brings the total to ~1Gb.
The storage for Recla genotype probabilities is roughly proportional to the number of individuals (261), number of markers (7739), and number of allele pairs (36), that is about 14b per individual, marker and allele pair.  

Internal to 'R', the FST genotype probability object `fpr` is 1.5Mb, and contains all the information about where to find the genotype probability information for individual chromosomes. One uses this object as one uses the genotype probability object `pr`, created by `qtl2::calc_genoprob`, using functions in 'qtl2'. However, it is much smaller. Typically, 'qtl2' functions will only need a small part of the genotype probabilities at any time. See documentation for [qtl2fst](https://kbroman.org/qtl2/assets/vignettes/qtl2fst.html) for more information.
