Selecting only a subset of individuals for writing the vcf but getting n-1 individuals

[sroy0802]

I have simulated my samples with SLiM and am trying to write the vcf file with only a subset of samples (100) with their sample ids. The issue is that I'm getting 99 (n-1) samples written in the vcf, everytime. Here is how I am doing it. I'm subsetting the tree and then writing the samples with them. Is there a more efficient way of doing it? Is there any particular reason I am getting n-1 samples instead?

writing a new vcf file with just 100 highly admixed samples

with open("path/to/sample_list", "r") as f:
subset_samples = [int(line.strip()) for line in f]

Assigning names to individuals

##individual_names = [f"Ind_{i}" for i in subset_samples]

Filtering the tree sequence to keep only selected samples

subset_ts = ts.simplify(samples=subset_samples, keep_unary=True)
print(subset_ts.samples())

with open("path/to/vcf", "w") as vcf_file:
subset_ts.write_vcf(vcf_file)

I would highly appreciate some help and suggestions.
Thanks,
Susweta

[molpopgen]

One possible problem may be that your "subset samples" are Individual ids and not Node ids? The latter is required for simplification. But it is hard to tell w/o knowing what the contents of the file that you are opening.

[sroy0802]

I tried verifying that and the sample_list are indeed sample ids and not node ids. It is a single column file containing a list of 100 samples. I have selected these ids from tspop output. If I understand it right I should be obtaining the node ids with: tspop.get_pop_ancestry(ts, census_time=n).ancestry_table
Please correct me if I am wrong.
I also verified if all the samples have genotype information and/or at least one ancestral allele and it seems like all of them does. Is there any other checks I can perform?

[molpopgen]

For a list of individual ids, you can access their rows in the individual table to extract nodes, and then concatenate ("flatten") the node ids into a 1-d array.

However, as @benjeffery points out, we need a minimal example that reproduces your issue. Without that we have to make guesses.

[benjeffery]

Please can you create a self-contained minimal example? It is hard to understand what you are doing with the snippets.

[sroy0802]

Yes, sure!
#Step 1: An example SLiM simulation
initialize()
{
initializeTreeSeq();
initializeMutationRate(0);
initializeMutationType("m1", 0.5, "f", 0.0);
initializeGenomicElementType("g1", m1, 1.0);
initializeGenomicElement(g1, 0, 100 * 1000000 - 1);
initializeRecombinationRate(1e-8);
}
// Here ancestral populations simulates for 1 generation
1 early()
{
sim.addSubpop("p1", 500);
sim.addSubpop("p2", 1000);
}

// In the 2nd generation, admixed population is generated entirely from migration from the ancestral population
2 early()
{
sim.addSubpop("p3", 500);
p3.setMigrationRates(c(p1, p2), c(0.75, 0.25));
// CENSUS TIME
sim.treeSeqRememberIndividuals(sim.subpopulations.individuals);
}

// From 3rd generation ancestral population is not needed.
2 late()
{
// removing the source subpopulations
p3.setMigrationRates(c(p1, p2), c(0.0, 0.0));
p1.setSubpopulationSize(0);
p2.setSubpopulationSize(0);
}

100 late()
{
sim.simulationFinished();
sim.treeSeqOutput("path/test.trees");
}

#Step2: using tspop to get the ancestry proportions
import tskit
import pyslim
import tspop
import numpy as np
import pandas as pd

ts = tskit.load("path/test.trees")
pa = tspop.get_pop_ancestry(ts, census_time=3)
ancestry_df = pa.ancestry_table

#Selecting a the first 20 samples from this table for this example
ancestry_df.sort_values(['sample', 'left'], inplace=True)
screened_df = ancestry_df.groupby("sample")["population"]
.apply(lambda x: ",".join(map(str, sorted(set(x))))).reset_index()
screened_df.head(20)to_csv(path/test_samples, sep="\t", index=False)

#Step3: rest of the code snippet provided

writing a new vcf file with just 100 highly admixed samples

with open("path/test_samples", "r") as f:
subset_samples = [int(line.strip()) for line in f]

Assigning names to individuals

##individual_names = [f"Ind_{i}" for i in subset_samples]

Filtering the tree sequence to keep only selected samples

subset_ts = ts.simplify(samples=subset_samples, keep_unary=True)
print(subset_ts.samples())

with open("path/to/vcf", "w") as vcf_file:
subset_ts.write_vcf(vcf_file)

[petrelharp]

Is the problem that you are expecting subset_samples to be a vector of individual IDs? The argument to simplify is a list of node IDs, where a "node" is one of the two chromosomes (for diploids). For more discussion, also see this bit of the pyslim docs.

[hyanwong]

What does len(np.unique(subset_samples)) and subset_ts.num_samples say?

[hyanwong]

Right, so one individual (id 80) is diploid, and the others are all haploid. That's why. You can see the node ids by subset_ts.individual(80)

Yan

[sroy0802]

Thank you so much Yan! Although, I am confused. I tried checking if my simulation is producing diploid or individuals with
samples = ts.samples()
print(f"Number of haplotypes (genomes): {len(samples)}")
print(f"Number of individuals: {ts.num_individuals}")

Turns out, it is producing diploid individuals but the other individual ids (except id 80 here) are being presented as haploids. So, how do I fix this issue? I have also provided an example code to mimic what I'm doing. I would highly appreciate any suggestions.

[hyanwong]

I'm afraid I'm no expert in SLiM, but I suspect that when you take a census, it is taken at some point in the past. You are not guaranteed that individuals in the past will have both their genomes as ancestral to individuals in the present. I.e. if you take a census, you are sampling genomes not individuals. You can't expect to get the full genome of individuals in the past unless you mark them for remembering in SLiM. And even then, those historical genomes will certainly contain regions that are not ancestral to anyone in the present.

EDIT - actually - looking at your code, I think it's not that you are sampling for a population in the past, but that you are selecting 100 sample genomes at random. Perhaps you want to sample 50 diploid individuals at random instead? Or maybe 100 diploid individuals (i.e. 200 sample genomes)?

I don't think you have show the code where you take the random sample, but you want to be sampling individuals, not genomes (sample nodes). You can then pass the nodes of all the individuals (which will be twice the size) to simplify()

[petrelharp]

I think this last point is the problem - see my answer above. Said another way, you're passing in a list of 100 nodes (i.e., chromosomes) to simplify; and by chance two of those happen to be in the same individual; so what happens is that you get 98 individuals with only one chromosome left and one individual with both chromosomes.

[hyanwong]

This sort of error happens often enough that I wonder if we should issue a (silenceable) UserWarning if simplify creates a set of individuals with non-matching ploidies? I imagine that's not the intention most of the time.

[sroy0802]

Thank you @petrelharp @hyanwong! that makes sense now. The pyslim documentation helps. I'm trying to see if I can make it work. Once again, this has been an extremely helpful discussion! I might come back and badger you with more queries in a few hours. Thank you so much for your patience.

[petrelharp]

This sort of error happens often enough that I wonder if we should issue a (silenceable) UserWarning if simplify creates a set of individuals with non-matching ploidies? I imagine that's not the intention most of the time.

That's a good idea. The thing I've had in mind that might address this is to add an individuals argument.