PCA tutorial by hannesbecher · Pull Request #301 · tskit-dev/tutorials

hannesbecher · 2026-03-25T12:11:53Z

Hi,
here is my proposed PCA tutorial. This covers:

PCA on genotypes vs branches
haplotypes vs diploid genotypes
simulated data only

Looking forward to opinions/comments/requests!

HELP please: for some reason the cell outputs of this tutorial don't render when I run make on my system (but they do for all other existing python tutorials). What I might be going on? I had written this as an ipynb and then converted to md.

Cheers,
Hannes

hannesbecher · 2026-03-25T15:56:15Z

OK, have figured it out. format_name in the md file has to be format_name: myst, it seems.

hyanwong · 2026-03-26T03:31:19Z

PCA.md

+
+```{code-cell} ipython3
+# load required libraries
+import msprime as msp


I think I would prefer simply import msprime, which matches the rest of the documentation

hyanwong · 2026-03-26T03:33:17Z

PCA.md

+fstmat
+```
+
+## PCA 'by hand'


Should we perhaps give a quick demo of ts.pca first, in case people just want to copy/paste the code? I know that doesn't quite fit with your logical flow, however.

I agree that would be better tutorial-wise, simple approach first, then details

hyanwong · 2026-03-26T03:34:39Z

A few comments, one minor, the other a bit more major (should we demo the built-in method quickly first?). Happy to merge either way, and we can change later.

One thing though: could you possibly squash all the commits into one? Thanks @hannesbecher !

gregorgorjanc · 2026-03-26T07:39:38Z

PCA.md

+                      demography=dmg,
+                      random_seed=1234,
+                      sequence_length=1e6,
+                      recombination_rate=1e-9)


Maybe increase to the same order as mut rate?

Good, makes the plots look clearer.

gregorgorjanc · 2026-03-26T07:41:43Z

README.md

 This command should do the trick:
 ```
-$ R -e 'install.packages(c("reticulate", "IRkernel")); IRkernel::installspec()'
+$ R -e 'options(repos=c(CRAN="http://cran.r-project.org")); install.packages(c("reticulate", "IRkernel")); IRkernel::installspec()'


Isn't this the default anyway?

Failed on my Mac without that addition. I guess it is not going to do any harm?

gregorgorjanc · 2026-03-26T07:46:23Z

PCA.md

+fstmat = np.zeros([nPop,nPop])
+for i in range(nPop-1):
+    for j in range(i+1,nPop):
+        fstmat[i,j] = tsm.Fst([range(i*nHap,(i+1)*nHap), range((i+1)*nHap,(i+2)*nHap)])


Would there be value in computing branch Fst too? I assume that below you will compute branch and site PCA ...

I would leave this for the stats tutorial(s).

PCA.md

gregorgorjanc · 2026-03-26T07:51:07Z

PCA.md

+
+```{code-cell} ipython3
+# haplotypes, each sample haplotype is ues by default
+hapBranchPca=ts.pca(num_components=10)


This is a nitpick, but above you used ht/gt to denote haplotypes and genotypes, so maybe use the same style for results too, htBranchPca etc.

gregorgorjanc · 2026-03-26T07:52:36Z

PCA.md

+                  c=np.repeat([1,2,3,4,5], [nHap] * nPop))
+axs[0, 0].set_title('Haplotypes (sites)')
+
+axs[0,0].set_ylabel("PC2")


Do you need to add PC1 axis label too?

Maybe not ... doing this in a browser;)

I think it's beautiful as it is, GG.

gregorgorjanc · 2026-03-26T07:54:47Z

PCA.md

+
+The plots on the left show one dot per haplotype. These have twice as many dots as the plots on the right, which show individuals. The colours indicate from which of the five islands a haplotype or individual was sampled. As expected with low geneflow, there is some grouping by island. Feel free to re-run with higher or lower values of `migRate` to see how the separations between the island samples changes.
+
+++


I quite like it. But happy should somebody send a PR...

PCA.md

gregorgorjanc · 2026-03-26T08:01:24Z

PCA.md

+## Empirical data
+Here, we demonstrated using simulated data how SNPs and ARG branches lead to equivalent PCA results. For empirical data, the ancestral states of variant sites are not known a priori, which will in practice often lead to polarisation differences. That may affect the outcome of PCA.
+
+**TODO:** Extend Tutorial to empirical data.


It might be simplest to demonstrate this by randomly? repolarising SNP matrices and repeating the PCA on these?

Or major allele -> ancestral, etc. Random as in 50:50, did not make much difference as far as I remember from the randPedPCA project. Needs more planning.

hannesbecher · 2026-03-26T09:34:46Z

A few comments, one minor, the other a bit more major (should we demo the built-in method quickly first?). Happy to merge either way, and we can change later.

One thing though: could you possibly squash all the commits into one? Thanks @hannesbecher !

Absolutely. I can do this with the merge button, though, correct, @hyanwong ?

Added tutorial on comparing branch and SNP-based PCA using msprime and tskit, including code examples for simulating ARGs and performing PCA. Add PCA file to the table of contents Modify R installation command in README Updated R package installation command to specify CRAN repository. Refactor PCA.md for clarity and organization Updated kernel specifications and code formatting for clarity. Added comments and organized code cells for better readability. Revise PCA tutorial and update Jupyter metadata Updated Jupyter Notebook metadata and improved PCA tutorial content. Update kernelspec display name and language in PCA.md Update PCA.md Co-authored-by: Gregor Gorjanc <gregor.gorjanc@gmail.com> Update PCA.md Co-authored-by: Gregor Gorjanc <gregor.gorjanc@gmail.com> Update PCA.md Co-authored-by: Gregor Gorjanc <gregor.gorjanc@gmail.com> PR comments implemented heading

hannesbecher marked this pull request as draft March 25, 2026 14:08

hannesbecher marked this pull request as ready for review March 25, 2026 15:54

hannesbecher mentioned this pull request Mar 25, 2026

PCA tutorial #300

Open

hyanwong approved these changes Mar 26, 2026

View reviewed changes

gregorgorjanc reviewed Mar 26, 2026

View reviewed changes

PCA.md Outdated Show resolved Hide resolved

gregorgorjanc reviewed Mar 26, 2026

View reviewed changes

PCA.md Outdated Show resolved Hide resolved

gregorgorjanc reviewed Mar 26, 2026

View reviewed changes

PCA.md Show resolved Hide resolved

gregorgorjanc reviewed Mar 26, 2026

View reviewed changes

bryo-han force-pushed the main branch from c5fd7cc to 7ec0128 Compare March 27, 2026 10:41


		The plots on the left show one dot per haplotype. These have twice as many dots as the plots on the right, which show individuals. The colours indicate from which of the five islands a haplotype or individual was sampled. As expected with low geneflow, there is some grouping by island. Feel free to re-run with higher or lower values of `migRate` to see how the separations between the island samples changes.

		+++

Conversation

hannesbecher commented Mar 25, 2026

Uh oh!

hannesbecher commented Mar 25, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hyanwong commented Mar 26, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hannesbecher Mar 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hannesbecher commented Mar 26, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

hannesbecher Mar 26, 2026 •

edited

Loading