R/qtl-related package demo
========================================================

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 .

```{r}
library(qtl)
library(qtlyeast)
summary(yeast.orf)
plot.map(yeast.orf)
```

This shows a summary of the Brem and Kruglyak (2001) yeast dataset, with `r sum(nmar(yeast.orf))` markers and `r nphe(yeast.orf)` phenotypes. 

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

## Invoke QTL hotspot and causal pair library.
library(qtlhot)
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]
```

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

```{r}
cand.mar <- which.min(abs(cand.pos - pull.map(yeast.orf)[[2]]))
geno.image(yeast.orf, 2, reorder = find.pheno(yeast.orf, tf.orf))
abline(v = cand.mar, lwd = 3, 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 1% permutation threshold.

```{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)
plot(scan.orf, chr = 2)
abline(v = cand.pos, lwd = 3, col = "green")
lod.thr <- c(summary(perm.orf, 0.01))
abline(h = lod.thr, lwd = 3, col = "red", lty = 2)
```

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 2, 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)
plot(hotsize.orf, chr = 2)
abline(v = cand.pos, lwd = 3, col = "green")
abline(h = hot.thr, lwd = 2, col = "red", lty = 2)
```

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 [TransFac](http://www.gene-regulation.com/pub/databases.html) database provides transcription factor 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 = 2, 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])
table(target.peak$pval <= 0.1, target.peak$call)
target.peak$gene <- yeast.annot[match(row.names(target.peak), yeast.annot$orf), "gene"]
target.peak[order(target.peak$pval),]
```

Only `r sum(target.peak$pval <= 0.1)` targets had a significant causal call, and all of those were for the independent model (pleiotropy). Let's look only at the alleged targets of `r tf.gene`. We use color to identify the calls for these `r sum(match.orf)` ORFs.

```{r}
plot(scan.orf, chr = 2)
for(target in row.names(target.peak)) {
  tmp <- scanone(yeast.orf, pheno.col = find.pheno(yeast.orf, target), method="hk")
  if(max(tmp)$lod > lod.thr) {
    col <- "gray"
    if(target.peak[target, "pval"] < 0.1)
      col <- c("red", "blue", "green","gray")[unclass(target.peak[target, "call"])]
    plot(tmp, chr = 2, add = TRUE, col = col)    
  }
}
plot(scan.orf, chr = 2, add = TRUE)
abline(v = cand.pos, lwd = 3, col = "green")
abline(h = lod.thr, lwd = 3, col = "red", lty = 2)
```

So, it appears that the `r tf.gene` is pleiotropic to targets with strong QTL signals.

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 <= 0.1])
library(qtlnet)
```
The following will create MCMC samples for the QTLnet, but it takes awhile. Instead we read already created data.
```
amn1.qtlnet <- mcmc.qtlnet(yeast.orf, pheno.col = find.pheno(yeast.orf, net.set),
                           threshold = lod.thr, nSamples = 1000, thinning = 20,
                           random.seed = 92387475)
save(amn1.qtlnet, file = "amn1.qtlnet.RData", compress = TRUE)
```
```{r}
load("amn1.qtlnet.RData")
amn1.graph <- igraph.qtlnet(amn1.qtlnet)
plot(amn1.graph)
```
The plot randomly arranges the traits and QTL. Alternatively, use the interactive
```
tkplot(amn1.graph)
```
This apparently did not work well, as the `r tf.orf` is shown as directly causal to one target and reactive to two others. More work would be needed to study this more carefully.
