
%
% NOTE -- ONLY EDIT THE .Rnw FILE!!!  The .tex file is
% likely to be overwritten.
%

%\VignetteIndexEntry{genetics of expression lab, 2010}
%\VignetteDepends{}
%\VignetteKeywords{genetics, gene expression, RNA-seq}
%\VignettePackage{GGtools}

\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{Genetics of gene expression: 2010 lab}
\author{VJ Carey}
\maketitle

\tableofcontents

\clearpage
\section{eQTL concepts and discovery tools}

The basic concern in the lab is the relationship between
structural variation in DNA and variation in mRNA abundance.
DNA variants of interest are primarily SNP as identified
through
\begin{itemize}
\item direct genotyping in the Sanger sequencing paradigm 
(yielding HapMap phase II genotypes, for example)
\item array-based genotyping (yielding HapMap phase III)
\item NGS-based variant calling (as provided for 1000 genomes (1KG))
\item hybrids of array-based and imputed genotypes (imputation to
the 1KG panel)
\end{itemize}
mRNA variation is typically characterized using gene expression microarrays,
but RNA-seq can also be considered.

A helpful schematic indicating possible impacts of structural variation
on gene transcription is given in \citet{Williams:2007p21}; see Figure \ref{schem}.

\begin{figure}
\begin{center}
\includegraphics{williamsSchem}
\end{center}
\caption{A schematic illustrating various nonexclusive mechanisms by which DNA
variants can affect transcript abundance \citep{Williams:2007p21}.}
\label{schem}
\end{figure}

\clearpage

To elucidate impacts of DNA variations on regulation of gene transcription
we need to be able to
\begin{itemize}
\item represent mRNA abundance and SNP genotypes (both observed
and imputed) in a coordinated way
for many samples
\item perform suitably focused and calibrated tests of association between mRNA
levels and genotype (e.g., allele count for a SNP)
\item provide significance assessments and functional contexts for
observed associations.
\end{itemize}
On the basis of excellent algorithms and software
in the \textit{snpMatrix} package of Clayton and Leung \citep{clayton2007}), the
\textit{GGtools} package provides tools to address all of these concerns.

\subsection{Checking for informative variants for a specified gene}

The following is a simple illustration of a focused workflow.
<<wk1>>=
library(GGtools)
if (!exists("hmceuB36.2021")) hmceuB36.2021 <- getSS("GGtools", c("20", "21"))
f1 = gwSnpTests(genesym("CPNE1")~male, hmceuB36.2021, chrnum(20))
topSnps(f1)
@


These findings can be visualized in the context of the chromosome in a
so-called Manhattan plot; see Figure \ref{manh}.

\begin{figure}
\begin{center}
<<f2,fig=TRUE>>=
plot(f1)
@
\end{center}
\caption{Manhattan plot for chromosome 20 eQTL for CPNE1.  The $y$ axis
measures the significance of association of SNP allele copy number with
expression for this gene.  The association is measured by a $\chi^2(1)$
statistic transformed to a $p$-value, which is logarithmically transformed
and then multiplied by $-10$.}
\label{manh}
\end{figure}

\subsection{Components of the workflow; refinements}

The primary workhorse thus far is \texttt{gwSnpTests}.  We need to understand
its interface and its return values.
<<lkgws>>=
showMethods("gwSnpTests")
@
<<lkf1>>=
f1
class(f1)
getClass(class(f1))
@

\subsubsection{Data representation: \texttt{smlSet}}

The \texttt{smlSet} container provides coordinated management of expression, genotype,
phenotype, and other metadata in a fashion broadly consistent with the \texttt{ExpressionSet}
class.  If \texttt{X} is an \texttt{smlSet} instance then
\begin{itemize}
\item \texttt{exprs(X)} and \texttt{pData(X)} have familiar roles of extracting a matrix of
expression values, and a \texttt{data.frame} of sample-level data
\item \texttt{X[,S]} subsets the samples according to predicate \texttt{S}
\item \texttt{X[probeId(P),]} subsets the expression data to probes enumerated in predicate \texttt{P}
\item \texttt{X[chrnum(C),]} subsets the genotype data to chromosome \texttt{C}
\item \texttt{annotation(X)} returns the annotation package used for mapping expression probe identifiers
\end{itemize}
Clearly there are two kinds of feature in use here: expression probes and SNPs.  The ambiguous
handling of the \texttt{X[G,]} idiom has not posed a problem thus far.  As the structure gains wider
use, more features may be added to simplify filtering through bracket-based operations.

Operations targeted on the genotype data component include:
\begin{itemize}
\item \texttt{smList(X)} returns a list of \texttt{snp.matrix} objects, using a byte to encode
each genotype
\item \texttt{MAFfilter(X,...)} will remove SNP that fail to meet a minor allele frequency criterion
\item \texttt{GTFfilter(X,...)} will remove SNP that fail to meet a minimum genotype frequency criterion
\end{itemize}

The \textit{GGtools} package provides \texttt{hmceuB36.2021} as a modest-sized exemplar of the
\texttt{smlSet} structure, with HapMap phase II genotype data from chromosomes 20 and 21, and GENEVAR
expression data for 90 CEU individuals.

<<lkhm>>=
hmceuB36.2021
names(pData(hmceuB36.2021))
table(hmceuB36.2021$isFounder)
@

The expression data is managed exactly as in an \texttt{ExpressionSet} from \textit{Biobase}.
The genotype data has a concise representation.
It and its coercions are defined in the \textit{snpMatrix} package.
<<lkg,keep.source=TRUE>>=
s20 = smList(hmceuB36.2021)[["20"]]
s20
as(s20, "matrix")[1:4,1:4]  # show raw
as(s20, "character")[1:4,1:4]
as(s20, "numeric")[1:4,1:4]
@


\textbf{Question 1:}  Use \verb+plot_EvG+ to visualize the specific relationship between expression
and allele copy number for the best result for CPNE1 shown above.

\textbf{Question 2:}  Create a new \texttt{smlSet} \texttt{hmff} limited to CEU founders
and with genotypes having minimum genotype frequency at least 5\%.  

\textbf{Question 3:}  Test for eQTL for CPNE1 in the filtered structure,
saving the result in \texttt{f1f} for later conversions.  Create relevant visualizations.


<<sol1,eval=TRUE,echo=FALSE>>=
#plot_EvG(genesym("CPNE1"), rsid("rs17093026"), hmceuB36.2021[chrnum(20),])
hmff = hmceuB36.2021[, which(hmceuB36.2021$isFounder)]
hmfg = GTFfilter( hmff, .05 )
f1f = gwSnpTests(genesym("CPNE1")~male, hmfg, chrnum(20))
plot_EvG(genesym("CPNE1"), rsid("rs11698601"), hmfg[chrnum(20),])
@

\subsubsection{Analysis: \texttt{snp.rhs.tests}; permutation testing}

The \textit{snpMatrix} package supplied various utilities for general analysis of
genome-wide association studies.  \texttt{gwSnpTests} is just a wrapper for the \texttt{snp.rhs.tests}
function.  To see it in action,
<<dosrs>>=
SS = snp.rhs.tests(exprs(hmceuB36.2021)["GI_23397697-A",]~male, data=pData(hmceuB36.2021),
   snp.data=s20, family="gaussian") -> SS
class(SS)
length(names(SS))
sum(is.na(p.value(SS)))
SS[1:4,]
SS["rs6060535",]
@
This shows that for a given phenotypic response, large numbers of GWAS tests can be performed
rapidly.

\textbf{Question 4:}  Why do so many tests have p-value NA?

<<lk,eval=FALSE,echo=FALSE>>=
table(as(smList(hmceuB36.2021)[[1]][,10], "character"))
@

Which of our findings are significant in the context of so many tests?  We can permute the expression
data, keeping the genotype data fixed.  Here is one way to obtain a sample from the permutation distribution
of the smallest p-value for a specific gene.

<<dot,eval=FALSE>>=
Mtests = sapply(1:100,  function(x) topSnps(gwSnpTests(genesym("CPNE1")~male, 
        sms=permEx(hmceuB36.2021), cnum=chrnum(20)))[[1]][1])
@

\subsubsection{Context: \textit{GenomicRanges} and \textit{rtracklayer}}

A major response to the transition to high-throughput sequencing assays is the \textit{IRanges}
infrastructure.  The \textit{GenomicRanges} package defines the \texttt{GRanges} class, to which
we can convert results of an eQTL search:
<<lk2>>=
gr1 = as(f1, "GRanges")
gr1[1:4,]
gr2 = as(f1f, "GRanges")
@

\textbf{Question 5:}  Export these GRanges instances as wig files and add to a UCSC browser
session as custom tracks.  You could do this with \textit{rtracklayer} if your searchlist had
a particular form.  You should see something like Figure \ref{trks}; R Track 1 is the search
on all SNP; R Track 2 follows the restriction to minimum genotype frequency 5\%.

\begin{figure}
\includegraphics{full20browseX2}
\caption{Browser tracks associated with eQTL screens \texttt{f1} and \texttt{f1f}.}
\label{trks}
\end{figure}

\clearpage
\section{Comprehensive eQTL surveys: performance issues}

A comprehensive survey of eQTL would consider possible relationships between 20000 or so
genes and 8 million or more SNP (the latter is the approximate size of the 1000 genomes
SNP panel).  Creating, storing, and interrogating the family of 160 billion results is
a significant undertaking.  We now illustrate a strategy for conducting comprehensive
surveys.  

\subsection{A study of genes coresident on chromosome 20}

We will select 50 genes at random from chromosome 20 and compute
all tests for association with SNP on chromosome 20.

First we select the genes and filter the \texttt{smlSet}:
<<getg>>=
set.seed(1234)
library(illuminaHumanv1.db)
g20 = get("20", revmap(illuminaHumanv1CHR))
samp = sample(g20, size=50, replace=FALSE)
hlit = hmceuB36.2021[ probeId(samp), ]
@

If you are using a POSIX-compliant system, you can speed up
this process by distributing over cores.  We will be writing data
to a folder called \texttt{dem50}.
<<domc,eval=FALSE>>=
try(system("rm -rf dem50"))
if (try(require(multicore))) {
 unix.time(e1 <- eqtlTests(hlit, ~male, geneApply = mclapply, targdir="dem50" ))
}
@
If you do not have \textit{multicore}, alter the code to:
<<dom2,eval=TRUE>>=
 try(system("rm -rf dem50"))
 unix.time(e1 <- eqtlTests(hlit, ~male, geneApply = lapply, targdir="dem50" ))
@
The resulting object is:
<<aaa>>=
e1
class(e1)
getClass(class(e1))
@
We are using the \textit{ff} package to use disk to maintain extremely
voluminous outputs.
<<lkff>>=
e1@fflist[[1]][1:4,1:4]
e1@shortfac
@
The information shown above is not intended for public consumption.
Instead we perform focused interrogation using brackets and casting:
<<doi>>=
g1 = probeNames(e1)[1]
r1 = rownames(e1@fflist[[1]])[1]
e1[rsid(r1), probeId(g1)]
@
Typical use case: finding a collection of strongest associations:
<<lku>>=
topFeats(probeId(g1), mgr=e1, ffind=1, anno="illuminaHumanv1.db")
topFeats(rsid(r1), mgr=e1, ffind=1, anno="illuminaHumanv1.db")
topFeats(genesym("SRC"), mgr=e1, ffind=1, anno="illuminaHumanv1.db")
@

\subsection{Expanding to genes on multiple chromosomes}

The analysis conducted just above uses genes on chromosome 20.  We
repeat with some genes on chromosome 21.  First sample the genes:
<<getg21>>=
g21 = get("21", revmap(illuminaHumanv1CHR))
samp21 = sample(g21, size=50, replace=FALSE)
hlit21 = hmceuB36.2021[ probeId(samp21), ]
@
Now perform the tests, using a different target folder:
<<dot>>=
try(system("rm -rf dem50.21"))
unix.time(e2 <- eqtlTests(hlit21, ~male, geneApply = lapply, targdir="dem50.21" ))
e2
@
To facilitate interrogation of the two sets of results simultaneously, the
\texttt{eqtlTestsManager} instances need to be coordinated by a
\texttt{cisTransDirector} instance.
<<lkdir>>=
try(system("rm -rf d2021_probetab.ff"))
try(system("rm -rf d2021_snptab.ff"))
d1 = mkCisTransDirector(list(e1,e2), "d2021", "snptab", "probetab", "illuminaHumanv1.db" )
d1
d1[ "rs6090120", "GI_42734342-S" ]
@
You are now managing and interrogating 
\Sexpr{ sum(sapply(d1@mgrs, function(x) sapply(x@fflist,length))) }
test results.  By increasing
the numbers of genes and SNP held in the \texttt{smlSet} instances
passed to \texttt{eqtlTests} you can increase this.  Interfaces need
further refinement.

\clearpage

\subsection{Exercises}

\textbf{Question 6.}  You can use \texttt{as.ram} from the ff package to work directly with
data in \texttt{ff} arrays.  You have to be careful or you will create large
objects and defeat the purpose of \texttt{ff}.
Show that the gene/SNP combination shown here has a special status in
this series of tests.
<<getmax,echo=FALSE>>=
s = rep(NA,4)
g = rep(NA,4)
# we search the 2 chromosomes by 2 gene sets for highest score
asram = ff:::as.ram
s[1] = which.max(apply(asram(e1@fflist[[1]]),1,max))
g[1] = which.max(apply(asram(e1@fflist[[1]]),2,max))
s[2] = which.max(apply(asram(e1@fflist[[2]]),1,max))
g[2] = which.max(apply(asram(e1@fflist[[2]]),2,max))
s[3] = which.max(apply(asram(e2@fflist[[1]]),1,max))
g[3] = which.max(apply(asram(e2@fflist[[1]]),2,max))
s[4] = which.max(apply(asram(e2@fflist[[2]]),1,max))
g[4] = which.max(apply(asram(e2@fflist[[2]]),2,max))
ind = which.max( c(
   e1@fflist[[1]][s[1],g[1]],
   e1@fflist[[2]][s[2],g[2]],
   e2@fflist[[1]][s[3],g[3]],
   e2@fflist[[2]][s[4],g[4]] ))
if (ind %in% c(1,2)) {  # max is in e1
     pr = colnames(e1@fflist[[ind]])[g[ind]]
     rs = rownames(e1@fflist[[ind]])[s[ind]]
     } else {    # max is in e2
     pr = colnames(e2@fflist[[ind-2]])[g[ind]]
     rs = rownames(e2@fflist[[ind-2]])[s[ind]] 
     }
<<lkdd,fig=TRUE,echo=FALSE>>=
if (ind %in% c(1,2)) plot_EvG(probeId(pr), rsid(rs), hlit) else
  plot_EvG(probeId(pr), rsid(rs), hlit21)
if (ind %in% c(1,3)) topFeats(probeId(pr), mgr=d1, ffind=1, anno="illuminaHumanv1.db") else
  topFeats(probeId(pr), mgr=d1, ffind=2, anno="illuminaHumanv1.db")
@

\textbf{Question 7.}  Use \texttt{getGRanges} to obtain a Manhattan plot
for the gene identified in question 6, displayed in the browser.  You should see
something like

\includegraphics{maxscan}

<<lkgg,eval=FALSE,echo=FALSE>>=
ggg = getGRanges(e2, 1, "GI_25777639-A", "chr20", 
   getNamedLocs(slpack="SNPlocs.Hsapiens.dbSNP.20100427", chr="ch20"))
@

\clearpage
\section{Investigating allele-specific expression with 
RNA-seq data}

We consider a study of read mapping bias in RNA-seq
\citep{Degner:2009p1333}.  We will consider how to identify
SNPs located in exons, how to analyze allele frequencies
for such SNP, and how to check findings against existing sequencing data
and eQTL statistics for a single individual and a HapMap population.

\subsection{The data}

We will use a samtools pileup representation of RNA-seq data that
were aligned using MAQ and distributed by Degner and colleagues through GEO.
We confine attention to chromosome 6 for individual NA19238,
and use the ``unmasked'' alignment files.  Data from two sequencing
runs was combined using samtools merge, and then the pileup was
generated against hg18.
<<rpup>>=
library(GGtools)
library(Rsamtools)
pup238_6 = readPileup(system.file("pup/combn238_chr6.pup", package="GGtools"), 
    variant="SNP")
pup238_6[1:4,]
@
It will be useful to have the specific pileup of calls as well.
<<getc>>=
getPupCalls = function(pupfn) {
 tmp = strsplit(readLines(pupfn), "\t")
 sapply(tmp, "[", 9)
}
callp = getPupCalls( system.file("pup/combn238_chr6.pup", package="GGtools") )
elementMetadata(pup238_6)$callp = callp
pup238_6[144:148,]
@

\subsection{Filtering to ``coding SNP''; checking for de novo variants}

We have derived a \texttt{GRanges} instance with addresses
of exons on chromosome 6, using the \textit{GenomicFeatures}
facilities for extracting feature data from UCSC tables.
<<getd>>=
#hg18.txdb = loadFeatures(system.file("extdata/hg18.txdb.sqlite",
#  package="anno4NGS"))
#if (.Platform$OS.type == "windows") 
#  data(ex6) else 
library(GenomicFeatures)
data(ex6)
ex6[1:4,]
#ex6 = exons(hg18.txdb, vals=list(exon_chrom="chr6"))
exp_6 = pup238_6[ which( ranges(pup238_6) %in% ranges(ex6) ) , ]
@

How many of the SNP variants in this individual reside in exons?
Restrict the pileup data to these.
<<lkmsnp>>=
sum( ranges(pup238_6) %in% ranges(ex6) )
isExonic = which( ranges(pup238_6) %in% ranges(ex6) )
pup238_6x = pup238_6[ isExonic, ]
@

How many of these exonic variants are already catalogued in
%the June 2008 dnSNP database?
the May 2009 dnSNP database?
<<lkdb>>=
library(SNPlocs.Hsapiens.dbSNP.20100427)
s6 = getSNPlocs("ch6")
s6[1:5,]
knownLocs = IRanges(s6$loc, s6$loc)
indbsnp = ranges(pup238_6x) %in% knownLocs
sum(indbsnp)
@

What are the heterozygosity frequencies at known and potentially novel SNP
locations?

<<lkhz>>=
p6xknown = pup238_6x[ which( indbsnp ), ] 
p6xnovel = pup238_6x[ -which( indbsnp ), ]
nKnown = length(p6xknown)
nhetKnown = sum( !(elementMetadata(p6xknown)$consensusBase %in% c("A", "C", "G", "T") ) )
nhetKnown/nKnown
nNovel = length(p6xnovel)
nhetNovel = sum( !(elementMetadata(p6xnovel)$consensusBase %in% c("A", "C", "G", "T") ) )
nhetNovel/nNovel
@

\subsection{Assessing allelic imbalance in transcripts harboring SNPs}

What exactly are we looking for?  A guide is available through the table
published by Degner et al.  
<<lktab>>=
data(degnerASE01)
antab = with(degnerASE01, degnerASE01[ chr=="chr6" & indiv=="GM19238", 1:8] )
antab
@
This shows 5 loci related to 4 genes for which there is evidence of
allele-specific expression, in that the individual is heterozygous
at the locus, but the transcript abundance measures are skewed towards
one of the two alleles present.

To check for this on the basis of our pileup data, we have a little
tabulation function:
<<mkta>>=
tabCalls = function (pup, ind) 
{
    ac = as.character
    empup = elementMetadata(pup)
    ref = ac(empup$referenceBase[ind])
    maqcall = ac(empup$consensusBase[ind])
    puptag = gsub("\\^.", "", empup$callp[ind])
    list(ref = ref, maqcall = maqcall, calls = table(strsplit(puptag, 
        "")))
}

@

In the following, we
tind overlaps between our exonic loci and the tabulated locations, and
tabulate the available calls.  First we build a little utility that
returns nucleotides given an ambiguous code.
<<zzz>>=
library(Biostrings) # has IUPAC_CODE_MAP
decodeIU = function(x) {
 if(length(x)>1) warning("only handling first entry")
 strsplit(IUPAC_CODE_MAP[toupper(x)], "")[[1]]
}
@

Now we identify those locations from the Degner table
that coincide with our pileup-based ranges.  There is
an idiosyncratic but intuitive matrix result wrapped
in S4.
<<lkra>>=
ol = findOverlaps(IRanges(antab$loc, width=1), ranges(pup238_6x) )
ol
@

In the following, we generate a little report to summarize
the comparisons.
<<lkcomp,keep.source=TRUE>>=
for (i in 1:nrow(ol@matchMatrix)) {  # ONEOFF
 cat("query", i, "\n")
 print(t1 <- antab[ ol@matchMatrix[i, "query"], ])
 t2 <- tabCalls( pup238_6x, ol@matchMatrix[i, "subject"])
 print(t2[["calls"]])
 dpref = t1[1,2]/(t1[1,2]+t1[1,3])
 cat("Degner refFreq =", round(dpref,4),"\n")
 dpuprefinds = which(toupper(names(t2[["calls"]])) %in% c(",", "."))
 targcodes = decodeIU( t2[["maqcall"]] )
 altcode = setdiff(targcodes, t2[["ref"]])
 dpupaltinds = which(toupper(names(t2[["calls"]])) == altcode)
 dpuppref = sum(t2[["calls"]][dpuprefinds])/(sum(t2[["calls"]][dpuprefinds])+
    sum(t2[["calls"]][dpupaltinds]))
 cat("Pileup refFreq =", round(dpuppref,4),"\n")
 cat("---\n")
}
@
The correspondence is very close.  Unfortunately there is relatively low
coverage for the most extreme imbalances.

\subsection{Checking consistency of findings with GENEVAR expression arrays}

The hmyriB36 package of Bioconductor's experimental data archive includes
expression and genotyping data on 90 YRI individuals including NA19238.
The expression data are distributed by the GENEVAR project \citep{Stranger:2007p1342}.
<<lkd>>=
library(hmyriB36)
if (!exists("hmyriB36")) data(hmyriB36)
library(illuminaHumanv1.db)
pidVEGFA = get("VEGFA", revmap(illuminaHumanv1SYMBOL))
ve = exprs(hmyriB36)[ pidVEGFA, "NA19238"]
@
The following code can be used to check for relationship between allele
copy number for a SNP and array-based expression values:
<<lkgev,fig=TRUE>>=
plot_EvG(genesym("VEGFA"), rsid("rs3025040"), hmyriB36)
abline(h=ve)
@


\subsection{Exercises}

\indent
\textbf{Question 8.}  Check the array-based expression configurations
for the other four genes with potential allele-specific expression.

\textbf{Question 9.}
Write code that scans all the exonic SNP on chr6 and identifies
loci with potential imbalance.  You might modify the
code snippet marked \texttt{ONEOFF} above.  What statistical tests and
multiple testing corrections should you use?

\textbf{Question 10.}  Acquire the data on NA19239 and reproduce some of the key
entries of the Degner table.


\clearpage
\section{Appendix: samtools pileup format, from samtools distribution site}

{\small
\begin{verbatim}
Pileup Format

Pileup format is first used by Tony Cox and Zemin Ning at the Sanger 
Institute. It describes the base-pair information at each chromosomal 
position. This format facilitates SNP/indel calling and brief alignment 
viewing by eyes.

The pileup format has several variants. The default output by SAMtools looks like this:

      seq1 272 T 24  ,.$.....,,.,.,...,,,.,..^+. <<<+;<<<<<<<<<<<=<;<;7<&
      seq1 273 T 23  ,.....,,.,.,...,,,.,..A <<<;<<<<<<<<<3<=<<<;<<+
      seq1 274 T 23  ,.$....,,.,.,...,,,.,...    7<7;<;<<<<<<<<<=<;<;<<6
      seq1 275 A 23  ,$....,,.,.,...,,,.,...^l.  <+;9*<<<<<<<<<=<<:;<<<<
      seq1 276 G 22  ...T,,.,.,...,,,.,....  33;+<<7=7<<7<&<<1;<<6<
      seq1 277 T 22  ....,,.,.,.C.,,,.,..G.  +7<;<<<<<<<&<=<<:;<<&<
      seq1 278 G 23  ....,,.,.,...,,,.,....^k.   %38*<<;<7<<7<=<<<;<<<<<
      seq1 279 C 23  A..T,,.,.,...,,,.,..... ;75&<<<<<<<<<=<<<9<<:<<

where each line consists of chromosome, 1-based coordinate, reference base, the number 
of reads covering the site, read bases and base qualities. At the read base 
column, a dot stands for a match to the reference base on the forward strand, a 
comma for a match on the reverse strand, `ACGTN' for a mismatch on the forward 
strand and `acgtn' for a mismatch on the reverse strand. A pattern `\+[0-9]+[ACGTNacgtn]+' 
indicates there is an insertion between this reference position and the next reference 
position.  The length of the insertion is given by the integer in the pattern, followed 
by the inserted sequence. Here is an example of 2bp insertions on three reads:

      seq2 156 A 11  .$......+2AG.+2AG.+2AGGG    <975;:<<<<<

Similarly, a pattern `-[0-9]+[ACGTNacgtn]+' represents a deletion from the reference. 
Here is an exmaple of a 4bp deletions from the reference, supported by two reads:

      seq3 200 A 20 ,,,,,..,.-4CACC.-4CACC....,.,,.^~. ==<<<<<<<<<<<::<;2<<

Also at the read base column, a symbol `^' marks the start of a read segment which is 
a contiguous subsequence 
on the read separated by `N/S/H' CIGAR operations. The ASCII of the character 
following `^' minus 33 gives the mapping quality. A symbol `$' marks the end of a read 
segment. Start and end markers of a read are largely inspired by Phil Green's CALF 
format. These markers make it possible to reconstruct the read sequences from pileup. 
 
\end{verbatim}
}


\section{Session information}

<<dos,results=tex>>=
toLatex(sessionInfo())
@

\clearpage
\bibliography{jsm2010}

\end{document}



