# R/qtl-related package demo

## Overview
This is a quick demo of R/qtl and related packages R/qtlhot, R/qtlnet and R/qtlyeast for workshops on causal networks. For information and tutorials on R and R/qtl, visit http://www.rqtl.org . Each package his its own vignettes.

## R/qtl

We first show some features of R/qtl. For more info, see http://rqtl.org/tutorials/rqtltour2.pdf . This shows a summary of the Brem and Kruglyak (2001) yeast dataset. 

```{r}
library(qtl)
library(qtlyeast)
summary(yeast.orf)
```
## R/qtl genetic map

```{r}
plot.map(yeast.orf)
```

## Examine one trait.

We now look at one particular trait, which is located on chr 15.

```{r}
tf.orf <- "YOL084W"
tf.gene <- yeast.annot[yeast.annot$orf == tf.orf, "gene"]
tf.chr <- yeast.annot[yeast.annot$orf == tf.orf, "chr"]

## Invoke QTL hotspot and causal pair library.
library(qtlhot, quietly = TRUE)
cand <- cis.cand.reg$cis.reg[cis.cand.reg$cis.reg$gene == tf.orf,]
cand.peak <- cand$peak.pos
cand.pos <- cis.cand.reg$cis.reg$phys.pos[cis.cand.reg$cis.reg$gene == tf.orf]
```

## Image map of chr `r tf.chr`.

We examine one particular open reading frame, `r tf.orf`, corresponding to gene, `r tf.gene`. Here is an image of the genotypes for chromosome 'r tf.chr' reordered by this ORF.

```{r}
cand.mar <- which.min(abs(cand.pos - pull.map(yeast.orf)[[tf.chr]]))
geno.image(yeast.orf, tf.chr, reorder = find.pheno(yeast.orf, tf.orf))
abline(v = cand.mar, lwd = 4, col = "green")
```

This ORF is located at `r cand.pos`cM; the nearest marker (number `r cand.mar`) is identified with a green line on the genotype image. Now we do a simple interval map scan, which finds the strongest evidence (LOD) for a QTL right at the gene location. The red line is the 0.1% permutation threshold.

## 1-D scan of genome

```{r}
yeast.orf <- calc.genoprob(yeast.orf, step = 2)
scan.orf <- scanone(yeast.orf, pheno.col = find.pheno(yeast.orf, tf.orf), method="hk")
plot(scan.orf)
```

## 1-D scan focused on `r tf.chr`

```{r}
plot(scan.orf, chr = tf.chr)
abline(v = cand.pos, lwd = 3, col = "green")
lod.thr <- c(summary(perm.orf, 0.001))
abline(h = lod.thr, lwd = 3, col = "red", lty = 2)
```

# R/qtlhot: hotspots

## R/qtlhot: hotspots

Now let's look at hotspots. The `vignette("qtlhot", "qtlhot")` walks through a simulated set of hotspots. For the yeast data, see `vignette("yeast_cmst", "qtlyeast")`. Basically, we use the saved high LOD scores (`highlod.orf`) and LOD permutation threshold (`lod.thr`) to find hotspot size across the genome (`hotsize.orf`). The summary shows several significant hotspots, with one of the largest on chr 'r tf.chr', right at `r tf.orf`. The red 5% threshold of 96 is for the number of traits with significant LOD at a position (Chaibub Neto et al. 2012).

```{r}
hotsize.orf <- hotsize(highlod.orf, lod.thr = lod.thr)
hot.thr <- 96
(hot.sum <- summary(hotsize.orf, threshold = hot.thr))
plot(hotsize.orf, chr = hot.sum$chr, bandcol = "gray70")
abline(h = hot.thr, lwd = 2, col = "red", lty = 2)
```
## Hotspot on chr `r tf.chr`
```{r}
plot(hotsize.orf, chr = tf.chr)
abline(v = cand.pos, lwd = 3, col = "green")
abline(h = hot.thr, lwd = 2, col = "red", lty = 2)
```

# R/qtlhot: causal pairs

## R/qtlhot: causal pairs

Find co-mapped traits and get joint-CMST p-values with Benjaminie-Hochberg adjustment (Chaibub Neto et al. 2013).

```{r}
cand.tf <- cand.reg[cand.reg$gene == tf.orf,, drop = FALSE]
comap.tf <- GetCoMappingTraits(highlod.orf, cand.tf)[[1]]
```

The knock-out experiments in Hughes et al. (2000) and Zhu et al. (2008) provide a database of experimentally validated targets, which are stored in `ko.list` in the `qtlyeast` package. Let's find which of the co-mapping traits are targets of the transcription factor `r tf.gene`, and do causal mapping tests for them.

```{r}
match.orf <- !is.na(match(comap.tf, ko.list[[tf.orf]]))
out.peak <- CMSTtests(yeast.orf, tf.orf, comap.tf[match.orf], Q.chr = tf.chr, Q.pos = cand.peak,
                 method="joint", penalty="bic")
causal.call <- function(x) {ordered(c("causal","react","indep","corr")[apply(x, 1, which.min)],
                                    c("causal","react","indep","corr"))}
target.peak <- data.frame(pval = apply(out.peak$pvals.j.BIC, 1, min),
                          call = causal.call(out.peak$pvals.j.BIC))
row.names(target.peak) <- factor(comap.tf[match.orf])
```

## Focus on Significant causal calls

```{r}
pval.thr <- 0.001
## Look at table by various significance levels.
tmp <- 1 + (target.peak$pval <= 0.05) + (target.peak$pval <= 0.005) + (target.peak$pval <= 0.001)
lev <- c("n.s.","0.05","0.005","0.001")
target.peak$signif <- ordered(lev[tmp], lev)
rm(tmp, lev)
table(target.peak$signif, target.peak$call)
## Now reduce to the more significant targets.
target.peak$gene <- yeast.annot[match(row.names(target.peak), yeast.annot$orf), "gene"]
target.peak <- target.peak[order(target.peak$pval),]
target.peak <- target.peak[target.peak$pval <= 0.005,]
target.peak
```

## 1-D scan of causal trait and targets

Only `r sum(target.peak$pval <= pval.thr)` targets had a highly significant causal call at the `r pval.thr` level. Here are their scans (in blue) along with the `r tf.orf` (in black).

```{r}
plot(scan.orf, chr = tf.chr, xlim = c(0,100))
targets <- row.names(target.peak)[target.peak$pval <= pval.thr]
scans <- scanone(yeast.orf, pheno.col = find.pheno(yeast.orf, targets), method="hk")
for(i in seq(targets))
  plot(scans, chr = tf.chr, lodcolumn = i, add = TRUE, col = "blue")
abline(v = cand.pos, lwd = 3, col = "green")
abline(h = lod.thr, lwd = 3, col = "red", lty = 2)
```

Thus, it appears that the `r tf.gene` is causal to several targets, suggesting the value of inferring a causal network.

# R/qtlnet: causal network

## R/qtlnet: causal network

Let's take the significant targets of `r tf.gene` and see if we can develop a causal network.

```{r}
net.set <- c(tf.orf, row.names(target.peak)[target.peak$pval <= pval.thr])
library(qtlnet, quietly = TRUE)
```
The following will create MCMC samples for the QTLnet, but it takes awhile. Instead we read already created data.
```
phm7.qtlnet <- mcmc.qtlnet(yeast.orf, pheno.col = find.pheno(yeast.orf, net.set),
                           threshold = lod.thr, nSamples = 1000, thinning = 20,
                           random.seed = 123, max.parents = 3)
save(phm7.qtlnet, file = "phm7.qtlnet.RData", compress = TRUE)
```

## Plot of causal network.

```{r}
load("phm7.qtlnet.RData")
phm7.graph <- igraph.qtlnet(phm7.qtlnet)
## Make causal trait yellow.
phm7.graph <- set.vertex.attribute(phm7.graph, "color", 1, "yellow")
plot(phm7.graph)
```
The plot randomly arranges the traits and QTL. Alternatively, use the interactive
```
tkplot(phm7.graph)
```
In accordance with the results of the causality tests, 'r tf.orf' shows up in the top of the transcriptional network, being causal for all other traits. Note the inclusion of other QTL.
