GenotypeDisplay outputted to Matrix/DataFrame

[zwallace425]

Hi,

Is there a simple way to get the output for sgkit.display_genotype into a dataframe or numpy array? I want each variant_id as a row, each sample_id as a column, and each entry to be the genotype (in may case 0/0, 1/0, 1/1, or ./.). The display_genotype method shows me exactly what I want, but I want that outputted to a matrix. Since the output is a 'GenotypeDisplay' object, typical python packages aren't letting me work with it appropriately.

Thank you for the help,
Zach

[tomwhite]

Hi Zach,

Calling sgkit.display_genotype(ds).df should return a pandas dataframe with what you want. The GenotypeDisplay object is usually used just for its string representation, but you should be able to get its underlying data even though that's not a documented part of the API.

To get the genotype data as an integer NumPy array you would call ds["call_genotype"].values. This is the preferred way to access genotype data for further processing.

Hope that helps,
Tom

[zwallace425]

Hi Tom,

Thank you! This is helps a lot! However, I see that display_genotypes(ds).df is WAY slower than ds["call_genotypes"].values, but I think the end product of display_genotypes is really nice considering it's a DataFrame with row and column labels that can make it straightforward for any downstream applications (ML/visualizations). That said, you mention that GenotypeDisplay is used for string representation, which makes sense considering the 0/0, 1/0, 1/1, ./. entries in the display data frame are all strings representations of the genotype. But is there a way in sgkit to get those strings represented as integers, say 0/0 = 0, 1/0 = 1, 1/1 = 2, and ./. = -1 without me having to hard code that from the results of the output matrix? No worries if not, I'm just generally curious in how to get these genotypes into machine readable data frames and I'm still trying to learn the components of sgkit!

Thank you for the help,
Zach

[tomwhite]

Hi Zach,

You are absolutely right that it's better to work with integers rather than strings for downstream analysis. Sgkit datastructures use xarray, so you can use its API for manipulation, including converting final results to pandas or pydata visualization libraries.

The call_genotypes variable has three dimensions: variants, samples, and ploidy, and what you want to do is aggregate across ploidy. Here's one way to do it:

In [1]: import sgkit as sg

In [2]: ds = sg.simulate_genotype_call_dataset(n_variant=10, n_sample=5, missing_pct=.2)

In [3]: ds
Out[3]: 
<xarray.Dataset>
Dimensions:             (variants: 10, alleles: 2, samples: 5, ploidy: 2)
Dimensions without coordinates: variants, alleles, samples, ploidy
Data variables:
    variant_contig      (variants) int64 0 0 0 0 0 0 0 0 0 0
    variant_position    (variants) int64 0 1 2 3 4 5 6 7 8 9
    variant_allele      (variants, alleles) |S1 b'T' b'A' b'G' ... b'C' b'C'
    sample_id           (samples) <U2 'S0' 'S1' 'S2' 'S3' 'S4'
    call_genotype       (variants, samples, ploidy) int8 0 0 1 0 1 ... 0 0 0 1 1
    call_genotype_mask  (variants, samples, ploidy) bool False False ... False
Attributes:
    contigs:  [0]
    source:   sgkit-0.6.1.dev2+gcc728043

In [4]: sg.display_genotypes(ds)
Out[4]: 
samples    S0   S1   S2   S3   S4
variants                         
0         0/0  1/0  1/0  ./1  1/0
1         0/1  0/0  1/0  ./1  ./0
2         1/.  0/.  1/0  1/1  1/1
3         ./0  1/.  0/0  1/0  ./1
4         1/1  1/1  0/0  0/0  0/1
5         1/0  1/0  ./0  0/0  1/0
6         ./1  ./0  0/0  0/0  0/1
7         0/1  0/1  1/0  ./0  ./0
8         0/1  0/0  1/.  0/.  1/0
9         1/1  0/0  1/0  0/0  1/1

In [5]: import xarray as xr

In [6]: call_g_mask = ds["call_genotype_mask"].any(dim="ploidy")

In [7]: call_g = xr.where(call_g_mask, -1, ds["call_genotype"].sum(dim="ploidy"))

In [8]: call_g.values
Out[8]: 
array([[ 0,  1,  1, -1,  1],
       [ 1,  0,  1, -1, -1],
       [-1, -1,  1,  2,  2],
       [-1, -1,  0,  1, -1],
       [ 2,  2,  0,  0,  1],
       [ 1,  1, -1,  0,  1],
       [-1, -1,  0,  0,  1],
       [ 1,  1,  1, -1, -1],
       [ 1,  0, -1, -1,  1],
       [ 2,  0,  1,  0,  2]])

The mask is needed to take into account missing calls. There are more examples like this in the sgkit docs at https://pystatgen.github.io/sgkit/latest/getting_started.html, and also https://pystatgen.github.io/sgkit/latest/examples/gwas_tutorial.html.

This operation to aggregate/collapse across the ploidy dimension seems fairly common. We do it in the sgkit code, but don't expose it as a user function. @timothymillar any thoughts about if this would be useful/what an API might look like if so?

[timothymillar]

I can see a couple of options. The first is to add a method which converts call_genotype or call_allele_count to call_dosage. However, I don't like this idea as call_dosage is only a suitable abstraction for biallelic data and it's intentionally flexible in its interpretation (could be the dose of the reference or alternate allele, can also represent probabilities as floats etc.).

The second option is to add a method which converts call_genotypes to their VCF-sort-order index. This is essentially a perfect integer hash for every unique genotpye of a given ploidy. Conveniently, it's just the count of the alternate allele for biallelic genotypes. We already have an implementation of the algorithm, so it would just be a matter of making a public method. The hard part will be naming that method...

[zwallace425]

Tom and Timothy, apologies for my delay, but thanks so much for your responses last week!

Tom - that second chunk of code for converting the genotypes to an integer matrix was exactly what I wanted and ended up working very efficiently for me when I ran it. I've been implementing it into my pipelines.

Timothy - I'm trying to clearly understand what you mean in the options your presented, as I'm somewhat confused. My data is indeed strictly biallelic, so how would this "call_dosage" method be aggregating genotypes as opposed to the -1, 0, 1, and 2 imputation that Tom showed above? It sounds like these could be floating point representations? And would either of these methods be more suitable for imputing the unknown genotypes (the -1's)?

[timothymillar]

Hi @zwallace425. Don't worry about my comment, it was just some ideas about how we could create a new function in sgkit that would also work for multiallelic data. If you are only have biallelic data then Toms code will be fine.

[tomwhite]

I opened #1048 to track the new function.

[hammer]

Should we also update the docs to make it easy for new users to find these methods?