March 25, 2015
Principal Components Analysis (PCA)
and visualization
Goal: a low-dimensional summary.
The first step to data analysis: look at the data.
How to look at millions of loci (/variables/dimensions)?
(keyword: "dimensionality reduction")
Overview:
- examples from the literature
- how to do it: overview
- properties of PCA
- how to do it: hands-on
Genes mirror geography within Europe, Novembre et al 2008 
New insights into the Tyrolean Iceman's origin and phenotype as inferred by whole-genome sequencing, Keller et al 2011 
Genomic divergence in a ring species complex, Alcaide et al 2014 
Anatomy of PCA
Your data:
| Â | SNP1 | SNP2 | SNP3 | SNP4 | SNP5 | SNP6 | SNP7 | SNP8 | SNP9 | SNP10 | SNP11 | SNP12 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| I1 | C/C | C/C | T/A | A/A | T/T | C/C | C/A | T/T | G/G | T/T | C/C | C/C |
| I2 | C/C | C/C | A/A | A/A | T/T | C/C | A/A | T/T | G/C | C/T | C/C | C/C |
| I3 | C/C | C/C | A/A | A/A | C/C | C/C | C/C | T/T | G/G | C/C | C/C | C/C |
| I4 | C/C | C/C | T/A | A/A | C/C | C/C | C/C | T/T | C/C | T/T | C/C | C/C |
| I5 | C/C | T/C | A/A | A/A | C/C | C/C | C/C | T/T | G/C | T/T | C/C | C/C |
First: change letters to numbers.
Translate alleles to "proportion of reference alleles":
| Â | SNP1 | SNP2 | SNP3 | SNP4 | SNP5 | SNP6 | SNP7 | SNP8 | SNP9 | SNP10 | SNP11 | SNP12 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| I1 | 1 | 0 | 0.5 | 0 | 1 | 1 | 0.5 | 1 | 1 | 0 | 1 | 0 |
| I2 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 1 | 0.5 | 0.5 | 1 | 0 |
| I3 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 |
| I4 | 1 | 0 | 0.5 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 1 | 0 |
| I5 | 1 | 0.5 | 0 | 0 | 0 | 1 | 1 | 1 | 0.5 | 0 | 1 | 0 |
note: arbitrary choice of "reference" allele.
Second: compute a similarity matrix
Default choice:
Trim tightly linked loci (e.g. with plink)
- reduces influence of large linked blocks
Subtract locus means (\(2 \times{}\) allele frequencies)
- makes result independent of reference allele choices
Estimate covariance between individuals
- methods available for missing data
Second: compute a similarity matrix
e.g. the covariance matrix, obtained by cov(x, use='pairwise.complete.obs') in R:
covmat <- cov(sweep(simdata$genotypes,1,rowMeans(simdata$genotypes),"-"))
| Â | I1 | I2 | I3 | I4 | I5 |
|---|---|---|---|---|---|
| I1 | 0.2579 | 0.1336 | -0.1603 | -0.1103 | -0.1209 |
| I2 | 0.1336 | 0.3427 | -0.1027 | -0.1891 | -0.1845 |
| I3 | -0.1603 | -0.1027 | 0.3306 | -0.07394 | 0.006364 |
| I4 | -0.1103 | -0.1891 | -0.07394 | 0.2715 | 0.1018 |
| I5 | -0.1209 | -0.1845 | 0.006364 | 0.1018 | 0.1973 |
note: you can apply this to any similarity matrix; we'll assume we're using the covariance matrix.
Third: eigendecomposition
pca <- eigen(covmat)
- Eigenvalues:
- proportion of variance explained
pca$values/sum(pca$values)
## PC1 PC2 PC3 PC4 PC5 ## 5.231685e-01 2.938457e-01 1.033388e-01 7.964699e-02 4.038463e-17
- Eigenvectors:
- the positions of the samples on the principal components
pca$vectors
## PC1 PC2 PC3 PC4 PC5 ## I1 0.4637239 -0.2550913 0.6865902 -0.2201873 0.4472136 ## I2 0.6188658 0.1126111 -0.5991463 0.2129496 0.4472136 ## I3 -0.2742593 0.8049064 0.1115687 -0.2538897 0.4472136 ## I4 -0.4090639 -0.4969971 -0.3615817 -0.5048953 0.4472136 ## I5 -0.3992665 -0.1654291 0.1625691 0.7660227 0.4472136
Fourth: making pretty pictures
The eigenvectors are the coordinates of the samples on the principal components: so just plot them against each other:
pairs( pca$vectors[,1:3] )
Properties of PCA: math
"The principal components are the axes that explain the most variance."
Technically, this means that the top \(k\) PCs give the best \(k\)-dimensional approximation to the simliarity matrix: with eigenvalues \(\lambda_i\), eigenvectors \(v_i\), and similiarity matrix \(\Sigma\):
\(\Sigma \approx \sum_{i=1}^k \lambda_i v_i v_i^T\)
More is true: if \(\tilde G\) is the matrix of genotypes, with row and column means subtracted off, then there are vectors of allele weights \(w_i\) so that
\(\tilde G \approx \sum_{i=1}^k \sqrt{\lambda_i} w_i v_i^T .\)
Another consequence: the PCs can be expressed as weighted sums of the genotypes. (so: admixed individuals show up intermediate to parentals)
Technical issues
- Normalization:
you can choose to upweight low frequency alleles (e.g. divide by \(1/\sqrt{p(1-p)}\)); this may or may not make a difference.
- Linkage:
large linkage blocks (e.g. inversions) can hijack the top PCs.
Solution: prune linked SNPs.
- Sample sizes:
if the sampling scheme is unbalanced, the top PCs may choose to explain the variance within the most heavily sampled population, rather than variation between populations.
Solution: downweight by sample size:
weighted.covmat <- covmat * weights[col(covmat)] * weights[row(covmat)] eigen(weighted.covmat)
Doing PCA: demonstrated
Simulated data:
nind <- 1000 nloci <- 1e4 simdata <- sim_data(nind=nind,nloci=nloci) plot( simdata$locs, col=cols, pch=pchs, main="sampling scheme" )
PCA on simulated data:
norm.genotypes <- sweep( simdata$genotypes, 1, rowMeans(simdata$genotypes), "-" ) covmat <- cov(norm.genotypes) pca <- eigen(covmat) plot(pca$values/sum(pca$values),ylab="proportion", main="proportion of variance explained", xlab="PC#")
PCA map:
layout(t(1:2)) plot(simdata$locs, col=cols, pch=pchs, main="locations" ) plot(pca$vectors[,1:2], col=cols, pch=pchs, main="PC map", xlab="PC1", ylab="PC2" )
More PCs:
pairs(pca$vectors[,1:5], col=cols, pch=pchs )
PCs on the map:
The PCs here are Fourier modes; so the PC plots are like Lissajous figures (discussion in Novembre & Stephens).
Doing PCA: hands-on
Your turn: the POPRES
Same data as in Genes mirror geography within Europe, but without subsetting:
popres/crossprod_all-covariance.tsv: covariance matrix for all chromosomespopres/crossprod_chr8-covariance.tsv: covariance matrix for just chromosome 8popres/indivinfo.tsv: sample informationsource(popres/indivinfo.R): for country abbreviations and colors
Your tasks:
- Find the eigendecomposition.
- What proportion of variance is explained by each PC?
- Plot the first two PCs against each other.
- Do the same thing after subsampling to even out number of samples from each country.
- If there's time, look at chromosome 8.
Continuous geography
Not all populations look like this
Continuous space
Everything is related to everything else, but near things are more related than distant things.
Waldo Tobler, 1970
Why?
Offspring tend to live near to their parents, and so neighbors tend to
- be more closely related,
- i.e. have fewer generations separting them from their common ancestors,
- i.e. have more similar genomes.
"Isolation by distance"
correlation of genetic and geographic distances
From a population viewpoint, local drift differentiates populations, while migration makes nearby ones more similar.
From the individual viewpoint, lower population sizes means more recent, nearby shared ancestors; a tendency counteracted by migration.
Isolation by distance
Isolation by distance, environment
Isolation by *
To study IBE we must account for IBD
Cartoon model
- Sample some alleles at two locations, \(x\) and \(y\),
- which descend from some ancestors a long time ago.
- Each ancestor's genotype was chosen independently, with frequency \(f\);
- samples sharing the same ancestor have the same allele;
- so the modern frequency is near \(f\), but has drifted.
- Nearby samples are more likely to share a recent common ancestor,
- so the allele frequencies at \(x\) and \(y\) are correlated
- with a strength that depends on their distance apart.
Aside on correlated Gaussians
- Multivariate Gaussian:
A vector \((X_1,\ldots,X_n)\) is multivariate Gaussian with covariance matrix \((\Sigma_{ij})_{i,j=1}^n\) if any linear combination \(a_1 X_1 + \cdots a_n X_n\) is Gaussian with variance \(\sum_{ij} a_i \Sigma_{ij} a_j\).
- Gaussian process:
A random function \(X\) on a region is a Gaussian process with covariance function \(\sigma(s,t)\) if for each collection of locations \(t_1, \ldots, t_n\), the vector \((X(t_1), \ldots, X(t_n))\) is multivariate Gaussian with covariance matrix \(\Sigma_{ij} = \sigma(t_i,t_j)\).
This means: a continuous-space model of correlated allele frequencies reduces to a covariance matrix.
Simulated Gaussian field
covfn <- function (s,t) { exp(-sqrt(sum((s-t)^2))) }
st <- expand.grid( s=(0:30)/30, t=(0:30)/30 )
covmat <- diag(nrow(st))
for (i in 1:(nrow(st)-1)) for (j in (i+1):nrow(st)) {
covmat[i,j] <- covmat[j,i] <- covfn(st[i,],st[j,])
}
X <- MASS::mvrnorm(n=1,mu=rep(0,nrow(st)),Sigma=covmat)
dim(X) <- c(31,31)
Gestalt model
The allele frequency at locus \(\ell\) in population \(k\) is \[f_{\ell k} = \text{( global mean )} + \text{( spatially autocorrelated noise )}\]
and it would be nice if the noise was Gaussian.
The sampled alleles are independent\({}^*\) draws from these frequencies.
But: need to transform from \((-\infty,\infty)\) (Gaussians) to \([0,1]\) (frequencies).
\({}^*\) actually, Beta-Binomial, to allow for population-specific inbreeding coefficients.
Mathematical model
The allele frequency at locus \(\ell\) in population \(k\) is
\(f_{\ell k} = \frac{1}{1+e^{-(\Theta_{\ell k}+\mu_\ell)}}\)
where \(\mu_\ell\) is the (transformed) global mean allele frequency
and \(\Theta\) are multivariate Gaussian deviations
\(\Theta_{\ell k} \sim N(0,\Omega)\)
- with covariance matrix \(\Omega_{}\).
A parametric model:
With
\(D_{i,j} = \text{(geographic distance between samples *i* and *j*)}\)
and
\(E_{i,j} = \text{(environmental distance between samples *i* and *j*)}\)
we use the flexible parametric form
\(\Omega_{i,j} = \frac{1}{\alpha_0} \exp\left\{-\left(\sqrt{ \alpha_D D_{ij}^2 + \alpha_E E_{ij}^2 }\right)^{\alpha_2} \right\}.\)
we call this model
Bayesian
Estimation of
Differentiation in
Alleles by
Spatial
Structure and
Local
Ecology
Estimation of
Differentiation in
Alleles by
Spatial
Structure and
Local
Ecology
a.k.a.
Using BEDASSLE
Method mechanics: input
- Sampled populations (or individuals)
- Data at sampled locations:
- lat/long, eco/enviro data
- Genetic data (numerically coded, as above)
Recap of model, and analysis goal
Populations covary in their allele frequencies.
Extent of covariance depends on the distances between them, geographic and ecological.
We want to find the relative contribution of these distances.
The main focus is to infer \(\alpha_E/\alpha_D\).
We do this by putting priors on the parameters, and finding the posterior with MCMC.
- \(\alpha_D\) : contribution of geographic distance
- \(\alpha_{E,i}\) : contribution of eco/enviromental distance \(i\)
- \(\alpha_0\), \(\alpha_2\) : scaling parameters
- \(\Phi_k\) : population-specific "inbreeding coefficients"
- \(\mu_\ell\), \(\Theta_{k,\ell}\) : allele-specific transformed frequencies
Lightning introduction to MCMC
Markov chain Monte Carlo is a class of algorithms for parameter inference: as applied here,
given
- some data
- a generative model of that data (with a computable likelihood function)
- prior distributions on the parameters of that model
it finds
- samples from the posterior distribution on the parameters
- which lets you estimate the posterior distribution.
Reminder:
\(\mathbb{P}(\text{parameters}\,|\,\text{data}) \propto \mathbb{P}(\text{parameters}) \mathbb{P}(\text{data}\,|\,\text{parameters})\)
Our goal is to move around parameter space in such a way that the proportion of time spent near some parameter combination is proprtional to the posterior probability of those parameters.
How MCMC works
Two ingredients:
- proposal distribution ("where should I think about moving to?")
- acceptance probabilities ("hm, is moving there a good idea?")
These are such that iterating (proposal, accept-or-reject) gives the desired distribution. (e.g. the higher the posterior distribution at a proposed location, the more likely it is to move there)
MCMC checklist:
- acceptance rates: should be around 20-40% for optimal efficiency; increasing step size decreases acceptance rate
- likelihood profile: increases at first, then levels off: run longer if it's still increasing
- mixing and stationarity: need to run it long enough that it mixes, i.e. samples the entire relevent set of parameters: run longer if it hasn't "gone everywhere several times"
Tune on short runs; then do some longer ones.
Results: the posterior distribution, summarized by marginals.
Empirical example: teosinte
data: J Ross-Ibarra
Teosinte:
- a wild relative of domesticated maize
- ongoing introgression
- … responsible for high-altitude adaptation?
Is teosinte locally adapted to altitude?
Empirical example: teosinte
data: J Ross-Ibarra
Empirical example: teosinte
Trace plots from the MCMC: 
Conclusion: 100 m elevation \(\approx\) 15 km distance
Next steps: for which reason?
Simulated data
nind <- 50; nloci <- 1000 locs <- data.frame( lat=runif(nind), lon=runif(nind) ) env <- ifelse( locs$lat>0.5, 1, 0 ) sim <- sim_data(nind=nind,nloci=nloci,locs=locs,env=env,aE=2)
Simulated data: IBDE
covmat <- cov(sweep(sim$genotypes,1,rowMeans(sim$genotypes),"-")) D <- fields::rdist(sim$locs) E <- fields::rdist(sim$env) genetic.dist <- fields::rdist(t(sim$genotypes))/nloci
PCA preview
tplot <- function (xy,col=1,...) { plot(xy,type='n',...); text(xy,labels=labels,col=col) }
pca <- eigen(covmat)
layout(t(1:2))
tplot( locs, main="location", col=1+E )
tplot( pca$vectors[,1:2], xlab="PC1", ylab="PC2", main="PC map", col=1+E )
PCA preview
tplot( locs, col=cols[cut(pca$vectors[,1],16)], pch=20, main="PC1" ); abline(v=0.5) tplot( locs, col=cols[cut(pca$vectors[,2],16)], pch=20, main="PC2" ); abline(v=0.5)
Running BEDASSLE
Method mechanics: input
t(head(sim$genotypes)) # counts
## SNP1 SNP2 SNP3 SNP4 SNP5 SNP6 ## I1 1 0 2 2 0 2 ## I2 0 0 2 2 1 2 ## I3 0 0 2 2 0 2 ## I4 0 0 2 1 0 2 ## I5 0 0 2 2 0 2 ## I6 0 0 2 2 0 2 ## I7 0 0 2 2 0 2 ## I8 0 0 2 2 0 2 ## I9 0 0 2 2 1 2 ## I10 0 0 2 2 0 2 ## I11 0 0 2 2 0 2 ## I12 0 1 2 2 0 2 ## I13 0 1 2 2 0 2 ## I14 0 0 2 2 0 2 ## I15 0 0 2 2 1 2 ## I16 0 0 2 2 0 2 ## I17 0 0 2 2 0 2 ## I18 0 0 2 2 0 2 ## I19 0 0 2 2 0 2 ## I20 0 0 2 2 1 2 ## I21 0 0 2 2 0 2 ## I22 0 0 2 2 1 2 ## I23 0 0 2 2 0 2 ## I24 0 0 2 2 0 2 ## I25 0 0 2 2 0 2 ## I26 0 0 2 2 1 2 ## I27 0 0 2 2 0 2 ## I28 0 0 2 2 1 2 ## I29 0 0 2 2 1 2 ## I30 0 0 2 2 0 2 ## I31 0 0 2 2 0 2 ## I32 0 0 2 2 1 2 ## I33 0 0 2 2 2 2 ## I34 0 0 2 2 1 2 ## I35 0 0 2 2 0 2 ## I36 0 0 2 2 1 2 ## I37 0 0 2 2 1 2 ## I38 0 0 2 2 0 2 ## I39 0 0 2 2 2 2 ## I40 0 0 2 2 0 2 ## I41 0 0 2 2 0 2 ## I42 0 0 2 2 0 2 ## I43 0 0 2 2 1 2 ## I44 0 0 2 2 0 2 ## I45 0 0 2 2 0 2 ## I46 0 0 2 2 0 2 ## I47 0 0 2 2 0 2 ## I48 0 0 1 2 1 2 ## I49 0 0 2 2 1 2 ## I50 0 0 2 2 0 2
t(head(matrix(2,nrow=nind,ncol=nloci))) # sample sizes
## [,1] [,2] [,3] [,4] [,5] [,6] ## [1,] 2 2 2 2 2 2 ## [2,] 2 2 2 2 2 2 ## [3,] 2 2 2 2 2 2 ## [4,] 2 2 2 2 2 2 ## [5,] 2 2 2 2 2 2 ## [6,] 2 2 2 2 2 2 ## [7,] 2 2 2 2 2 2 ## [8,] 2 2 2 2 2 2 ## [9,] 2 2 2 2 2 2 ## [10,] 2 2 2 2 2 2 ## [11,] 2 2 2 2 2 2 ## [12,] 2 2 2 2 2 2 ## [13,] 2 2 2 2 2 2 ## [14,] 2 2 2 2 2 2 ## [15,] 2 2 2 2 2 2 ## [16,] 2 2 2 2 2 2 ## [17,] 2 2 2 2 2 2 ## [18,] 2 2 2 2 2 2 ## [19,] 2 2 2 2 2 2 ## [20,] 2 2 2 2 2 2 ## [21,] 2 2 2 2 2 2 ## [22,] 2 2 2 2 2 2 ## [23,] 2 2 2 2 2 2 ## [24,] 2 2 2 2 2 2 ## [25,] 2 2 2 2 2 2 ## [26,] 2 2 2 2 2 2 ## [27,] 2 2 2 2 2 2 ## [28,] 2 2 2 2 2 2 ## [29,] 2 2 2 2 2 2 ## [30,] 2 2 2 2 2 2 ## [31,] 2 2 2 2 2 2 ## [32,] 2 2 2 2 2 2 ## [33,] 2 2 2 2 2 2 ## [34,] 2 2 2 2 2 2 ## [35,] 2 2 2 2 2 2 ## [36,] 2 2 2 2 2 2 ## [37,] 2 2 2 2 2 2 ## [38,] 2 2 2 2 2 2 ## [39,] 2 2 2 2 2 2 ## [40,] 2 2 2 2 2 2 ## [41,] 2 2 2 2 2 2 ## [42,] 2 2 2 2 2 2 ## [43,] 2 2 2 2 2 2 ## [44,] 2 2 2 2 2 2 ## [45,] 2 2 2 2 2 2 ## [46,] 2 2 2 2 2 2 ## [47,] 2 2 2 2 2 2 ## [48,] 2 2 2 2 2 2 ## [49,] 2 2 2 2 2 2 ## [50,] 2 2 2 2 2 2 ## [51,] 2 2 2 2 2 2 ## [52,] 2 2 2 2 2 2 ## [53,] 2 2 2 2 2 2 ## [54,] 2 2 2 2 2 2 ## [55,] 2 2 2 2 2 2 ## [56,] 2 2 2 2 2 2 ## [57,] 2 2 2 2 2 2 ## [58,] 2 2 2 2 2 2 ## [59,] 2 2 2 2 2 2 ## [60,] 2 2 2 2 2 2 ## [61,] 2 2 2 2 2 2 ## [62,] 2 2 2 2 2 2 ## [63,] 2 2 2 2 2 2 ## [64,] 2 2 2 2 2 2 ## [65,] 2 2 2 2 2 2 ## [66,] 2 2 2 2 2 2 ## [67,] 2 2 2 2 2 2 ## [68,] 2 2 2 2 2 2 ## [69,] 2 2 2 2 2 2 ## [70,] 2 2 2 2 2 2 ## [71,] 2 2 2 2 2 2 ## [72,] 2 2 2 2 2 2 ## [73,] 2 2 2 2 2 2 ## [74,] 2 2 2 2 2 2 ## [75,] 2 2 2 2 2 2 ## [76,] 2 2 2 2 2 2 ## [77,] 2 2 2 2 2 2 ## [78,] 2 2 2 2 2 2 ## [79,] 2 2 2 2 2 2 ## [80,] 2 2 2 2 2 2 ## [81,] 2 2 2 2 2 2 ## [82,] 2 2 2 2 2 2 ## [83,] 2 2 2 2 2 2 ## [84,] 2 2 2 2 2 2 ## [85,] 2 2 2 2 2 2 ## [86,] 2 2 2 2 2 2 ## [87,] 2 2 2 2 2 2 ## [88,] 2 2 2 2 2 2 ## [89,] 2 2 2 2 2 2 ## [90,] 2 2 2 2 2 2 ## [91,] 2 2 2 2 2 2 ## [92,] 2 2 2 2 2 2 ## [93,] 2 2 2 2 2 2 ## [94,] 2 2 2 2 2 2 ## [95,] 2 2 2 2 2 2 ## [96,] 2 2 2 2 2 2 ## [97,] 2 2 2 2 2 2 ## [98,] 2 2 2 2 2 2 ## [99,] 2 2 2 2 2 2 ## [100,] 2 2 2 2 2 2 ## [101,] 2 2 2 2 2 2 ## [102,] 2 2 2 2 2 2 ## [103,] 2 2 2 2 2 2 ## [104,] 2 2 2 2 2 2 ## [105,] 2 2 2 2 2 2 ## [106,] 2 2 2 2 2 2 ## [107,] 2 2 2 2 2 2 ## [108,] 2 2 2 2 2 2 ## [109,] 2 2 2 2 2 2 ## [110,] 2 2 2 2 2 2 ## [111,] 2 2 2 2 2 2 ## [112,] 2 2 2 2 2 2 ## [113,] 2 2 2 2 2 2 ## [114,] 2 2 2 2 2 2 ## [115,] 2 2 2 2 2 2 ## [116,] 2 2 2 2 2 2 ## [117,] 2 2 2 2 2 2 ## [118,] 2 2 2 2 2 2 ## [119,] 2 2 2 2 2 2 ## [120,] 2 2 2 2 2 2 ## [121,] 2 2 2 2 2 2 ## [122,] 2 2 2 2 2 2 ## [123,] 2 2 2 2 2 2 ## [124,] 2 2 2 2 2 2 ## [125,] 2 2 2 2 2 2 ## [126,] 2 2 2 2 2 2 ## [127,] 2 2 2 2 2 2 ## [128,] 2 2 2 2 2 2 ## [129,] 2 2 2 2 2 2 ## [130,] 2 2 2 2 2 2 ## [131,] 2 2 2 2 2 2 ## [132,] 2 2 2 2 2 2 ## [133,] 2 2 2 2 2 2 ## [134,] 2 2 2 2 2 2 ## [135,] 2 2 2 2 2 2 ## [136,] 2 2 2 2 2 2 ## [137,] 2 2 2 2 2 2 ## [138,] 2 2 2 2 2 2 ## [139,] 2 2 2 2 2 2 ## [140,] 2 2 2 2 2 2 ## [141,] 2 2 2 2 2 2 ## [142,] 2 2 2 2 2 2 ## [143,] 2 2 2 2 2 2 ## [144,] 2 2 2 2 2 2 ## [145,] 2 2 2 2 2 2 ## [146,] 2 2 2 2 2 2 ## [147,] 2 2 2 2 2 2 ## [148,] 2 2 2 2 2 2 ## [149,] 2 2 2 2 2 2 ## [150,] 2 2 2 2 2 2 ## [151,] 2 2 2 2 2 2 ## [152,] 2 2 2 2 2 2 ## [153,] 2 2 2 2 2 2 ## [154,] 2 2 2 2 2 2 ## [155,] 2 2 2 2 2 2 ## [156,] 2 2 2 2 2 2 ## [157,] 2 2 2 2 2 2 ## [158,] 2 2 2 2 2 2 ## [159,] 2 2 2 2 2 2 ## [160,] 2 2 2 2 2 2 ## [161,] 2 2 2 2 2 2 ## [162,] 2 2 2 2 2 2 ## [163,] 2 2 2 2 2 2 ## [164,] 2 2 2 2 2 2 ## [165,] 2 2 2 2 2 2 ## [166,] 2 2 2 2 2 2 ## [167,] 2 2 2 2 2 2 ## [168,] 2 2 2 2 2 2 ## [169,] 2 2 2 2 2 2 ## [170,] 2 2 2 2 2 2 ## [171,] 2 2 2 2 2 2 ## [172,] 2 2 2 2 2 2 ## [173,] 2 2 2 2 2 2 ## [174,] 2 2 2 2 2 2 ## [175,] 2 2 2 2 2 2 ## [176,] 2 2 2 2 2 2 ## [177,] 2 2 2 2 2 2 ## [178,] 2 2 2 2 2 2 ## [179,] 2 2 2 2 2 2 ## [180,] 2 2 2 2 2 2 ## [181,] 2 2 2 2 2 2 ## [182,] 2 2 2 2 2 2 ## [183,] 2 2 2 2 2 2 ## [184,] 2 2 2 2 2 2 ## [185,] 2 2 2 2 2 2 ## [186,] 2 2 2 2 2 2 ## [187,] 2 2 2 2 2 2 ## [188,] 2 2 2 2 2 2 ## [189,] 2 2 2 2 2 2 ## [190,] 2 2 2 2 2 2 ## [191,] 2 2 2 2 2 2 ## [192,] 2 2 2 2 2 2 ## [193,] 2 2 2 2 2 2 ## [194,] 2 2 2 2 2 2 ## [195,] 2 2 2 2 2 2 ## [196,] 2 2 2 2 2 2 ## [197,] 2 2 2 2 2 2 ## [198,] 2 2 2 2 2 2 ## [199,] 2 2 2 2 2 2 ## [200,] 2 2 2 2 2 2 ## [201,] 2 2 2 2 2 2 ## [202,] 2 2 2 2 2 2 ## [203,] 2 2 2 2 2 2 ## [204,] 2 2 2 2 2 2 ## [205,] 2 2 2 2 2 2 ## [206,] 2 2 2 2 2 2 ## [207,] 2 2 2 2 2 2 ## [208,] 2 2 2 2 2 2 ## [209,] 2 2 2 2 2 2 ## [210,] 2 2 2 2 2 2 ## [211,] 2 2 2 2 2 2 ## [212,] 2 2 2 2 2 2 ## [213,] 2 2 2 2 2 2 ## [214,] 2 2 2 2 2 2 ## [215,] 2 2 2 2 2 2 ## [216,] 2 2 2 2 2 2 ## [217,] 2 2 2 2 2 2 ## [218,] 2 2 2 2 2 2 ## [219,] 2 2 2 2 2 2 ## [220,] 2 2 2 2 2 2 ## [221,] 2 2 2 2 2 2 ## [222,] 2 2 2 2 2 2 ## [223,] 2 2 2 2 2 2 ## [224,] 2 2 2 2 2 2 ## [225,] 2 2 2 2 2 2 ## [226,] 2 2 2 2 2 2 ## [227,] 2 2 2 2 2 2 ## [228,] 2 2 2 2 2 2 ## [229,] 2 2 2 2 2 2 ## [230,] 2 2 2 2 2 2 ## [231,] 2 2 2 2 2 2 ## [232,] 2 2 2 2 2 2 ## [233,] 2 2 2 2 2 2 ## [234,] 2 2 2 2 2 2 ## [235,] 2 2 2 2 2 2 ## [236,] 2 2 2 2 2 2 ## [237,] 2 2 2 2 2 2 ## [238,] 2 2 2 2 2 2 ## [239,] 2 2 2 2 2 2 ## [240,] 2 2 2 2 2 2 ## [241,] 2 2 2 2 2 2 ## [242,] 2 2 2 2 2 2 ## [243,] 2 2 2 2 2 2 ## [244,] 2 2 2 2 2 2 ## [245,] 2 2 2 2 2 2 ## [246,] 2 2 2 2 2 2 ## [247,] 2 2 2 2 2 2 ## [248,] 2 2 2 2 2 2 ## [249,] 2 2 2 2 2 2 ## [250,] 2 2 2 2 2 2 ## [251,] 2 2 2 2 2 2 ## [252,] 2 2 2 2 2 2 ## [253,] 2 2 2 2 2 2 ## [254,] 2 2 2 2 2 2 ## [255,] 2 2 2 2 2 2 ## [256,] 2 2 2 2 2 2 ## [257,] 2 2 2 2 2 2 ## [258,] 2 2 2 2 2 2 ## [259,] 2 2 2 2 2 2 ## [260,] 2 2 2 2 2 2 ## [261,] 2 2 2 2 2 2 ## [262,] 2 2 2 2 2 2 ## [263,] 2 2 2 2 2 2 ## [264,] 2 2 2 2 2 2 ## [265,] 2 2 2 2 2 2 ## [266,] 2 2 2 2 2 2 ## [267,] 2 2 2 2 2 2 ## [268,] 2 2 2 2 2 2 ## [269,] 2 2 2 2 2 2 ## [270,] 2 2 2 2 2 2 ## [271,] 2 2 2 2 2 2 ## [272,] 2 2 2 2 2 2 ## [273,] 2 2 2 2 2 2 ## [274,] 2 2 2 2 2 2 ## [275,] 2 2 2 2 2 2 ## [276,] 2 2 2 2 2 2 ## [277,] 2 2 2 2 2 2 ## [278,] 2 2 2 2 2 2 ## [279,] 2 2 2 2 2 2 ## [280,] 2 2 2 2 2 2 ## [281,] 2 2 2 2 2 2 ## [282,] 2 2 2 2 2 2 ## [283,] 2 2 2 2 2 2 ## [284,] 2 2 2 2 2 2 ## [285,] 2 2 2 2 2 2 ## [286,] 2 2 2 2 2 2 ## [287,] 2 2 2 2 2 2 ## [288,] 2 2 2 2 2 2 ## [289,] 2 2 2 2 2 2 ## [290,] 2 2 2 2 2 2 ## [291,] 2 2 2 2 2 2 ## [292,] 2 2 2 2 2 2 ## [293,] 2 2 2 2 2 2 ## [294,] 2 2 2 2 2 2 ## [295,] 2 2 2 2 2 2 ## [296,] 2 2 2 2 2 2 ## [297,] 2 2 2 2 2 2 ## [298,] 2 2 2 2 2 2 ## [299,] 2 2 2 2 2 2 ## [300,] 2 2 2 2 2 2 ## [301,] 2 2 2 2 2 2 ## [302,] 2 2 2 2 2 2 ## [303,] 2 2 2 2 2 2 ## [304,] 2 2 2 2 2 2 ## [305,] 2 2 2 2 2 2 ## [306,] 2 2 2 2 2 2 ## [307,] 2 2 2 2 2 2 ## [308,] 2 2 2 2 2 2 ## [309,] 2 2 2 2 2 2 ## [310,] 2 2 2 2 2 2 ## [311,] 2 2 2 2 2 2 ## [312,] 2 2 2 2 2 2 ## [313,] 2 2 2 2 2 2 ## [314,] 2 2 2 2 2 2 ## [315,] 2 2 2 2 2 2 ## [316,] 2 2 2 2 2 2 ## [317,] 2 2 2 2 2 2 ## [318,] 2 2 2 2 2 2 ## [319,] 2 2 2 2 2 2 ## [320,] 2 2 2 2 2 2 ## [321,] 2 2 2 2 2 2 ## [322,] 2 2 2 2 2 2 ## [323,] 2 2 2 2 2 2 ## [324,] 2 2 2 2 2 2 ## [325,] 2 2 2 2 2 2 ## [326,] 2 2 2 2 2 2 ## [327,] 2 2 2 2 2 2 ## [328,] 2 2 2 2 2 2 ## [329,] 2 2 2 2 2 2 ## [330,] 2 2 2 2 2 2 ## [331,] 2 2 2 2 2 2 ## [332,] 2 2 2 2 2 2 ## [333,] 2 2 2 2 2 2 ## [334,] 2 2 2 2 2 2 ## [335,] 2 2 2 2 2 2 ## [336,] 2 2 2 2 2 2 ## [337,] 2 2 2 2 2 2 ## [338,] 2 2 2 2 2 2 ## [339,] 2 2 2 2 2 2 ## [340,] 2 2 2 2 2 2 ## [341,] 2 2 2 2 2 2 ## [342,] 2 2 2 2 2 2 ## [343,] 2 2 2 2 2 2 ## [344,] 2 2 2 2 2 2 ## [345,] 2 2 2 2 2 2 ## [346,] 2 2 2 2 2 2 ## [347,] 2 2 2 2 2 2 ## [348,] 2 2 2 2 2 2 ## [349,] 2 2 2 2 2 2 ## [350,] 2 2 2 2 2 2 ## [351,] 2 2 2 2 2 2 ## [352,] 2 2 2 2 2 2 ## [353,] 2 2 2 2 2 2 ## [354,] 2 2 2 2 2 2 ## [355,] 2 2 2 2 2 2 ## [356,] 2 2 2 2 2 2 ## [357,] 2 2 2 2 2 2 ## [358,] 2 2 2 2 2 2 ## [359,] 2 2 2 2 2 2 ## [360,] 2 2 2 2 2 2 ## [361,] 2 2 2 2 2 2 ## [362,] 2 2 2 2 2 2 ## [363,] 2 2 2 2 2 2 ## [364,] 2 2 2 2 2 2 ## [365,] 2 2 2 2 2 2 ## [366,] 2 2 2 2 2 2 ## [367,] 2 2 2 2 2 2 ## [368,] 2 2 2 2 2 2 ## [369,] 2 2 2 2 2 2 ## [370,] 2 2 2 2 2 2 ## [371,] 2 2 2 2 2 2 ## [372,] 2 2 2 2 2 2 ## [373,] 2 2 2 2 2 2 ## [374,] 2 2 2 2 2 2 ## [375,] 2 2 2 2 2 2 ## [376,] 2 2 2 2 2 2 ## [377,] 2 2 2 2 2 2 ## [378,] 2 2 2 2 2 2 ## [379,] 2 2 2 2 2 2 ## [380,] 2 2 2 2 2 2 ## [381,] 2 2 2 2 2 2 ## [382,] 2 2 2 2 2 2 ## [383,] 2 2 2 2 2 2 ## [384,] 2 2 2 2 2 2 ## [385,] 2 2 2 2 2 2 ## [386,] 2 2 2 2 2 2 ## [387,] 2 2 2 2 2 2 ## [388,] 2 2 2 2 2 2 ## [389,] 2 2 2 2 2 2 ## [390,] 2 2 2 2 2 2 ## [391,] 2 2 2 2 2 2 ## [392,] 2 2 2 2 2 2 ## [393,] 2 2 2 2 2 2 ## [394,] 2 2 2 2 2 2 ## [395,] 2 2 2 2 2 2 ## [396,] 2 2 2 2 2 2 ## [397,] 2 2 2 2 2 2 ## [398,] 2 2 2 2 2 2 ## [399,] 2 2 2 2 2 2 ## [400,] 2 2 2 2 2 2 ## [401,] 2 2 2 2 2 2 ## [402,] 2 2 2 2 2 2 ## [403,] 2 2 2 2 2 2 ## [404,] 2 2 2 2 2 2 ## [405,] 2 2 2 2 2 2 ## [406,] 2 2 2 2 2 2 ## [407,] 2 2 2 2 2 2 ## [408,] 2 2 2 2 2 2 ## [409,] 2 2 2 2 2 2 ## [410,] 2 2 2 2 2 2 ## [411,] 2 2 2 2 2 2 ## [412,] 2 2 2 2 2 2 ## [413,] 2 2 2 2 2 2 ## [414,] 2 2 2 2 2 2 ## [415,] 2 2 2 2 2 2 ## [416,] 2 2 2 2 2 2 ## [417,] 2 2 2 2 2 2 ## [418,] 2 2 2 2 2 2 ## [419,] 2 2 2 2 2 2 ## [420,] 2 2 2 2 2 2 ## [421,] 2 2 2 2 2 2 ## [422,] 2 2 2 2 2 2 ## [423,] 2 2 2 2 2 2 ## [424,] 2 2 2 2 2 2 ## [425,] 2 2 2 2 2 2 ## [426,] 2 2 2 2 2 2 ## [427,] 2 2 2 2 2 2 ## [428,] 2 2 2 2 2 2 ## [429,] 2 2 2 2 2 2 ## [430,] 2 2 2 2 2 2 ## [431,] 2 2 2 2 2 2 ## [432,] 2 2 2 2 2 2 ## [433,] 2 2 2 2 2 2 ## [434,] 2 2 2 2 2 2 ## [435,] 2 2 2 2 2 2 ## [436,] 2 2 2 2 2 2 ## [437,] 2 2 2 2 2 2 ## [438,] 2 2 2 2 2 2 ## [439,] 2 2 2 2 2 2 ## [440,] 2 2 2 2 2 2 ## [441,] 2 2 2 2 2 2 ## [442,] 2 2 2 2 2 2 ## [443,] 2 2 2 2 2 2 ## [444,] 2 2 2 2 2 2 ## [445,] 2 2 2 2 2 2 ## [446,] 2 2 2 2 2 2 ## [447,] 2 2 2 2 2 2 ## [448,] 2 2 2 2 2 2 ## [449,] 2 2 2 2 2 2 ## [450,] 2 2 2 2 2 2 ## [451,] 2 2 2 2 2 2 ## [452,] 2 2 2 2 2 2 ## [453,] 2 2 2 2 2 2 ## [454,] 2 2 2 2 2 2 ## [455,] 2 2 2 2 2 2 ## [456,] 2 2 2 2 2 2 ## [457,] 2 2 2 2 2 2 ## [458,] 2 2 2 2 2 2 ## [459,] 2 2 2 2 2 2 ## [460,] 2 2 2 2 2 2 ## [461,] 2 2 2 2 2 2 ## [462,] 2 2 2 2 2 2 ## [463,] 2 2 2 2 2 2 ## [464,] 2 2 2 2 2 2 ## [465,] 2 2 2 2 2 2 ## [466,] 2 2 2 2 2 2 ## [467,] 2 2 2 2 2 2 ## [468,] 2 2 2 2 2 2 ## [469,] 2 2 2 2 2 2 ## [470,] 2 2 2 2 2 2 ## [471,] 2 2 2 2 2 2 ## [472,] 2 2 2 2 2 2 ## [473,] 2 2 2 2 2 2 ## [474,] 2 2 2 2 2 2 ## [475,] 2 2 2 2 2 2 ## [476,] 2 2 2 2 2 2 ## [477,] 2 2 2 2 2 2 ## [478,] 2 2 2 2 2 2 ## [479,] 2 2 2 2 2 2 ## [480,] 2 2 2 2 2 2 ## [481,] 2 2 2 2 2 2 ## [482,] 2 2 2 2 2 2 ## [483,] 2 2 2 2 2 2 ## [484,] 2 2 2 2 2 2 ## [485,] 2 2 2 2 2 2 ## [486,] 2 2 2 2 2 2 ## [487,] 2 2 2 2 2 2 ## [488,] 2 2 2 2 2 2 ## [489,] 2 2 2 2 2 2 ## [490,] 2 2 2 2 2 2 ## [491,] 2 2 2 2 2 2 ## [492,] 2 2 2 2 2 2 ## [493,] 2 2 2 2 2 2 ## [494,] 2 2 2 2 2 2 ## [495,] 2 2 2 2 2 2 ## [496,] 2 2 2 2 2 2 ## [497,] 2 2 2 2 2 2 ## [498,] 2 2 2 2 2 2 ## [499,] 2 2 2 2 2 2 ## [500,] 2 2 2 2 2 2 ## [501,] 2 2 2 2 2 2 ## [502,] 2 2 2 2 2 2 ## [503,] 2 2 2 2 2 2 ## [504,] 2 2 2 2 2 2 ## [505,] 2 2 2 2 2 2 ## [506,] 2 2 2 2 2 2 ## [507,] 2 2 2 2 2 2 ## [508,] 2 2 2 2 2 2 ## [509,] 2 2 2 2 2 2 ## [510,] 2 2 2 2 2 2 ## [511,] 2 2 2 2 2 2 ## [512,] 2 2 2 2 2 2 ## [513,] 2 2 2 2 2 2 ## [514,] 2 2 2 2 2 2 ## [515,] 2 2 2 2 2 2 ## [516,] 2 2 2 2 2 2 ## [517,] 2 2 2 2 2 2 ## [518,] 2 2 2 2 2 2 ## [519,] 2 2 2 2 2 2 ## [520,] 2 2 2 2 2 2 ## [521,] 2 2 2 2 2 2 ## [522,] 2 2 2 2 2 2 ## [523,] 2 2 2 2 2 2 ## [524,] 2 2 2 2 2 2 ## [525,] 2 2 2 2 2 2 ## [526,] 2 2 2 2 2 2 ## [527,] 2 2 2 2 2 2 ## [528,] 2 2 2 2 2 2 ## [529,] 2 2 2 2 2 2 ## [530,] 2 2 2 2 2 2 ## [531,] 2 2 2 2 2 2 ## [532,] 2 2 2 2 2 2 ## [533,] 2 2 2 2 2 2 ## [534,] 2 2 2 2 2 2 ## [535,] 2 2 2 2 2 2 ## [536,] 2 2 2 2 2 2 ## [537,] 2 2 2 2 2 2 ## [538,] 2 2 2 2 2 2 ## [539,] 2 2 2 2 2 2 ## [540,] 2 2 2 2 2 2 ## [541,] 2 2 2 2 2 2 ## [542,] 2 2 2 2 2 2 ## [543,] 2 2 2 2 2 2 ## [544,] 2 2 2 2 2 2 ## [545,] 2 2 2 2 2 2 ## [546,] 2 2 2 2 2 2 ## [547,] 2 2 2 2 2 2 ## [548,] 2 2 2 2 2 2 ## [549,] 2 2 2 2 2 2 ## [550,] 2 2 2 2 2 2 ## [551,] 2 2 2 2 2 2 ## [552,] 2 2 2 2 2 2 ## [553,] 2 2 2 2 2 2 ## [554,] 2 2 2 2 2 2 ## [555,] 2 2 2 2 2 2 ## [556,] 2 2 2 2 2 2 ## [557,] 2 2 2 2 2 2 ## [558,] 2 2 2 2 2 2 ## [559,] 2 2 2 2 2 2 ## [560,] 2 2 2 2 2 2 ## [561,] 2 2 2 2 2 2 ## [562,] 2 2 2 2 2 2 ## [563,] 2 2 2 2 2 2 ## [564,] 2 2 2 2 2 2 ## [565,] 2 2 2 2 2 2 ## [566,] 2 2 2 2 2 2 ## [567,] 2 2 2 2 2 2 ## [568,] 2 2 2 2 2 2 ## [569,] 2 2 2 2 2 2 ## [570,] 2 2 2 2 2 2 ## [571,] 2 2 2 2 2 2 ## [572,] 2 2 2 2 2 2 ## [573,] 2 2 2 2 2 2 ## [574,] 2 2 2 2 2 2 ## [575,] 2 2 2 2 2 2 ## [576,] 2 2 2 2 2 2 ## [577,] 2 2 2 2 2 2 ## [578,] 2 2 2 2 2 2 ## [579,] 2 2 2 2 2 2 ## [580,] 2 2 2 2 2 2 ## [581,] 2 2 2 2 2 2 ## [582,] 2 2 2 2 2 2 ## [583,] 2 2 2 2 2 2 ## [584,] 2 2 2 2 2 2 ## [585,] 2 2 2 2 2 2 ## [586,] 2 2 2 2 2 2 ## [587,] 2 2 2 2 2 2 ## [588,] 2 2 2 2 2 2 ## [589,] 2 2 2 2 2 2 ## [590,] 2 2 2 2 2 2 ## [591,] 2 2 2 2 2 2 ## [592,] 2 2 2 2 2 2 ## [593,] 2 2 2 2 2 2 ## [594,] 2 2 2 2 2 2 ## [595,] 2 2 2 2 2 2 ## [596,] 2 2 2 2 2 2 ## [597,] 2 2 2 2 2 2 ## [598,] 2 2 2 2 2 2 ## [599,] 2 2 2 2 2 2 ## [600,] 2 2 2 2 2 2 ## [601,] 2 2 2 2 2 2 ## [602,] 2 2 2 2 2 2 ## [603,] 2 2 2 2 2 2 ## [604,] 2 2 2 2 2 2 ## [605,] 2 2 2 2 2 2 ## [606,] 2 2 2 2 2 2 ## [607,] 2 2 2 2 2 2 ## [608,] 2 2 2 2 2 2 ## [609,] 2 2 2 2 2 2 ## [610,] 2 2 2 2 2 2 ## [611,] 2 2 2 2 2 2 ## [612,] 2 2 2 2 2 2 ## [613,] 2 2 2 2 2 2 ## [614,] 2 2 2 2 2 2 ## [615,] 2 2 2 2 2 2 ## [616,] 2 2 2 2 2 2 ## [617,] 2 2 2 2 2 2 ## [618,] 2 2 2 2 2 2 ## [619,] 2 2 2 2 2 2 ## [620,] 2 2 2 2 2 2 ## [621,] 2 2 2 2 2 2 ## [622,] 2 2 2 2 2 2 ## [623,] 2 2 2 2 2 2 ## [624,] 2 2 2 2 2 2 ## [625,] 2 2 2 2 2 2 ## [626,] 2 2 2 2 2 2 ## [627,] 2 2 2 2 2 2 ## [628,] 2 2 2 2 2 2 ## [629,] 2 2 2 2 2 2 ## [630,] 2 2 2 2 2 2 ## [631,] 2 2 2 2 2 2 ## [632,] 2 2 2 2 2 2 ## [633,] 2 2 2 2 2 2 ## [634,] 2 2 2 2 2 2 ## [635,] 2 2 2 2 2 2 ## [636,] 2 2 2 2 2 2 ## [637,] 2 2 2 2 2 2 ## [638,] 2 2 2 2 2 2 ## [639,] 2 2 2 2 2 2 ## [640,] 2 2 2 2 2 2 ## [641,] 2 2 2 2 2 2 ## [642,] 2 2 2 2 2 2 ## [643,] 2 2 2 2 2 2 ## [644,] 2 2 2 2 2 2 ## [645,] 2 2 2 2 2 2 ## [646,] 2 2 2 2 2 2 ## [647,] 2 2 2 2 2 2 ## [648,] 2 2 2 2 2 2 ## [649,] 2 2 2 2 2 2 ## [650,] 2 2 2 2 2 2 ## [651,] 2 2 2 2 2 2 ## [652,] 2 2 2 2 2 2 ## [653,] 2 2 2 2 2 2 ## [654,] 2 2 2 2 2 2 ## [655,] 2 2 2 2 2 2 ## [656,] 2 2 2 2 2 2 ## [657,] 2 2 2 2 2 2 ## [658,] 2 2 2 2 2 2 ## [659,] 2 2 2 2 2 2 ## [660,] 2 2 2 2 2 2 ## [661,] 2 2 2 2 2 2 ## [662,] 2 2 2 2 2 2 ## [663,] 2 2 2 2 2 2 ## [664,] 2 2 2 2 2 2 ## [665,] 2 2 2 2 2 2 ## [666,] 2 2 2 2 2 2 ## [667,] 2 2 2 2 2 2 ## [668,] 2 2 2 2 2 2 ## [669,] 2 2 2 2 2 2 ## [670,] 2 2 2 2 2 2 ## [671,] 2 2 2 2 2 2 ## [672,] 2 2 2 2 2 2 ## [673,] 2 2 2 2 2 2 ## [674,] 2 2 2 2 2 2 ## [675,] 2 2 2 2 2 2 ## [676,] 2 2 2 2 2 2 ## [677,] 2 2 2 2 2 2 ## [678,] 2 2 2 2 2 2 ## [679,] 2 2 2 2 2 2 ## [680,] 2 2 2 2 2 2 ## [681,] 2 2 2 2 2 2 ## [682,] 2 2 2 2 2 2 ## [683,] 2 2 2 2 2 2 ## [684,] 2 2 2 2 2 2 ## [685,] 2 2 2 2 2 2 ## [686,] 2 2 2 2 2 2 ## [687,] 2 2 2 2 2 2 ## [688,] 2 2 2 2 2 2 ## [689,] 2 2 2 2 2 2 ## [690,] 2 2 2 2 2 2 ## [691,] 2 2 2 2 2 2 ## [692,] 2 2 2 2 2 2 ## [693,] 2 2 2 2 2 2 ## [694,] 2 2 2 2 2 2 ## [695,] 2 2 2 2 2 2 ## [696,] 2 2 2 2 2 2 ## [697,] 2 2 2 2 2 2 ## [698,] 2 2 2 2 2 2 ## [699,] 2 2 2 2 2 2 ## [700,] 2 2 2 2 2 2 ## [701,] 2 2 2 2 2 2 ## [702,] 2 2 2 2 2 2 ## [703,] 2 2 2 2 2 2 ## [704,] 2 2 2 2 2 2 ## [705,] 2 2 2 2 2 2 ## [706,] 2 2 2 2 2 2 ## [707,] 2 2 2 2 2 2 ## [708,] 2 2 2 2 2 2 ## [709,] 2 2 2 2 2 2 ## [710,] 2 2 2 2 2 2 ## [711,] 2 2 2 2 2 2 ## [712,] 2 2 2 2 2 2 ## [713,] 2 2 2 2 2 2 ## [714,] 2 2 2 2 2 2 ## [715,] 2 2 2 2 2 2 ## [716,] 2 2 2 2 2 2 ## [717,] 2 2 2 2 2 2 ## [718,] 2 2 2 2 2 2 ## [719,] 2 2 2 2 2 2 ## [720,] 2 2 2 2 2 2 ## [721,] 2 2 2 2 2 2 ## [722,] 2 2 2 2 2 2 ## [723,] 2 2 2 2 2 2 ## [724,] 2 2 2 2 2 2 ## [725,] 2 2 2 2 2 2 ## [726,] 2 2 2 2 2 2 ## [727,] 2 2 2 2 2 2 ## [728,] 2 2 2 2 2 2 ## [729,] 2 2 2 2 2 2 ## [730,] 2 2 2 2 2 2 ## [731,] 2 2 2 2 2 2 ## [732,] 2 2 2 2 2 2 ## [733,] 2 2 2 2 2 2 ## [734,] 2 2 2 2 2 2 ## [735,] 2 2 2 2 2 2 ## [736,] 2 2 2 2 2 2 ## [737,] 2 2 2 2 2 2 ## [738,] 2 2 2 2 2 2 ## [739,] 2 2 2 2 2 2 ## [740,] 2 2 2 2 2 2 ## [741,] 2 2 2 2 2 2 ## [742,] 2 2 2 2 2 2 ## [743,] 2 2 2 2 2 2 ## [744,] 2 2 2 2 2 2 ## [745,] 2 2 2 2 2 2 ## [746,] 2 2 2 2 2 2 ## [747,] 2 2 2 2 2 2 ## [748,] 2 2 2 2 2 2 ## [749,] 2 2 2 2 2 2 ## [750,] 2 2 2 2 2 2 ## [751,] 2 2 2 2 2 2 ## [752,] 2 2 2 2 2 2 ## [753,] 2 2 2 2 2 2 ## [754,] 2 2 2 2 2 2 ## [755,] 2 2 2 2 2 2 ## [756,] 2 2 2 2 2 2 ## [757,] 2 2 2 2 2 2 ## [758,] 2 2 2 2 2 2 ## [759,] 2 2 2 2 2 2 ## [760,] 2 2 2 2 2 2 ## [761,] 2 2 2 2 2 2 ## [762,] 2 2 2 2 2 2 ## [763,] 2 2 2 2 2 2 ## [764,] 2 2 2 2 2 2 ## [765,] 2 2 2 2 2 2 ## [766,] 2 2 2 2 2 2 ## [767,] 2 2 2 2 2 2 ## [768,] 2 2 2 2 2 2 ## [769,] 2 2 2 2 2 2 ## [770,] 2 2 2 2 2 2 ## [771,] 2 2 2 2 2 2 ## [772,] 2 2 2 2 2 2 ## [773,] 2 2 2 2 2 2 ## [774,] 2 2 2 2 2 2 ## [775,] 2 2 2 2 2 2 ## [776,] 2 2 2 2 2 2 ## [777,] 2 2 2 2 2 2 ## [778,] 2 2 2 2 2 2 ## [779,] 2 2 2 2 2 2 ## [780,] 2 2 2 2 2 2 ## [781,] 2 2 2 2 2 2 ## [782,] 2 2 2 2 2 2 ## [783,] 2 2 2 2 2 2 ## [784,] 2 2 2 2 2 2 ## [785,] 2 2 2 2 2 2 ## [786,] 2 2 2 2 2 2 ## [787,] 2 2 2 2 2 2 ## [788,] 2 2 2 2 2 2 ## [789,] 2 2 2 2 2 2 ## [790,] 2 2 2 2 2 2 ## [791,] 2 2 2 2 2 2 ## [792,] 2 2 2 2 2 2 ## [793,] 2 2 2 2 2 2 ## [794,] 2 2 2 2 2 2 ## [795,] 2 2 2 2 2 2 ## [796,] 2 2 2 2 2 2 ## [797,] 2 2 2 2 2 2 ## [798,] 2 2 2 2 2 2 ## [799,] 2 2 2 2 2 2 ## [800,] 2 2 2 2 2 2 ## [801,] 2 2 2 2 2 2 ## [802,] 2 2 2 2 2 2 ## [803,] 2 2 2 2 2 2 ## [804,] 2 2 2 2 2 2 ## [805,] 2 2 2 2 2 2 ## [806,] 2 2 2 2 2 2 ## [807,] 2 2 2 2 2 2 ## [808,] 2 2 2 2 2 2 ## [809,] 2 2 2 2 2 2 ## [810,] 2 2 2 2 2 2 ## [811,] 2 2 2 2 2 2 ## [812,] 2 2 2 2 2 2 ## [813,] 2 2 2 2 2 2 ## [814,] 2 2 2 2 2 2 ## [815,] 2 2 2 2 2 2 ## [816,] 2 2 2 2 2 2 ## [817,] 2 2 2 2 2 2 ## [818,] 2 2 2 2 2 2 ## [819,] 2 2 2 2 2 2 ## [820,] 2 2 2 2 2 2 ## [821,] 2 2 2 2 2 2 ## [822,] 2 2 2 2 2 2 ## [823,] 2 2 2 2 2 2 ## [824,] 2 2 2 2 2 2 ## [825,] 2 2 2 2 2 2 ## [826,] 2 2 2 2 2 2 ## [827,] 2 2 2 2 2 2 ## [828,] 2 2 2 2 2 2 ## [829,] 2 2 2 2 2 2 ## [830,] 2 2 2 2 2 2 ## [831,] 2 2 2 2 2 2 ## [832,] 2 2 2 2 2 2 ## [833,] 2 2 2 2 2 2 ## [834,] 2 2 2 2 2 2 ## [835,] 2 2 2 2 2 2 ## [836,] 2 2 2 2 2 2 ## [837,] 2 2 2 2 2 2 ## [838,] 2 2 2 2 2 2 ## [839,] 2 2 2 2 2 2 ## [840,] 2 2 2 2 2 2 ## [841,] 2 2 2 2 2 2 ## [842,] 2 2 2 2 2 2 ## [843,] 2 2 2 2 2 2 ## [844,] 2 2 2 2 2 2 ## [845,] 2 2 2 2 2 2 ## [846,] 2 2 2 2 2 2 ## [847,] 2 2 2 2 2 2 ## [848,] 2 2 2 2 2 2 ## [849,] 2 2 2 2 2 2 ## [850,] 2 2 2 2 2 2 ## [851,] 2 2 2 2 2 2 ## [852,] 2 2 2 2 2 2 ## [853,] 2 2 2 2 2 2 ## [854,] 2 2 2 2 2 2 ## [855,] 2 2 2 2 2 2 ## [856,] 2 2 2 2 2 2 ## [857,] 2 2 2 2 2 2 ## [858,] 2 2 2 2 2 2 ## [859,] 2 2 2 2 2 2 ## [860,] 2 2 2 2 2 2 ## [861,] 2 2 2 2 2 2 ## [862,] 2 2 2 2 2 2 ## [863,] 2 2 2 2 2 2 ## [864,] 2 2 2 2 2 2 ## [865,] 2 2 2 2 2 2 ## [866,] 2 2 2 2 2 2 ## [867,] 2 2 2 2 2 2 ## [868,] 2 2 2 2 2 2 ## [869,] 2 2 2 2 2 2 ## [870,] 2 2 2 2 2 2 ## [871,] 2 2 2 2 2 2 ## [872,] 2 2 2 2 2 2 ## [873,] 2 2 2 2 2 2 ## [874,] 2 2 2 2 2 2 ## [875,] 2 2 2 2 2 2 ## [876,] 2 2 2 2 2 2 ## [877,] 2 2 2 2 2 2 ## [878,] 2 2 2 2 2 2 ## [879,] 2 2 2 2 2 2 ## [880,] 2 2 2 2 2 2 ## [881,] 2 2 2 2 2 2 ## [882,] 2 2 2 2 2 2 ## [883,] 2 2 2 2 2 2 ## [884,] 2 2 2 2 2 2 ## [885,] 2 2 2 2 2 2 ## [886,] 2 2 2 2 2 2 ## [887,] 2 2 2 2 2 2 ## [888,] 2 2 2 2 2 2 ## [889,] 2 2 2 2 2 2 ## [890,] 2 2 2 2 2 2 ## [891,] 2 2 2 2 2 2 ## [892,] 2 2 2 2 2 2 ## [893,] 2 2 2 2 2 2 ## [894,] 2 2 2 2 2 2 ## [895,] 2 2 2 2 2 2 ## [896,] 2 2 2 2 2 2 ## [897,] 2 2 2 2 2 2 ## [898,] 2 2 2 2 2 2 ## [899,] 2 2 2 2 2 2 ## [900,] 2 2 2 2 2 2 ## [901,] 2 2 2 2 2 2 ## [902,] 2 2 2 2 2 2 ## [903,] 2 2 2 2 2 2 ## [904,] 2 2 2 2 2 2 ## [905,] 2 2 2 2 2 2 ## [906,] 2 2 2 2 2 2 ## [907,] 2 2 2 2 2 2 ## [908,] 2 2 2 2 2 2 ## [909,] 2 2 2 2 2 2 ## [910,] 2 2 2 2 2 2 ## [911,] 2 2 2 2 2 2 ## [912,] 2 2 2 2 2 2 ## [913,] 2 2 2 2 2 2 ## [914,] 2 2 2 2 2 2 ## [915,] 2 2 2 2 2 2 ## [916,] 2 2 2 2 2 2 ## [917,] 2 2 2 2 2 2 ## [918,] 2 2 2 2 2 2 ## [919,] 2 2 2 2 2 2 ## [920,] 2 2 2 2 2 2 ## [921,] 2 2 2 2 2 2 ## [922,] 2 2 2 2 2 2 ## [923,] 2 2 2 2 2 2 ## [924,] 2 2 2 2 2 2 ## [925,] 2 2 2 2 2 2 ## [926,] 2 2 2 2 2 2 ## [927,] 2 2 2 2 2 2 ## [928,] 2 2 2 2 2 2 ## [929,] 2 2 2 2 2 2 ## [930,] 2 2 2 2 2 2 ## [931,] 2 2 2 2 2 2 ## [932,] 2 2 2 2 2 2 ## [933,] 2 2 2 2 2 2 ## [934,] 2 2 2 2 2 2 ## [935,] 2 2 2 2 2 2 ## [936,] 2 2 2 2 2 2 ## [937,] 2 2 2 2 2 2 ## [938,] 2 2 2 2 2 2 ## [939,] 2 2 2 2 2 2 ## [940,] 2 2 2 2 2 2 ## [941,] 2 2 2 2 2 2 ## [942,] 2 2 2 2 2 2 ## [943,] 2 2 2 2 2 2 ## [944,] 2 2 2 2 2 2 ## [945,] 2 2 2 2 2 2 ## [946,] 2 2 2 2 2 2 ## [947,] 2 2 2 2 2 2 ## [948,] 2 2 2 2 2 2 ## [949,] 2 2 2 2 2 2 ## [950,] 2 2 2 2 2 2 ## [951,] 2 2 2 2 2 2 ## [952,] 2 2 2 2 2 2 ## [953,] 2 2 2 2 2 2 ## [954,] 2 2 2 2 2 2 ## [955,] 2 2 2 2 2 2 ## [956,] 2 2 2 2 2 2 ## [957,] 2 2 2 2 2 2 ## [958,] 2 2 2 2 2 2 ## [959,] 2 2 2 2 2 2 ## [960,] 2 2 2 2 2 2 ## [961,] 2 2 2 2 2 2 ## [962,] 2 2 2 2 2 2 ## [963,] 2 2 2 2 2 2 ## [964,] 2 2 2 2 2 2 ## [965,] 2 2 2 2 2 2 ## [966,] 2 2 2 2 2 2 ## [967,] 2 2 2 2 2 2 ## [968,] 2 2 2 2 2 2 ## [969,] 2 2 2 2 2 2 ## [970,] 2 2 2 2 2 2 ## [971,] 2 2 2 2 2 2 ## [972,] 2 2 2 2 2 2 ## [973,] 2 2 2 2 2 2 ## [974,] 2 2 2 2 2 2 ## [975,] 2 2 2 2 2 2 ## [976,] 2 2 2 2 2 2 ## [977,] 2 2 2 2 2 2 ## [978,] 2 2 2 2 2 2 ## [979,] 2 2 2 2 2 2 ## [980,] 2 2 2 2 2 2 ## [981,] 2 2 2 2 2 2 ## [982,] 2 2 2 2 2 2 ## [983,] 2 2 2 2 2 2 ## [984,] 2 2 2 2 2 2 ## [985,] 2 2 2 2 2 2 ## [986,] 2 2 2 2 2 2 ## [987,] 2 2 2 2 2 2 ## [988,] 2 2 2 2 2 2 ## [989,] 2 2 2 2 2 2 ## [990,] 2 2 2 2 2 2 ## [991,] 2 2 2 2 2 2 ## [992,] 2 2 2 2 2 2 ## [993,] 2 2 2 2 2 2 ## [994,] 2 2 2 2 2 2 ## [995,] 2 2 2 2 2 2 ## [996,] 2 2 2 2 2 2 ## [997,] 2 2 2 2 2 2 ## [998,] 2 2 2 2 2 2 ## [999,] 2 2 2 2 2 2 ## [1000,] 2 2 2 2 2 2
D[1:5,1:5] # geographic distance
## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.0000000001 0.8114599329 0.6722476885 0.1654538663 0.6848265372 ## [2,] 0.8114599329 0.0000000001 0.2074216669 0.6495206211 0.4030935748 ## [3,] 0.6722476885 0.2074216669 0.0000000001 0.5068578162 0.2102744527 ## [4,] 0.1654538663 0.6495206211 0.5068578162 0.0000000001 0.5313442655 ## [5,] 0.6848265372 0.4030935748 0.2102744527 0.5313442655 0.0000000001
E[1:5,1:5] # ecological distance
## [,1] [,2] [,3] [,4] [,5] ## [1,] 1e-10 1e+00 1e-10 1e-10 1e-10 ## [2,] 1e+00 1e-10 1e+00 1e+00 1e+00 ## [3,] 1e-10 1e+00 1e-10 1e-10 1e-10 ## [4,] 1e-10 1e+00 1e-10 1e-10 1e-10 ## [5,] 1e-10 1e+00 1e-10 1e-10 1e-10
Initiating an MCMC run
MCMC(
counts = t(sim$genotypes),
sample_sizes = matrix(2,nrow=nind,ncol=nloci),
D = D, # geographic distances
E = E, # environmental distances
k = nind, loci = nloci, # dimensions of the data
delta = 0.0001, # a small, positive, number
aD_stp = 0.01, # step sizes for the MCMC
aE_stp = 0.01,
a2_stp = 0.025,
thetas_stp = 0.2,
mu_stp = 0.35,
ngen = 1000, # number of steps (2e6)
printfreq = 10000, # print progress (10000)
savefreq = 1000, # save out current state
samplefreq = 5, # record current state for posterior (2000)
prefix = "bedassle-ex/sim/CGW_1_", # filename prefix
continue=FALSE,
continuing.params=NULL)
## Warning in MCMC(counts = t(sim$genotypes), sample_sizes = matrix(2, nrow = ## nind, : the value chosen for delta may too large compared to the ## ecological distances between the sampled populations
## [1] "Output 1000 runs to bedassle-ex/sim/CGW_1_MCMC_output*.Robj ."
Examining a run
load("bedassle-ex/sim/CGW_1_MCMC_output1.Robj")
plot(aD, main="aD trace"); plot(aD_accept/aD_moves, main="aD acceptance", ylim=c(0,1))
plot(as.vector(aE), main="aE trace"); plot(aE_accept/aE_moves, main="aE acceptance", ylim=c(0,1))
plot(a2, main="a2 trace"); plot(a2_accept/a2_moves, main="a2 acceptance", ylim=c(0,1))
Examining a run
plot(Prob, xlab="MCMC generation", main="log likelihood") plot(mu_accept/mu_moves, xlab="MCMC generation", ylab="", main="mu acceptance", ylim=c(0,1) ) plot(thetas_accept/thetas_moves, xlab="MCMC generation", ylab="", main="thetas acceptance", ylim=c(0,1) )
Adjusting the parameters
make.continuing.params("bedassle-ex/sim/CGW_1_MCMC_output1.Robj",file.name="bedassle-ex/sim/CGW_1_MCMC_final1.Robj")
MCMC( counts = t(sim$genotypes), sample_sizes = matrix(2,nrow=nind,ncol=nloci),
D = D, E = E, k = nind, loci = nloci, delta = 0.0001,
aD_stp = 0.1, # increased from 0.01
aE_stp = 0.1, # increased from 0.01
a2_stp = 0.025, # kept the same
thetas_stp = 0.2, # kept the same
mu_stp = 0.35, # kept the same
ngen = 1000, printfreq = 10000, savefreq = 1000, samplefreq = 5,
prefix = "bedassle-ex/sim/CGW_2_", # NEW filename prefix
continue=TRUE, # CONTINUE from previous run
continuing.params="bedassle-ex/sim/CGW_1_MCMC_final1.Robj")
## Warning in MCMC(counts = t(sim$genotypes), sample_sizes = matrix(2, nrow = ## nind, : the value chosen for delta may too large compared to the ## ecological distances between the sampled populations
## [1] "Output 1000 runs to bedassle-ex/sim/CGW_2_MCMC_output*.Robj ."
Examining, again
Examining, again
Start again
make.continuing.params("bedassle-ex/sim/CGW_2_MCMC_output1.Robj",file.name="bedassle-ex/sim/CGW_2_MCMC_final1.Robj")
run.time <- system.time( MCMC( counts = t(sim$genotypes), sample_sizes = matrix(2,nrow=nind,ncol=nloci),
D = D, E = E, k = nind, loci = nloci, delta = 0.0001,
aD_stp = 0.1, # kept the same
aE_stp = 0.1, # kept the same
a2_stp = 0.02, # decreased from 0.025
thetas_stp = 0.2, # kept the same
mu_stp = 0.25, # decreased from 0.35
ngen = 1000, printfreq = 10000, savefreq = 1000, samplefreq = 5,
prefix = "bedassle-ex/sim/CGW_3_", # NEW filename prefix
continue=TRUE, # CONTINUE from previous run
continuing.params="bedassle-ex/sim/CGW_2_MCMC_final1.Robj") )
## Warning in MCMC(counts = t(sim$genotypes), sample_sizes = matrix(2, nrow = ## nind, : the value chosen for delta may too large compared to the ## ecological distances between the sampled populations
Examining, take three
Examining, take three
Ooops: step it back 
Find a good sampling rate
The autocorrelation function of the traces: correlation of parameter value as a function of number of intervening steps ("lag"):
acf(a0,lag.max=150,xlab='',ylab='')
Result: looks like sampling every \(50 \times 5 = 250\) steps would be reasonable.
Looks good
\(1000 \times 5\) iterations took 42.38 seconds … so let's set up a longer run:
make.continuing.params("bedassle-ex/sim/CGW_3_MCMC_output1.Robj",
file.name="bedassle-ex/sim/CGW_3_MCMC_final1.Robj")
save(sim, D, E, file="bedassle-ex/sim/sim-data.Robj")
Create this file, and run it (in the background, with Rscript), e.g. Rscript example-bedassle-script.R &
#!/usr/bin/Rscript
library(BEDASSLE)
setwd("bedassle-ex/sim")
load("sim-data.Robj")
MCMC( counts = t(sim$genotypes), sample_sizes = matrix(2,nrow=ncol(sim$genotypes),ncol=nrow(sim$genotypes)),
D = D, E = E, k = ncol(sim$genotypes), loci = nrow(sim$genotypes), delta = 0.0001,
aD_stp = 0.1, aE_stp = 0.1, a2_stp = 0.02, thetas_stp = 0.2, mu_stp = 0.25,
ngen = 1e6, # should take about 2.3 hours (on this machine)
printfreq = 10000, # no need to write out all the time
savefreq = 1e5, # save every 0.23 hours
samplefreq = 250, # as determined from acf plots
prefix = "CGW_4_", # NEW filename prefix
continue=TRUE, # CONTINUE from previous run
continuing.params="CGW_3_MCMC_final1.Robj")
Your turn
The results are saved as bedassle-ex/sim/CGW_4_MCMC_output1.Robj.
Suggestions:
- look at the run
- start a new (short!) run from this starting point
For results, look at
Summary and wrap-up
Spatial structure:
- Isolation by distance is driven by local migration and history;
- it is not clear how to translate island model results to continuous geography.
- First, look at the data.
PCA:
- Great for quick visualization,
- but beware of over-interpretation,
- as minor details are not stable, and major details might not be either.
BEDASSLE and family:
- Model-based, and closer to reality,
- but not a mechanistic model.
- Potentially much more information in the locus-specific residuals.
Technical appendices
Computational speed
Both BEDASSLE and PCA will slow down dramatically as the number of individuals grows.
They mostly do lots of linear algebra; this is faster if you use all of your computer's processors, instead of just one.
Do link R to multithreaded linear algebra libraries; as described on this page.
Trim linked SNPs and recode as numeric
Using plink:
# Trim down to pairwise r2 no more than 0.25: plink --file $INFILE --indep-pairwise 50 5 0.25 --out $INFILE_prunesnps # Turn result to numeric matrix plink --file $INFILE --extract $INFILE_prunesnps.prune.in --recodeA --out $INFILE_pruned
Estimate covariance matrix with R
In the presence of missing data:
cov(genotypes,use="pairwise.complete.obs")