-
Notifications
You must be signed in to change notification settings - Fork 0
/
README.Rmd
113 lines (83 loc) · 2.9 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
---
output: md_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
:warning: this package is under development - use with care
# `R/qtl2helper`
`qtl2helper` provides with some helper functions to work with objects created by several functions in the [`R/qtl2`](https://kbroman.org/qtl2/) package.
To install the package run:
```{r, eval=FALSE}
# install the remotes package if not installed already
if (!("remotes" %in% rownames(installed.packages()))){
install.packages("remotes")
}
# install the package directly from github
remotes::install_github("tavareshugo/qtl2helper")
```
## Quick overview
Let's start with an example dataset from the [`qtl2` user guide](https://kbroman.org/qtl2/assets/vignettes/user_guide.html):
```{r}
library(qtl2)
library(qtl2helper)
# read example data
iron <- read_cross2( system.file("extdata", "iron.zip", package="qtl2") )
map <- insert_pseudomarkers(iron$gmap, step=1)
pr <- calc_genoprob(iron, map, error_prob=0.002)
Xcovar <- get_x_covar(iron)
# run QTL scan
out <- scan1(pr, iron$pheno, Xcovar=Xcovar)
# get QTL effects
c2eff <- scan1coef(pr[,"2"], iron$pheno[,"liver"])
```
`qtl2helper` provides methods for the `tidy` function (from the [`broom`]() package), which convert several `R/qtl2` objects to (long) data frames.
```{r}
# convert scan1 object to a tibble
tidy(out)
```
You can provide a _map_ of markers to `tidy()` to get marker positions in the output:
```{r}
# convert scan1 object to a tibble
tidy(out, map)
```
This makes it ideal to do further plotting, for example with `ggplot2`:
```{r}
library(ggplot2)
out_tbl <- tidy(out, map)
ggplot(out_tbl, aes(x = pos, y = LOD, colour = pheno)) +
geom_line() +
facet_grid(. ~ chrom, scales = "free_x", space = "free_x") +
theme_classic() +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
labs(x = "Chromosome", y = "LOD score", colour = "Tissue")
```
Here's the same idea applied to the estimated QTL effects:
```{r}
# tidy the scan1coef object
c2eff_tbl <- tidy(c2eff, map)
# plot
ggplot(c2eff_tbl, aes(x = pos, y = estimate, colour = coef)) +
geom_line() +
labs(x = "Chromosome 2 position", y = "QTL effects", colour = "Genotype")
```
Finally, here's an example of tidying permutations to obtain a tibble of thresholds:
```{r}
# run permutation scans
operm <- scan1perm(pr, iron$pheno, Xcovar=Xcovar, n_perm=10)
# tidy output
operm_tbl <- tidy(operm)
operm_tbl
```
These can now be added to the scan plot, for example using `geom_hline()`:
```{r}
ggplot(out_tbl, aes(x = pos, y = LOD, colour = pheno)) +
geom_line() +
geom_hline(data = operm_tbl,
aes(yintercept = threshold, colour = pheno),
linetype = "dotted") +
facet_grid(. ~ chrom, scales = "free_x", space = "free_x") +
theme_classic() +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
labs(x = "Chromosome", y = "LOD score", colour = "Tissue")
```