physio_trait.Rmd

---
author: "Brian Yandell"
date: "`r format(Sys.time(), '%d %B %Y')`"
output: html_document
params:
  physioshort: "cpep_ins"
  physioname: "log10_cpep_ins_molar_ratio_15_14wk"
  genename: "Wfdc21"
  isoformid: ""
  chrname: "11"
  chrwidth: 1
  lod_drop: 3
  snp36: TRUE
  normscore: TRUE
  datadir: "data"
---
---
title: "`r params$physioshort` with `r params$genetrait` on chr `r params$chr`"
---

```{r setup, include=FALSE, message=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```
From: Tara Price <trprice@wisc.edu>
Date: Wednesday, December 13, 2023 at 4:02 PM
To: Brian Yandell <brian.yandell@wisc.edu>
Subject: Puzzling Cholesterol QTL
Hi Brian, 
Attaching a powerpoint with QTL scans and BLUPs for the cholesterol QTL on Chromosome 5. The first 5 slides are the QTL and their allele effects and how we found a candidate gene based on matching allele effects (old DO) and biological plausibility. Then, we ran the eQTL for Scarb1 using liver RNA from the new, diet DO study and the allele effects did not match the phenotype or old DO. 

The latter half of the slides are just pictures I generated from running the QTL scans/BLUPs when the two diets were separated. Essentially, there were slight differences in the SNP plots between the two diets, but the BLUPs looked very similar. I also added Scarb1 as a covariate to my scans, but that had very little effect on the cholesterol QTL. 

This slide deck focuses on the Chromosome 5 QTL, but I am also having difficulty identifying a plausible candidate on Chromosome 1 (despite Apoa2, a known cholesterol protein, being at that same locus). 

Any guidance/insight/suggestions you can provide will be GREATLY appreciated.
Please let me know if you need other information from me. I’m also happy to arrange a time and talk through my triage steps, if that’s useful.
```

```
From: MARK P KELLER <mark.keller@wisc.edu>
Date: Friday, December 8, 2023 at 3:52 PM
To: Brian Yandell <brian.yandell@wisc.edu>
Subject: percent_fat_17wk

Based on literature, the most plausible gene under the QTL is: Glb1.
 
There are only two isoforms for Glb1:
ENSMUST00000063042
ENSMUST00000217583
 
The one highlighted yellow has a cis-eQTL with a similar allele pattern as the fast mass QTL; the other does not.
```

# Data Setup

```{r}
library(dplyr)
library(tidyr)
library(ggplot2)
library(qtl2)
library(qtl2fst)
library(qtl2ggplot)
```

```{r}
chrname <- params$chrname
physioshort <- params$physioshort
physioname <- params$physioname
genename <- params$genename
isoformid <- params$isoformid
lod_drop <- params$lod_drop
chrwidth <- params$chrwidth
snp36 <- as.logical(params$snp36) # use 36 genotypes for SNPs
normscore <- as.logical(params$normscore) # use normal scores (nqrank)
datadir <- as.character(params$datadir)
```

**Clinical trait `r physioname` (`r physioshort`) with covariate `r genename` (optional isoform `r isoformid`) on chr `r chrname`.
Window half-width of `r chrwidth` with SNP LOD drop of `r lod_drop` based on
`r ifelse(snp36, "36 allele pairs", "8 alleles")`.
Data are `r ifelse(normscore, "", "not")` transformed by normal scores.
Data is folder `datadir`.**
Source code at
[physio_trait.Rmd](https://github.com/byandell/DO_Diet/blob/master/physio_trait.Rmd).

# Data Setup

```{r}
source("R/setup.R")
```

# Trait Scan

```{r}
trait_addscan_chr <- qtl2::scan1(aprobs_chr, physio_trait, kinship_chr, addcovar,
                             cores = 0)

trait_intscan_chr <- qtl2::scan1(aprobs_chr, physio_trait, kinship_chr, addcovar,
                             intcovar = intcovar, cores = 0)

plot(trait_intscan_chr, pmap, col = "red")
plot(trait_addscan_chr, pmap, add = TRUE)
```

```{r}
(trait_peaks <- find_peaks(trait_intscan_chr, pmap, threshold = 9, drop = 2))
```

## Zoom in to small window

```{r}
plot(trait_intscan_chr, pmap, col = "red",
     xlim = trait_peaks$pos + c(-1,1) * chrwidth)
plot(trait_addscan_chr, pmap, add = TRUE,
     xlim = trait_peaks$pos + c(-1,1) * chrwidth)
abline(v = trait_peaks$pos, col = "gray", lwd = 2)
```

# Strain Coefficients

```{r}
trait_nocoef <- scan1coef(aprobs_chr, physio_trait, kinship_chr)
plot_coef(trait_nocoef, pmap[chrname], columns = 1:8)
```

```{r}
trait_addcoef <- scan1coef(aprobs_chr, physio_trait, kinship_chr, addcovar)
plot_coef(trait_addcoef, pmap[chrname], columns = 1:8,
          xlim = trait_peaks$pos + c(-1,1))
abline(v = trait_peaks$pos, col = "grey", lwd = 2)
title(paste(physioshort, "on chr", chrname))
```

The following plot is flawed as it only presents part of the solution.

```{r}
trait_intcoef <- scan1coef(aprobs_chr, physio_trait, kinship_chr, addcovar,
                           intcovar = intcovar)

plot_coef(trait_intcoef, pmap[chrname], columns = 1:8,
          xlim = trait_peaks$pos + c(-1,1))
abline(v = trait_peaks$pos, col = "grey", lwd = 2)
title(paste(physioshort, "on chr", chrname, "additive part of diet interaction"))
```

Following does separate coefficients by diet. 
It is complicated by fact that `kinship` calculation does not include one mouse.

```{r}
wh_HF <- (addcovar[,"Diet.nameHF"] == 1)
m <- match(rownames(addcovar), rownames(kinship_chr[[1]]))
#rownames(addcovar)[is.na(m)]
id_HF <- rownames(addcovar)[wh_HF & !is.na(m)]
id_HC <- rownames(addcovar)[!wh_HF & !is.na(m)]

trait_coef_HF <- scan1coef(
  subset(aprobs_chr, ind = id_HF),
  physio_trait[id_HF,, drop = FALSE],
  kinship_chr, addcovar, intcovar = intcovar)

plot_coef(trait_coef_HF, pmap[chrname], columns = 1:8,
          xlim = trait_peaks$pos + c(-1,1))
abline(v = trait_peaks$pos, col = "grey", lwd = 2)
title(paste("HF Diet", physioshort, "on chr", chrname))

trait_coef_HC <- scan1coef(
  subset(aprobs_chr, ind = id_HC),
  physio_trait[id_HC,, drop = FALSE],
  kinship_chr, addcovar, intcovar = intcovar)

plot_coef(trait_coef_HC, pmap[chrname], columns = 1:8,
          xlim = trait_peaks$pos + c(-1,1))
abline(v = trait_peaks$pos, col = "grey", lwd = 2)
title(paste("HC Diet", physioshort, "on chr", chrname))
```

# SNP Associations

SNP genotype probabilities are calculated as 2-level additive (using `aprobs_chr` from the 8-allele genotype probablities) or 3-level full (using `probs_chr` from the 36-allele-pair genotype probabilities).
The latter are used with `snp36` is `TRUE` (default) and the 36-pair
probabilities are available.

```{r}
library(ggplot2)
library(qtl2ggplot)
```

```{r}
chrpos <- trait_peaks$pos
```

```{r}
snpdb_file <- "data/cc_variants.sqlite"
query_variant <- qtl2::create_variant_query_func(snpdb_file)
snpinfo <- query_variant(chrname,
                         chrpos - chrwidth,
                         chrpos + chrwidth)
```

Need to change names of alleles

```{r}
names(snpinfo)[8:15] <- names(qtl2::CCcolors)
names(snpinfo)[1] <- "snp_id"
```

Perform SNP association mapping.
It is possible to use `qtl2::scan1snps` instead, which bundles these three routines,
but we want to have the SNP probabilities for later use.

```{r}
snpinfo <- qtl2::index_snps(pmap, snpinfo)
```

```{r}
if(snp36 & file.exists(file.path(datadir, "fst_genoprob36"))) {
  snppr <- qtl2::genoprob_to_snpprob(
    subset(
      readRDS(file.path(datadir, "fst_genoprob36", "probs_fstindex.rds")),
      chr = chrname),
    snpinfo)
} else {
  snppr <- qtl2::genoprob_to_snpprob(aprobs_chr, snpinfo)
}
```

```{r}
scan_snppr <- qtl2::scan1(snppr, physio_trait, kinship_chr, addcovar,
                          intcovar = intcovar)
```

```{r}
autoplot(scan_snppr, snpinfo, drop_hilit = lod_drop)
```

```{r}
autoplot(scan_snppr, snpinfo, patterns="hilit", drop_hilit=lod_drop, cex=2)
```

```{r}
autoplot(scan_snppr, snpinfo, patterns="all", drop_hilit=lod_drop, cex=2)
```


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

Find top snps for distinct SDPs within `r `lod_drop` of peak.

```{r}
(top_snps_tbl <-
   qtl2pattern::top_snps_pattern(scan_snppr, snpinfo,
                                 drop = lod_drop,
                                 show_all_snps = FALSE) |>
  select(snp_id, pos, sdp, index, lod) |>
  distinct(sdp, .keep_all = TRUE) |>
  arrange(desc(lod)))
```

# Use best SNPs as covariate

```{r}
snp_id <- top_snps_tbl$snp_id
snp_df <- dim(snppr[[1]])[2] - 1
snp_geno <- matrix(NA, nrow(addcovar), length(snp_id) * snp_df)
rownames(snp_geno) <- rownames(addcovar)
if(snp_df == 1) {
  colnames(snp_geno) <- snp_id
} else {
  colnames(snp_geno) <- paste(rep(snp_id, each = 2), seq(snp_df), sep = "_")
}
snp_geno[dimnames(snppr[[1]])[[1]],] <- snppr[[1]][, -(1 + snp_df), snp_id]
```

```{r}
addcovar_geno <- cbind(addcovar, snp_geno)
addcovar_geno1 <- cbind(addcovar, snp_geno[,seq(snp_df), drop = FALSE])
```

```{r}
trait_genoscan_chr <- qtl2::scan1(aprobs_chr, physio_trait,
                                  kinship_chr, addcovar_geno,
                                  intcovar = intcovar, cores = 10)
trait_geno1scan_chr <- qtl2::scan1(aprobs_chr, physio_trait,
                                  kinship_chr, addcovar_geno1,
                                  intcovar = intcovar, cores = 10)
```

```{r}
bind_rows(
  add = qtl2::find_peaks(trait_addscan_chr, pmap),
  int = qtl2::find_peaks(trait_intscan_chr, pmap),
  geno1 = qtl2::find_peaks(trait_geno1scan_chr, pmap),
  geno = qtl2::find_peaks(trait_genoscan_chr, pmap),
  .id = "scan")
```

```{r}
plot(trait_intscan_chr, pmap, col = "red",
     xlim = trait_peaks$pos + c(-1,1))
plot(trait_genoscan_chr, pmap, col = "blue", add = TRUE,
     xlim = trait_peaks$pos + c(-1,1))
plot(trait_geno1scan_chr, pmap, col = "purple", add = TRUE,
     xlim = trait_peaks$pos + c(-1,1))
plot(trait_addscan_chr, pmap, col = "green", add = TRUE,
     xlim = trait_peaks$pos + c(-1,1))
abline(v = trait_peaks$pos, col = "gray", lwd = 2)
abline(v = snp_peak$pos, col = "blue", lwd = 2)
```

# Add gene and isoform as covariate

```{r}
dat <- as.data.frame(cbind(physio_trait, gene_trait, isoform_trait))
names(dat) <- c(physioshort, genename, isoformid)
ggplot(dat) +
  aes(.data[[physioshort]], .data[[genename]]) +
  geom_point() +
  geom_smooth(se = FALSE)
```

```{r}
ggplot(dat) +
  aes(.data[[physioshort]], .data[[isoformid]]) +
  geom_point() +
  geom_smooth(se = FALSE)
```

```{r}
ggplot(dat) +
  aes(.data[[genename]], .data[[isoformid]]) +
  geom_point() +
  geom_smooth(se = FALSE)
```

```{r}
addcovar_gene <- cbind(addcovar, gene_trait)
addcovar_isoform <- cbind(addcovar, isoform_trait)
addcovar_gene_isoform <- cbind(addcovar, gene_trait, isoform_trait)
```

```{r}
trait_genescan_chr <- qtl2::scan1(aprobs_chr, physio_trait,
                                  kinship_chr, addcovar_gene,
                                  intcovar = intcovar, cores = 10)
trait_isoformscan_chr <- qtl2::scan1(aprobs_chr, physio_trait,
                                  kinship_chr, addcovar_isoform,
                                  intcovar = intcovar, cores = 10)
trait_gene_isoscan_chr <- qtl2::scan1(aprobs_chr, physio_trait,
                                  kinship_chr, addcovar_gene_isoform,
                                  intcovar = intcovar, cores = 10)
```

```{r}
bind_rows(
  add = qtl2::find_peaks(trait_addscan_chr, pmap),
  int = qtl2::find_peaks(trait_intscan_chr, pmap),
  gene = qtl2::find_peaks(trait_genescan_chr, pmap),
  isoform = qtl2::find_peaks(trait_isoformscan_chr, pmap),
  gene_iso = qtl2::find_peaks(trait_gene_isoscan_chr, pmap),
  .id = "scan")
```

```{r}
plot(trait_intscan_chr, pmap, col = "red",
     xlim = trait_peaks$pos + c(-1,1))
plot(trait_genescan_chr, pmap, col = "blue", add = TRUE,
     xlim = trait_peaks$pos + c(-1,1))
plot(trait_isoformscan_chr, pmap, col = "purple", add = TRUE,
     xlim = trait_peaks$pos + c(-1,1))
plot(trait_gene_isoscan_chr, pmap, col = "black", add = TRUE,
     xlim = trait_peaks$pos + c(-1,1))
plot(trait_addscan_chr, pmap, col = "green", add = TRUE,
     xlim = trait_peaks$pos + c(-1,1))
abline(v = trait_peaks$pos, col = "gray", lwd = 2)
```

### Plot and Analysis of Best SNP by Trait

Proceed here assuming that we are working with `SNP36`.
If not, get those.

```{r}
if(snp36) {
  asnppr <- qtl2::genoprob_to_snpprob(aprobs_chr, snpinfo)
} else {
  asnppr <- snppr
  snppr <- qtl2::genoprob_to_snpprob(
    subset(
      readRDS(file.path(datadir, "fst_genoprob36", "probs_fstindex.rds")),
        chr = chrname),
      snpinfo)
}
```

```{r}
asnp_df <- dim(asnppr[[1]])[2] - 1
asnp_geno <- matrix(NA, nrow(addcovar), length(snp_id) * asnp_df)
rownames(asnp_geno) <- rownames(addcovar)
if(asnp_df == 1) {
  colnames(asnp_geno) <- snp_id
} else {
  colnames(asnp_geno) <- paste(rep(snp_id, each = 2), seq(asnp_df), sep = "_")
}
asnp_geno[dimnames(asnppr[[1]])[[1]],] <- asnppr[[1]][, -(1 + asnp_df), snp_id]
```

Combine physio trait (`physio_trait`), gene trait (`gene_trait`)
and genotype (`snp_geno` or `asnp_geno`) along
with covariates. Some machinations to impute genotypes from genoprobs.

```{r}
dat <- as.data.frame(cbind(physio_trait, gene_trait, round(snp_geno, 3)))
names(dat) <- c(physioshort, genename, "snp1", "snp2")
dat$snp3 <- round(1 - dat$snp1 - dat$snp2, 3)
dat <- dplyr::bind_cols(dat, dplyr::select(Physio, Sex,Diet.name)) |>
  tidyr::unite(SexDiet, Sex, Diet.name, remove = FALSE)
dat$geno <- unlist(apply(dat[, paste0("snp",1:3)], 1, function(x) {
  y <- which.max(x)
  if(length(y) == 0) y <- NA
  if(!all(is.na(y)) & all(x[y] < 0.75)) y <- NA
  y
  }))
dat <- dplyr::mutate(dat,
                     geno = ifelse(is.na(geno), NA, c("A","H","B")[geno]),
                     geno = factor(geno, c("A","H","B")),
                     asnp = round(asnp_geno[,1], 3),
                     ageno = round(2 * asnp),
                     ageno = ifelse(abs(2 * asnp - ageno) > 0.25, NA, ageno),
                     ageno = ifelse(is.na(ageno), NA, c("A","H","B")[3 - ageno]),
                     ageno = factor(ageno, c("A","H","B")),
                     wave = Physio$Gen_Lit)
```

Examine disagreements of 2-level vs 3-level SNPs

```{r}
dat |> dplyr::filter(geno != ageno) |> dplyr::select(-(1:2))
```

Find missing data for 3-level SNPs.

```{r}
dat |> dplyr::filter(is.na(geno)) |> dplyr::select(-(1:2))
```

Fill in missing `geno` where possible using `ageno`.

```{r}
dat <- dplyr::mutate(dat,
                     geno = ifelse(is.na(geno), as.character(ageno), as.character(geno)),
                     geno = factor(geno, c("A","H","B")),
                     add = unclass(geno) - 2,
                     dom = abs(add))
```

Plot faceting by sex and diet, with colors for waves.

```{r warning = FALSE}
ggplot2::ggplot(dat |> dplyr::filter(!is.na(geno))) +
  ggplot2::aes(x = .data[[genename]], y = .data[[physioshort]]) +
  ggplot2::geom_jitter(aes(col = geno)) +
  ggplot2::geom_smooth(col = "black", method = "lm", se = FALSE,
                       formula = "y~x") +
  ggplot2::facet_grid(SexDiet ~ geno) +
  ggplot2::theme(legend.position = "none")
```

Formal linear model fit without gene or `r genename`.

```{r}
fit <- lm(formula(paste(physioshort, "~ wave + Sex * Diet.name")),
          dat)
```

```{r}
drop1(fit, fit, test = "F")
```

Formal linear model fit, allowing for interaction of `r genename` with sex and diet.

```{r}
fit <- lm(formula(paste(physioshort, "~",
                        "+ wave + Sex * Diet.name *", genename)),
          dat)
```

```{r}
drop1(fit, fit, test = "F")
```

Formal linear model fit with predictor `r genename`,
allowing for interaction of SNP (`geno`) with sex and diet.

```{r}
fit <- lm(formula(paste(physioshort, "~", genename,
                        "+ wave + Sex * Diet.name * geno")),
          dat)
```

Examine significance of `r genename` as first predictor.

```{r}
anova(fit)
```

Type III adjusted tests.

```{r}
drop1(fit,fit, test="F")
```

Formal analysis to examine `add`itive and `dom`inance effects.

```{r}
fit <- lm(formula(paste(physioshort, "~", genename,
                        "+ wave + Sex * Diet.name * (add + dom)")),
          dat)
```

```{r}
drop1(fit,fit, test="F")
```

```{r}
knitr::knit_exit()
```

# Mediation

```{r}
medscan <- FALSE # mediation scan is very slow!
source("R/mediate.R")
```

```{r}
if(medscan) ggplot2::autoplot(med_scan)
```

```{r}
summary(med_test) |> dplyr::select(id,chr,pos,triad,pvalue,LR_mediation,LR_CMST)
```

```{r}
med_effect <- intermediate::mediation_effect(med_test, "symbol")
ggplot2::autoplot(med_effect)
```
## qtl2pattern

scan1pattern does not pass on incovar yet.
It also does not seem to work properly.
It should use the SDP in snpinfo identified to index snp to impute SNP probabilities
for pattern. It may be set up to do this with pmap but have some sort of disconnect in code.

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

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

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

```{r}
scan_pat <- qtl2pattern::scan1pattern(snppr, physio_trait, kinship_chr, addcovar,
                                      map = pmap, patterns = patterns)
```


```{r}
autoplot(scan_pat$scan, snpinfo, patterns = "hilit")
```

```{r}
summary(scan_pat$scan, pmap, chr = chrname)
```

```{r}
qtl2::find_peaks()
```

```{r}
tmp <- snpinfo[!duplicated(snpinfo$snp_id), "pos"]
names(tmp) <- unique(snpinfo$snp_id)
tmp <- list("19" = tmp)
```

```{r}
qtl2pattern::ggplot_scan1pattern(scan_pat, tmp)
```