\documentclass{article}

\usepackage{fullpage}
\usepackage{hyperref}
\marginparwidth 0pt
\oddsidemargin 0pt
\evensidemargin 0pt
\topmargin 0pt
\textwidth 16cm
\textheight 21cm

%\VignetteIndexEntry{Genome Reconstruction of Diversity Outbred Mice}

\usepackage{Sweave}

\begin{document}

\SweaveOpts{concordance=TRUE}

\title{Genome Reconstruction of Diversity Outbred Mice}
\author{Daniel M. Gatti}
\date{03 December 2013}
\maketitle

\section{Introduction}

\href{http://jaxmice.jax.org/strain/009376.html}{Diversity Outbred} (DO) mice are dervied from eight inbred founders: 
  \href{http://jaxmice.jax.org/strain/000646.html}{A/J},
  \href{http://jaxmice.jax.org/strain/000664.html}{C57BL/6J},
  \href{http://jaxmice.jax.org/strain/002448.html}{129S1/SvImJ},
  \href{http://jaxmice.jax.org/strain/001976.html}{NOD/ShiLtJ},
  \href{http://jaxmice.jax.org/strain/002105.html}{NZO/HlLtJ},
  \href{http://jaxmice.jax.org/strain/000928.html}{CAST/EiJ},
  \href{http://jaxmice.jax.org/strain/003715.html}{PWK/PhJ},
  \href{http://jaxmice.jax.org/strain/001145.html}{WSB/EiJ}.
DO genomes are a heterozygous mixture of these eight founders. Their genomes can be 
reconstructed using a hidden Markov model (HMM) in which the hidden states are the 36 possible 
unphased genotypes (8 homozygotes and 28 heterozygotes).  In order to estimate the genotype
probabilities, the HMM requires input data that allows it to estimate which founders contributed
to each sample at each marker.

Note that you will need an internet connection to run this vignette.

\section{Data Setup}

Below, we load in data from an example data package. The files contain X and Y intensities and allele calls.

<<label=setup>>=
library(DOQTL)
wd = tempdir()
library(MUGAExampleData)
data(x)
data(y)
data(geno)
data(pheno)
@

\section{Intensity Based Genome Reconstruction}

Genome reconstruction can be performed using the X and Y allele intensities from the MUGA or 
MegaMUGA. We have found that, using the MUGA or MegaMUGA in the DO, array intensities produce 
better genome reconstructions than the allele calls. The reason for this is that the 
intensities can contain more than three genotype groups and this aids in narrowing down the 
genotype categories.

When you have allele intensities, the overall X and Y chromosome intensities may be used
to assign the sex of each mouse. This is a good quality control check and these results
should be compared with the sex recorded during your experiment.

<<label=predict_sex,include=F>>=
load(url("ftp://ftp.jax.org/MUGA/muga_snps.Rdata"))
sex = sex.predict(x = x, y = y, snps = muga_snps, plot = T)
@

\begin{figure}
\begin{center}
<<label=fig1,fig=TRUE,echo=FALSE>>=
<<predict_sex>>
@
\end{center}
\caption{Mean X and Y chromosome intensity for DO samples, colored by sex. The females (black) have a higher mean X intensity than the males (red). Conversely, the males have higher mean Y chromosome intensity. The single female at (0.63, 0.43) may be an XO female. }
\label{fig:predict_sex}
\end{figure}

When DO mice are genotyped using the MUGA or MegaMUGA, DOQTL has access to the SNP locations and 
founder genotype data. You must supply the sample intensities, sex and DO outbreeding generation
in the \texttt{data} argument. 

data is a named list containing:
data\begin{description}
  \item{x: numeric matrix containing X allele intensities with samples in rows and markers in columns. Sample names in rownames and SNP IDs in colnames.}
  \item{y: numeric matrix containing Y allele intensities with samples in rows and markers in columns. Sample names in rownames and SNP IDs in colnames.}
  \item{sex: character vector containing either \texttt{M} or \texttt{F}. Sample names in names.}
  \item{gen: character vector containing 'DO' followed by teh DO outbreeding generation, with no space between them (i.e. DO4). Sample names in names.}
\end{description}

The DO outbreding generation should be requested when ordering mice.

<<label=calc_genoprob_intensity,fig=FALSE,results=hide>>=
gen = paste("DO", gsub("^G|L[12]$", "", pheno$Gen), sep = "")
names(gen) = rownames(pheno)
gen = gen[names(gen) %in% names(sex)]
gen = gen[match(names(sex), names(gen))]
stopifnot(all(rownames(x) == names(sex)))
stopifnot(all(rownames(x) == names(gen)))
data = list(x = x, y = y, sex = sex, gen = gen)
rm(x, y, sex, gen)
gc()
calc.genoprob(data = data, chr = c(1:19, "X"), output.dir = wd, array = "muga")
@

\texttt{calc.genoprob} uses the X and Y intensity data to reconstruct the DO genomes in terms 
of the 8 founder haplotypes. It then writes out one file for each sample that contains the 36
estimated genotype probabilities at each marker. It also writes out the model parameters and 
the 8 founder haplotype contributions for each sample at each locus.

\vspace{5 mm}

The output directory will now contain a single file for each sample containing the genotype probabilities at each marker. It will also contain a file with founder haplotype contributions for all samples and all markers called \texttt{model.probs.Rdata}. This is the file that will be used for QTL mapping.

<<label=intensity_genoplot,include=FALSE>>=
load(paste(wd, "F15.genotype.probs.Rdata", sep = "/"))
plot(x = prsmth, snps = muga_snps, main = "F15")
@

\begin{figure}
\begin{center}
<<label=fig2,fig=TRUE,echo=FALSE>>=
<<intensity_genoplot>>
@
\end{center}
\caption{Genome reconstruction of DO mouse using array intensities. }
\label{fig:intensity_genoplot}
\end{figure}

We can summarize the recombinations using \texttt{summarize.genotype.transitions()}. 

<<label=summarize_intensity,results=hide>>=
recomb = summarize.genotype.transitions(path = wd, snps = muga_snps)
@

This uses the maximum marginal probability at each locus to determine the recombinations in each sample. The function summarizes the recombination locations and genotypes for each sample.

<<>>=
head(recomb[[1]])
@

We can look at the average number of recombinations per sample.

<<>>=
mean(sapply(recomb, nrow))
@

\section{Allele Call Based Genome Reconstruction}

When DO mice are genotyped using the MUGA or MegaMUGA, DOQTL has access to the SNP locations and 
founder genotype data. You must supply the sample intensities, sex and DO outbreeding generation
in the \texttt{data} argument. 

data is a named list containing:
data\begin{description}
  \item{geno: character matrix containing allele calls with samples in rows and markers in columns. Sample names in rownames and SNP IDs in colnames.}
  \item{sex: character vector containing either \texttt{M} or \texttt{F}. Sample names in names.}
  \item{gen: character vector containing 'DO' followed by the DO outbreeding generation, with no space between them (i.e. DO4). Sample names in names.}
\end{description}

The DO outbreding generation should be requested when ordering mice.

<<label=calc_genoprob_allele,fig=FALSE,results=hide>>=
sex = as.character(pheno$Sex)
names(sex) = rownames(pheno)
gen = paste("DO", gsub("^G|L[12]$", "", pheno$Gen), sep = "")
names(gen) = rownames(geno)
sex = sex[names(sex) %in% rownames(geno)]
gen = gen[names(gen) %in% rownames(geno)]
data = list(geno = geno, sex = sex, gen = gen)
calc.genoprob(data = data, chr = c(1:19, "X"), output.dir = wd, array = "muga", 
              method = "allele")
@

The output directory will now contain a single file for each sample containing the genotype probabilities at each marker. It will also contain a file with founder haplotype contributions for all samples and all markers called \texttt{model.probs.Rdata}. This is the file that will be used for QTL mapping.

<<label=allele_genoplot,include=FALSE>>=
load(url("ftp://ftp.jax.org/MUGA/muga_snps.Rdata"))
load(paste(wd, "F15.genotype.probs.Rdata", sep = "/"))
plot(x = prsmth, snps = muga_snps, main = "F15")
@

\begin{figure}
\begin{center}
<<label=fig3,fig=TRUE,echo=FALSE>>=
<<allele_genoplot>>
@
\end{center}
\caption{Genome reconstruction of DO mouse using array allele calls. }
\label{fig:allele_genoplot}
\end{figure}

Note that the genotype resonstructions look quite different. We can summarize the number of recombinations as above.

<<label=summarize_intensity,results=hide>>=
recomb = summarize.genotype.transitions(path = wd, snps = muga_snps)
@

<<>>=
mean(sapply(recomb, nrow))
@

Whereas there were an average of 250 recombinations per sample above, now there are 133. This is consistent with the observations of several groups that the array intensities on the MUGA provide better genome reconstructions than that allele calls.

\end{document}