%\VignetteIndexEntry{ENCODE/hg18 track demo: DNaseI hypersensitivity}
%\VignetteDepends{GGtools, GGdata}
%\VignetteKeywords{Genetics}
%\VignettePackage{encoDnaseI}


%
% NOTE -- ONLY EDIT THE .Rnw FILE!!!  The .tex file is
% likely to be overwritten.
%
\documentclass[12pt]{article}

\usepackage{amsmath,pstricks}
\usepackage[authoryear,round]{natbib}
\usepackage{hyperref}


\textwidth=6.2in
\textheight=8.5in
%\parskip=.3cm
\oddsidemargin=.1in
\evensidemargin=.1in
\headheight=-.3in

\newcommand{\scscst}{\scriptscriptstyle}
\newcommand{\scst}{\scriptstyle}


\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\textit{#1}}}
\newcommand{\Rmethod}[1]{{\texttt{#1}}}
\newcommand{\Rfunarg}[1]{{\texttt{#1}}}
\newcommand{\Rclass}[1]{{\textit{#1}}}

\textwidth=6.2in

\bibliographystyle{plainnat} 
 
\begin{document}
%\setkeys{Gin}{width=0.55\textwidth}

\title{Using the DNaseI hypersensitivity data from encode
in R}
\author{VJ Carey}

\maketitle

\section{Introduction}

Annotation tracks from UCSC hg18 can be used with Bioconductor
to help establish genomic contexts of events or alterations.
The CD4-based hypersensitivity assays are collected in the
structure rawCD4 in package encoDnaseI:
<<lkc>>=
library(encoDnaseI)
data(rawCD4)
rawCD4
@

At present, we can subset the data by casting a chromosome
number:
<<lkd>>=
c19 = rawCD4[ chrnum(19) ]
c19
@

And we can get a trace of values along the chromosome:
<<lklk,fig=TRUE>>=
c19xy = getTrkXY(c19)
plot(c19xy)
@

\section{Coupling the DnaseI series to genetics of gene expression}

We would like to subset a racExSet from GGdata and look at
snps that are in regions of high DNaseI sensitivity.
Some infrastructure to help with this is:

<<doinf>>=
clipSnps = function( rexset, chrmeta, lo, hi ) {
 allp = pos(chrmeta)
 ok = allp >= lo & allp <= hi
 rid = names(allp[ok])
 rexset[ snpID(rid), ]
}

rangeX = function(htrk) {
 range(getTrkXY(htrk)$x)
}
@

So we get the information on expression and SNPs in chr19 and
filter:
<<dogg>>=
library(GGtools)
library(GGdata)
data(chr19GGceuRMA)
data(chr19meta)
rs19 = rangeX(c19)
c19f = clipSnps( chr19GGceuRMA, chr19meta, rs19[1], rs19[2] )
c19f
@

A gene-specific screen can be computed as follows:
<<lkmxi1,fig=TRUE>>=
smxi1 = snpScreen(c19f, chr19meta, genesym("MXI1"), ~., fastAGMfitter )
plot_mlp(smxi1, chr19meta, "10")
@

We'd like to look at the SNP screen results juxtaposed with the
DnaseI results.

<<doj,fig=TRUE>>=
#juxtaPlot = function( trk, ssr ) {
# sy = abs(ssr$trat)
# sx = ssr@locs
# txy = getTrkXY( trk )
# df = data.frame( y = sy, x = sx, type="absT" )
# df = rbind(df, data.frame(y = txy$y, x = txy$x, type = "dnaseI"))
# df$type = factor(df$type, levels=c("dnaseI", "absT"))
# require(lattice)
# xyplot( y~x| type, data=df, layout=c(1,2), 
#    scales=list(relation=list(y="free")))
#}
print(juxtaPlot( c19, smxi1 ))
@
 
Another example:

<<donex,fig=TRUE>>=
sOSR2 = snpScreen(c19f, chr19meta, genesym("OSR2"), ~., fastAGMfitter )
print(juxtaPlot( c19, sOSR2 ))
@
 
We can score the highly associated snps for closeness to a highly DnaseI sensitive
region using ALICOR:

<<dow>>=
ALICOR( sOSR2, c19 )
ALICOR( smxi1, c19 )
@

To do this in the context of a transcriptome wide screen,
we filter the genes in the expression set non-specifically.

The following code takes about 15 minutes to run on a
decent linux box, and is suppressed from the standard build.

<<donsf>>=
if (interactive()) {
if (!exists("mads")) mads = apply(exprs(c19f), 1, mad)
if (interactive()) fn = featureNames(c19f)[ which(mads > quantile(mads, .6)) ]
if (!interactive()) fn = featureNames(c19f)[ which(mads > quantile(mads, .97)) ]
n19 = c19f[ exFeatID(fn), ]
if (file.exists("tw19.rda")) load("tw19.rda")
if (!exists("tw19")) tw19 = twSnpScreen( n19, chr19meta, ~., fastAGMfitter)
if (!file.exists("tw19.rda")) save(tw19, file="tw19.rda")
if (file.exists("allscor.rda")) load("allscor.rda")
if (!exists("allscor")) allscor = sapply(tw19, function(x) {if (inherits(x, "try-error")) return(NA) else return(ALICOR(x, c19))})
if (!file.exists("allscor.rda")) save(allscor, file="allscor.rda")
}

@

With these scores, we can find gene-snp combinations for which
association is at least partly synchronized with DHS.  Algorithms 
for systematically assessing synchronicity are in development.


\end{document}
