rs13107325

rs13107325.rmd

---
title: Analysis of rs13107325 heterogeneity in BMI - osteoarthritis MR analysis
author: Gibran Hemani
date: "`r Sys.time()"
---

```{r}
library(tryx)
library(TwoSampleMR)
library(ieugwasr)
library(ggplot2)
library(RadialMR)
library(rms2)
library(mr.ios)
```

Questions: Is rs13107325 having an influence on osteoarthritis via a pathway other than BMI

Analyses:

1. Heterogeneity analysis (does it contribute substantially as an outlier)
2. Index of suspicion (does it disproportionately associate with other traits)
3. rms2 (are there non-BMI pathways to the outcome)
4. tryx (is outlier status explained by pleiotropic pathways)

Other things

1. Does it have different eQTL profiles in different tissuesq / cell lines? 
2. Check individual level data


## Heterogeneity using ieu-a-2

```{r}
dat <- make_dat("ieu-a-2", "ebi-a-GCST007092")
datr <- dat_to_RadialMR(dat)[[1]]
resr <- ivw_radial(datr)
resr
```

```{r}
resr$outliers %>%
  mutate(p.value.adj = p.value * nrow(datr))
```


```{r}
plot_radial(resr)
```


Use the following GWASs:

```{r}
gwasinfo("ieu-b-40") %>% str()
```

```{r}
gwasinfo("ebi-a-GCST007092") %>% str()
```

Generate summary set

```{r}
dat <- make_dat("ieu-b-40", "ebi-a-GCST007092")
str(dat)
```

Check that rs13107325 is present among the instruments

```{r}
"rs13107325" %in% dat$SNP
```

Perform MR

```{r}
mr(dat)
```

Heterogeneity is high

```{r}
mr_heterogeneity(dat)
```

Is rs13107325 an outlier?

```{r}
res <- mr(dat, method="mr_ivw")
mr_scatter_plot(res, dat)[[1]] +
geom_label(data=subset(dat, SNP == "rs13107325"), aes(x=beta.exposure, y=beta.outcome, label=SNP))
```

Test outlier status more formally

```{r}
datr <- dat_to_RadialMR(dat)
resr <- ivw_radial(datr[[1]])
resr$outliers %>% filter(SNP == "rs13107325") %>%
  mutate(p.value.adj = p.value * nrow(dat))
```

This would not survive multiple testing


## Index of suspicion

```{r}
id_bg <- background_ids(id_exp="ieu-a-2", id_out="ebi-a-GCST007092", type="default")
# make background dataset
d <- make_dat("ieu-a-2", "ebi-a-GCST007092")
bg_dat <- make_background(exp = d, id_bg = id_bg)
```

```{r}
ios_dat <- ios(exp=d, bg=bg_dat)
ios_dat %>%
  tidyr::pivot_longer(!SNP) %>%
  tidyr::separate(name, sep="_", into=c("method", "dist")) %>%
  ggplot(., aes(y=value, x=SNP)) +
  geom_point(aes(colour=dist)) +
  facet_grid(method ~ ., scale="free_y") +
  theme(axis.text.x=element_text(angle=90))
```


## Use TRYX to analyse the potential pleiotropic paths

```{r}
tr <- Tryx$new(dat)
tr$get_outliers()
# Override to just use rs13107325
tr$output$outliers <- "rs13107325"

gi <- gwasinfo()
idlist <- gi$id

tr$set_candidate_traits(idlist)
tr$scan()
tr$output$sig_search$outcome

tr$output$sig_search <- tr$output$sig_search[c(1, 2, 3, 4, 5, 6, 15, 16, 22, 33:56, 58:67, 69, 72, 77, 78, 80, 82:86, 91, 98:126, 128:135, 140, 141, 144, 145, 153:178, 183:196, 198:201), ]

tr$output$sig_search <- tr$output$sig_search %>% 
	filter(!duplicated(outcome))

tr$output$sig_search$outcome
x <- tr
x$extractions()
x$harmonise()
x$mr()
x$tryx.sig()
x$analyse()



x$analyse.mv(lasso=FALSE)
x$output$analysis$adj_full %>% select(candidate, adj.qi, )

x$output$candidate_outcome_mr[x$output$candidate_outcome_mr$sig, ]

x$output$analysis$adj_full$id.exposure
```
Variability of qi

```{r}
x$output$analysis$adj_full %>% 
  ggplot(., aes(x=candidate, y=adj.qi)) +
  geom_point() +
  geom_hline(aes(yintercept=qi)) +
  theme(axis.text.x=element_text(angle=90, hjust=1))

```


```{r}
idlist <- c("ukb-d-30600_raw", "ukb-b-1489", "met-d-Total_CE", "ieu-b-110", "ebi-a-GCST90002412", "ukb-d-30690_irnt", "ukb-d-30600_irnt", "ukb-b-2393", "ukb-b-16489", "ukb-a-397", "ukb-a-9", "ubm-a-161", "ubm-a-547", "met-d-Sphingomyelins", "ieu-a-1239", "ukb-d-1468_4", "ukb-b-4424", "ubm-a-3041", "ebi-a-GCST006572", "ukb-d-30690_raw", "ukb-d-30680_raw", "ukb-d-30680_irnt")
expdat <- mv_extract_exposures(idlist)
od <- extract_outcome_data(expdat$SNP, "ebi-a-GCST007092")
d <- mv_harmonise_data(expdat, od)
mv_multiple(d)$result %>% select(exposure, pval)

ed <- extract_instruments(idlist)

ed %>% group_by(id.exposure) %>%
	summarise(n()) %>%
	as.data.frame()
snp <- ed %>% arrange(pval.exposure) %>% filter(!duplicated(SNP)) %>% select(SNP, pval.exposure) %>% clump_data()

od <- extract_outcome_data(snp$SNP, idlist)
table(od$id.outcome)

```


## RMS2 analysis

```{r}


```