Assistance with filtering and relatatedness analysis

[Molly Barlow]

Hi all, I'm attempting some simple relatedness and inbreeding analysis of a fenced population of marsupial, but my relatedness networks are indicating a much higher degree of relatedness than they should be (i.e. suggesting they're all identical). I was also hoping for some advice on my filtering procedure. 

I've sent my code directly to Luis. 

Cheers,

Molly 

[Jose Luis Mijangos]

Hi Molly,

Thanks for sending your code, your filtering already looks great. Just a few comments.

- Remove individuals with low call rate first, with gl.filter.callrate(x, method = "ind"). Doing this first keeps SNPs you would otherwise lose.
- Use gl.filter.rdepth(x, lower = 5, upper = 1000). Loci with low depth are often miscalled as homozygotes, so this gives you more accurate heterozygote calls.
- Use gl.filter.maf(x, threshold = 0.05). Rare alleles add noise to relatedness estimates and are frequently sequencing errors.

On the relatedness estimates themselves, the thing to keep in mind is that relatedness is relative, not absolute. It is estimated against a reference population, which in practice is your own dataset, and estimators assume that reference population is outbred and randomly mating. If several samples in the dataset are related or inbred, estimates are biased upwards. If the dataset has samples from two differentiated populations, they are biased downwards.

Methods handle this differently. gl.grm builds a genomic relationship matrix (GRM) centred on your sample, so 0 is not an unrelated pair but a pair of average relatedness for your dataset. In a dataset like yours, descended from a few founders, the average pair already shares a lot of ancestry, so a GRM value of 0 already sits above true unrelatedness.

I would use EMIBD9 instead, via gl.run.EMIBD9 with Inbreed = TRUE. It estimates allele frequencies and relatedness coefficients jointly by maximum likelihood, which partly corrects for the problem above.

Also note that relatedness and kinship have different definitions. Relatedness is the expected proportion of alleles that are identical by descent between two individuals, so full siblings have a theoretical relatedness value of 0.5. Kinship is the probability that two alleles, one drawn at random from each individual, are identical by descent, so full siblings have a theoretical kinship value of 0.25 (i.e. kinship is half of relatedness).

 EMIBD9 reports kinship, and you can feed it straight into gl.grm.network, as shown below.

library(dartRverse)
library(viridis)
programs_path <- "path_to_folder_containing_EMIBD9"

t1 <- platypus.gl
EMIBD9 <- gl.run.EMIBD9(t1,
                        emibd9.path = programs_path,
                        Inbreed = TRUE)
grm_sib <- gl.grm.network(G = EMIBD9$rel,
                          x = t1,
                          categorise = TRUE,
                          node.size = 8,
                          palette_discrete = viridis::viridis_pal(),
                          kinship.threshold = 0.125,
                          legend.title = "Sampling locations",
                          title = "")

You can download the EMIBD9 binary from here: https://github.com/green-striped-gecko/dartRverse/tree/main/binaries

Happy to help further if anything is unclear.

Cheers,
Luis

[Jose Luis Mijangos]

Hi Molly,

A correction to my last reply. I had the two bias directions swapped.

Many close relatives in the dataset bias estimates for related pairs downwards, not upwards. The alleles a family shares become common in the dataset, so estimators treat that sharing as ordinary rather than as evidence of relatedness. EMIBD9 corrects this one.

Mixing differentiated populations biases estimates within each population upwards, not downwards. Alleles common in one population but rare in the combined dataset make individuals from that population look like relatives. EMIBD9 does not correct this one, so run it separately for each population or sampling location.

Cheers,
Luis

[Molly Barlow]

Thanks Luis. I incorporated your suggestions on my read-depth and MAF filtering. Thanks for those suggestions. I also attempted to remove individuals with low call rates first as per your suggestion: I filtered at 0.95 (the default) and retained only 10 of my 94 individuals. I tried a few different/ lower thresholds but found it was biased in its removal, i.e. I’d lose most of my founder individuals, which are individuals that are important for further analysis.

Thanks also for the advice regarding gl.grm v EMIBD9. I generated a network using EMIBD9 and although it still looks highly related, it looks a bit more like what we expect that the gl.grm network 

I'm also playing around with Ho, uHe and FIS estimates. Should I be filtering differently to how I have for kinship/ relatedness estimates. I.e. is it better to not filter MAF? Or include lower read depth? 

Thanks

[Jose Luis Mijangos]

Hi Molly,

Yes, Ho, uHe and FIS require different filtering. The main recommendations are:

- Don’t use a MAF filter, because removing rare alleles can inflate heterozygosity estimates.
- Use the read-depth filter. A minimum depth of >8 is a good starting point, as low depth can cause heterozygotes to be miscalled as homozygotes.
- Use a strict locus call-rate threshold, ideally 1, because loci with missing data can give unreliable heterozygosity estimates.
- Use equal subsamples when comparing groups, because heterozygosity estimates depend on sample size. You can test this with `gl.report.heterozygosity(x, subsample.pop = TRUE)`.

Also, many factors influence heterozygosity estimates, so they are best compared within your dataset rather than across different studies, datasets, or species.

Here are a couple of useful papers with more detail:

- Schmidt, T. L., Jasper, M.-E., Weeks, A. R., & Hoffmann, A. A. (2021). Unbiased population heterozygosity estimates from genome-wide sequence data. *Methods in Ecology and Evolution*, 12(10), 1888–1898. https://doi.org/10.1111/2041-210X.13659
- Sopniewski, J., & Catullo, R. A. (2024). Estimates of heterozygosity from single nucleotide polymorphism markers are context-dependent and often wrong. *Molecular Ecology Resources*, 24(4), e13947. https://doi.org/10.1111/1755-0998.13947

Cheers,
Luis