%\VignetteIndexEntry{Notes for eSet developers}
%\VignetteDepends{}
%\VignetteKeywords{Expression Analysis}
%\VignettePackage{Biobase}


%
% 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{Comprehensive search for trans-eQTL, multicohort case}
\author{VJ Carey}
\maketitle

\section{Introduction}

The search is conducted as follows.  A chromosome $c$ is identified and
all $S_c$ SNP on it will be analyzed.  The expression probes are divided
up by chromosome.  The top $K$ eQTL scores for each SNP on chromosome $c$
are collected, consisting of $K$ largest $\chi^2$ statistics and identifiers
for the $K$ probes on chromosome 1 corresponding to each SNP.  Thus two $S_c \times K$ matrices
are formed on the basis of associations for SNP on chromosome $c$ and probes
for chromosome 1.  The process is repeated, but for probes on chromosome 2.
The top $K$ probes (over probes on chromosomes 1 and 2) for each SNP are then
retained.  Ultimately we save only $S_c \times K$ scores and identifiers representing a transcriptome-wide search.

When multiple cohorts are involved, the sum of association tests over
cohorts is used to represent SNP-probe association.

\section{Illustration}

We will use HapMap genotype data (phase 2, 4 million SNP per individual) for
CEU and YRI populations.
<<getd>>=
library(GGdata)
library(hmyriB36)
@
We can get a representation of expression and genotype data as follows:
<<lkd>>=
ceu1 = getSS("GGdata", "1", renameChrs="chr1")
ceu1
@

\section{One cohort}

We can conduct a limited trans search for one cohort as follows.
We need to be able to determine which combinations are 'cis' and exclude
these from our collections, so some distance information is required.

Here we collect two chromosomes worth of recently curated SNP locations.
<<lksn>>=
library(SNPlocs.Hsapiens.dbSNP.20101109)
if (file.exists("c1s.rda")) load("c1s.rda")
if (!exists("c1s")) c1s = getSNPlocs("ch1", as.GRanges=TRUE)
if ("ch1" %in% seqlevels(c1s)) seqlevels(c1s) = gsub("ch", "chr", seqlevels(c1s))
if (!file.exists("c1s.rda")) save(c1s,file="c1s.rda")

if (file.exists("c2s.rda")) load("c2s.rda")
if (!exists("c2s")) c2s = getSNPlocs("ch2", as.GRanges=TRUE)
if ("ch2" %in% seqlevels(c2s)) seqlevels(c2s) = gsub("ch", "chr", seqlevels(c2s))
if (!file.exists("c2s.rda")) save(c2s,file="c2s.rda")
@

Probe related locations are a bit more involved.  For
these datasets, illumina human v1 is the chip.
<<getglocs>>=
library(GenomicRanges)
library(illuminaHumanv1.db)
if (!exists("gra")) {
 if (file.exists("gra.rda")) load("gra.rda")
 else {
  allst = as.list(illuminaHumanv1CHRLOC)
  allen = as.list(illuminaHumanv1CHRLOCEND)
  pchrs = sapply(allst, function(x)names(x)[1])
  bad = which(sapply(pchrs,length)==0)
  if (length(bad)>0) {
   pchrs = pchrs[-bad]
   pchrs = sapply(pchrs, function(x)x)
   pchrs = paste("chr", pchrs, sep="")
   allst = allst[-bad]
   allen = allen[-bad]
   }
 st = sapply(allst, function(x) abs(x)[1])
 en = sapply(allen, function(x) abs(x)[1])
 gra = GRanges(seqnames=pchrs, IRanges(st,en))
 names(gra) = names(allst)
 gra = split(gra, seqnames(gra))
 lkg = getSS("GGdata", "5")
 fn5 = featureNames(lkg)
 gra = lapply(gra, function(x)x[ intersect(names(x), fn5) ])
 save(gra, file="gra.rda")
 }
}
@

Now we can invoke the basic procedure.
<<lkdt,eval=FALSE>>=
#library(multicore)
#options(cores=10)
options(verbose=TRUE)
options(error=recover)
if (!exists("t1")) {
 if (file.exists("t1.rda")) load("t1.rda")
 else {
   system("rm -rf tscochr1")
   t1 = transScores("GGdata", snpchr = "1", ~male, K=10, geneApply=lapply,
     renameChrs="chr1", probesToKeep=1:300, batchsize=500, snpRanges=c1s,
     geneRanges=gra[["chr1"]], genegran=5)
   save(t1, file="t1.rda")
   }
}
@
This will obtain scores for the 10 best probes for all SNP on chr1,
out of the 300 probes 'kept' for this simple illustration.
We now investigate the resulting structure, which will be characteristic
for any \texttt{transScores} run:
<<loo,echo=FALSE>>=
if (!exists("t1")) load("t1.rda")
<<lkstr>>=
names(t1)
dim(t1$scores)
dim(t1$inds)
length(t1$guniv)
length(t1$snpnames)
t1$scores[101:115,1:5]
t1$inds[101:115,1:5]
@
The scores are 100 times the $\chi^2(1)$ statistic computed for the additive genetic
model relating mean expression for gene identified in the \texttt{inds} element to
number of 'B' alleles in the SNP.
We will now verify one of the reports.  We see a reasonably high value in the
14th row of scores, so identify the SNP-gene combination:
<<diver,keep.source=TRUE>>=
t1$snpnames[114]
t1$guniv[242]  # corresponding inds entry
library(GGtools)
N1 = getSS("GGdata", "1")
m1 = gwSnpTests(probeId("GI_11037056-S")~male, N1, chrnum("1")) # fit model with 
                               # the identified gene
m1@.Data[[1]]["rs2465126",]    # use snpStats structures to recover test result
@
It will be noted that there are lots of monomorphic SNPs.  For these, the test
statistics are not computable, but to avoid handling special codes for NA, we just
impute a random $\chi^2(1)$.

\section{Multiple cohorts}

<<doit>>=
if (file.exists("mt1.rda")) load("mt1.rda")
if (!exists("mt1")) {
mt1 = mtransScores(c("GGdata", "hmyriB36"), snpchr = "1", 
     c(~male, ~male), K=10, geneApply=lapply,
     renameChrs="chr1", probesToKeep=1:300, batchsize=500, snpRanges=c1s,
     geneRanges=gra[["chr1"]], genegran=5)
save(mt1,file="mt1.rda")
}
@


%#transScores = function (smpack, snpchr = "chr1", rhs, K = 10, targdirpref = "tsco",
%#    geneApply = mclapply, chrnames = paste("chr", as.character(1:22), sep=""),
%#    geneRanges = NULL, snpRanges = NULL, radius = 2e+06, renameChrs=NULL, probesToKeep=NULL)
%#

\end{document}
