PCA tutorial

[bryo-han]

@hyanwong (in #240) suggested somebody write a PCA tutorial.

I propose I'll send a PCA PR soon, once I have figured out how to use Jupyter Book. Does not look too hard.

Hannes

[hyanwong]

Amazing, thanks. I was just looking at this. In particular, I was interested in how to take a PCA result (e.g. obtained by averaging over diploid individuals, using ts.pca(individuals=XXX), and then taking a different sample node from the tree sequence (or even from a tsinfer matching path) and seeing where it landed in PCA space.

I've done this as follows, but I haven't really checked it through:

def get_PCs(focal_node, ts, individuals_used_in_pca, pca_output):
    comp = tts.genetic_relatedness(
        sample_sets=[[f]] + [list(x) for x in ts.individuals_nodes[individuals_used_in_pca,:]],
        indexes=[[0, i] for i in range(1, len(select_ind) + 1)],
        mode="branch", span_normalise=False)

## e.g.

select_ind = np.arange(ts.num_individuals - 1)
pca_output = ts.pca(num_components=2, individuals=select_ind)
plt.scatter(*pca_output.factors.T)
# Add an "X" for a different focal node
focal_nodes = ts1.individuals_nodes[-1,:]
plt.plot(* get_PCs(focal_nodes[0], ts, select_ind, pca_output), "x", c=popn_colours[-1])

[hannesbecher]

Not sure I understand. Thinking about a genotype matrix, one output of PCA/SVD/eigendecomposition is a matrix of 'loadings'. When the original (centred) genotypes are multiplied by the loadings matrix, we get the PC scores, of which we usually plot the first two cols.

Do you mean to transform previously unseen individuals into PC space? Based on variant data, you could get the covariance matrix of the makers (not individuals), get eigenvectors (the loadings matrix), and do GTs %*% eigenvectors to get the principal component scores (I think). But I don't know what to do with a tree sequence in that context. Probably best to rope in Peter R, maybe better to discuss on Slack?

[hyanwong]

Yes, transforming previously unseen individuals (or even haploid genomes) into the space. This could be done either by site measures, or by branch measures. In the latter case, of course, we need the branches of the new individual, so these could be obtain either because we have the ARG, but the original PCA was only run on a subset, or by threading the new sample into an existing ARG.

[hannesbecher]

Makes sense. Not sure how that would be done with branches.

Have sent a PR (#301), rendering works now.