%\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{Analysis of association of gene expression with rare variants}
\author{VJ Carey et al}

\section{Data}

We will use the CEU thousand genomes data, and focus on chromosome 15.
<<getd>>=
library(ceu1kg)
if (!exists("ceu1kg.sml")) data(ceu1kg.sml)
c6 = ceu1kg.sml[6]
c6
c6 = c6[[1]]
@
We will characterize the distribution of MAF.
<<lkdD>>=
csum = col.summary(c6)
dim(csum)
@

We have metadata on the locations of the chr6 SNP:
<<getm>>=
data(ceu1kgMeta_6)
ceu1kgMeta_6 = updateObject(ceu1kgMeta_6)
ceu1kgMeta_6[1:2]
all.equal(names(ceu1kgMeta_6), colnames(c6))
@

We will use MAF = 2.5\% as the threshold for rare variant declaration.
<<lkr>>=
isRare = which(csum[,"MAF"] > 0 & csum[,"MAF"] <= 0.025)
rareSNP_6 = c6[, isRare]
rareLocs_6 = ceu1kgMeta_6[ colnames(rareSNP_6) ]
@
The exons in chromosome 6 possessing rare variants are found via
<<lkfo,keep.source=TRUE>>=
data(ex6)  # excerpt from hg19 transcriptDb via exons(), saved for GGtools demos
exwr = findOverlaps( rareLocs_6, ex6 )
dim(exwr@matchMatrix)
@
The frequency of rare variants within exons is
<<lkfrr>>=
tt = table(table(exwr@matchMatrix[,2]))
tt
@
We can check these claims:
<<lkcla,keep.source=TRUE>>=
numr = max(as.numeric(names(tt))) # only names of the table is relevant
ind = as.numeric(names(which(table(exwr@matchMatrix[,2])==numr)))[1]
snin = which(ranges(rareLocs_6) %in% ranges(ex6[ind]))
ex6[ind]
rareLocs_6[snin]
@
So we have a (somewhat longwinded) mechanism for isolating the exons with rare variants.
<<getex>>=
exWithRareInd = unique(exwr@matchMatrix[,2])
exWithRare = ex6[ exWithRareInd ]
@
Now we want to split up the variants by exons harboring them in an identifiable way.

<<getspl>>=
vine = findOverlaps(exWithRare, rareLocs_6) 
enames = elementMetadata(exWithRare)$exon_id[ vine@matchMatrix[,1] ]
snames = names(rareLocs_6)[ vine@matchMatrix[,2] ]
vbye = split(snames,enames)
@
This would be useful if we had exon-targeting probes for phenotype
measurement.  We don't, so we'll aggregate by 'genes'.  Linking the
exon ids to genes is tedious.  This code assumes availability of the
hg19 transcript database, so is not evaluated in this vignette.
<<getlink,eval=TRUE>>=
library(GenomicFeatures)
#hg19.txdb = makeTranscriptDbFromUCSC("hg19")
#saveFeatures(hg19.txdb, "hg19.txdb.sqlite")
hg19.txdb = loadFeatures("~/hg19.txdb.sqlite")
exb = exonsBy(hg19.txdb, "gene")
library(org.Hs.eg.db)
egG6 = get("6", revmap(org.Hs.egCHR)) 
exb6 = exb[ names(exb) %in% egG6 ]
exb6namesByGenes = lapply(exb6, function(x) elementMetadata(x)$exon_id)
exb6namesByGenes[1:10]
ebg6 = lapply(exb6namesByGenes, as.character)
vbyg = lapply(ebg6, function(x) unlist(vbye[x]))
table(lens <- sapply(vbyg,length))
vbyg = vbyg[ lens > 0 ]
@
We have now grouped rare SNP by genes.  We can subset
our expression plus genotype data to include only the
relevant chromosome and genes.
<<doex,keep.source=TRUE>>=
if (!exists("ceu1kg")) data(ceu1kg)
ceu6 = ceu1kg[ chrnum("chr6"), ] # focus on chr6 snp
library(illuminaHumanv1.db)
ginplay = names(vbyg)
probes = mget(ginplay, revmap(illuminaHumanv1ENTREZID))
probes = unique(unlist(probes))
fn = featureNames(ceu6)
ceu6 = ceu6[ probeId(intersect(fn, probes)), ]
ceu6
@
We like to recode the genotype data so that the individuals
homozygous for the common allele have code 0.
<<dor>>=
remap = function(x) {
        n2 = sum(x==2);
         if (n2 > sum(x==0)) {
            tmp = x;
            tmp[x==0] = 2;
            tmp[x==2] = 0;
            return(tmp)
         } else x
   }
gt = as(smList(ceu6)[[1]], "numeric")
gt = gt[, as.character(unlist(vbyg)) ]
gt = apply(gt, 2, remap)
gt[1:5,1:5]
table(apply(gt,2,sum))
@
Finally, we want to accumulate within genes, the count of
rare alleles in the rare SNPs, for each individual.
Our groups are formed using entrez gene identifiers; need
to map to probe ids:
<<domo>>=
en = mget(names(vbyg), revmap(illuminaHumanv1ENTREZID))
table(sapply(en, length))
en1 = sapply(en, "[", 1)
vbygi = vbyg
names(vbygi) = en1
nn = which(is.na(names(vbygi)))
vbygi = vbygi[-nn]
@
Now we simply form vectors of accumulated rare allele counts,
and exclude genes with no variation in counts over individuals.
<<getrc>>=
rc = lapply(vbygi, function(sn) apply(gt[,sn,drop=FALSE],1,sum))
rcs = sapply(rc, sum)
rc = rc[which(rcs>0)]
@
We can test the specific linear models as follows:
<<lklm>>=
lms = list()
gn = names(rc)
for (i in 1:length(rc)) {
  lms[[i]] = lm(exprs(ceu6)[gn[i],]~rc[[i]])
}
names(lms) = gn
plms = sapply(lms, function(x)summary(x)$coef[2,4])
top = sort(plms)[1:10]
top
@
Let's look at some 'good hits':
<<lkda,fig=TRUE>>=
par(mfrow=c(3,3))
gd = names(top)
for (i in 1:9)
 plot( exprs(ceu6)[gd[i], ] ~ rc[[gd[i]]], main=gd[i],
   ylab = get(gd[i], illuminaHumanv1SYMBOL), xlab="# rare SNP minor alleles in gene" ) 
@






\end{document}
