How do I …#

Create a test dataset?#

Call simulate_genotype_call_dataset() to create a test xarray.Dataset:

In [1]: import sgkit as sg

In [2]: ds = sg.simulate_genotype_call_dataset(n_variant=100, n_sample=50, n_contig=23, missing_pct=.1)

Look at the dataset summary?#

Print using the xarray.Dataset repr:

In [3]: ds
Out[3]: 
<xarray.Dataset>
Dimensions:             (contigs: 23, variants: 100, alleles: 2, samples: 50,
                         ploidy: 2)
Dimensions without coordinates: contigs, variants, alleles, samples, ploidy
Data variables:
    contig_id           (contigs) <U2 '0' '1' '2' '3' ... '19' '20' '21' '22'
    variant_contig      (variants) int64 0 0 0 0 0 1 1 ... 21 21 21 22 22 22 22
    variant_position    (variants) int64 0 1 2 3 4 0 1 2 3 ... 3 0 1 2 3 0 1 2 3
    variant_allele      (variants, alleles) |S1 b'T' b'C' b'C' ... b'T' b'A'
    sample_id           (samples) <U3 'S0' 'S1' 'S2' 'S3' ... 'S47' 'S48' 'S49'
    call_genotype       (variants, samples, ploidy) int8 0 0 1 0 1 ... 1 0 1 0 0
    call_genotype_mask  (variants, samples, ploidy) bool False False ... False
Attributes:
    contigs:  ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', ...
    source:   sgkit-unknown

Get the values for a variable in a dataset?#

Call xarray.Variable.values:

In [4]: ds.variant_contig.values
Out[4]: 
array([ 0,  0,  0,  0,  0,  1,  1,  1,  1,  1,  2,  2,  2,  2,  2,  3,  3,
        3,  3,  3,  4,  4,  4,  4,  4,  5,  5,  5,  5,  5,  6,  6,  6,  6,
        6,  7,  7,  7,  7,  7,  8,  8,  8,  8,  9,  9,  9,  9, 10, 10, 10,
       10, 11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 14,
       15, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 17, 18, 18, 18, 18, 19,
       19, 19, 19, 20, 20, 20, 20, 21, 21, 21, 21, 22, 22, 22, 22])

In [5]: ds["variant_contig"].values # equivalent alternative
Out[5]: 
array([ 0,  0,  0,  0,  0,  1,  1,  1,  1,  1,  2,  2,  2,  2,  2,  3,  3,
        3,  3,  3,  4,  4,  4,  4,  4,  5,  5,  5,  5,  5,  6,  6,  6,  6,
        6,  7,  7,  7,  7,  7,  8,  8,  8,  8,  9,  9,  9,  9, 10, 10, 10,
       10, 11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 14,
       15, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 17, 18, 18, 18, 18, 19,
       19, 19, 19, 20, 20, 20, 20, 21, 21, 21, 21, 22, 22, 22, 22])

Warning

Calling values materializes a variable’s data in memory, so is only suitable for small datasets.

Find the definition for a variable in a dataset?#

Use the comment attribute on the variable:

In [6]: ds.variant_contig.comment
Out[6]: 'Index corresponding to contig name for each variant. In some less common\nscenarios, this may also be equivalent to the contig names if the data\ngenerating process used contig names that were also integers.'

All the variables defined in sgkit are documented on the Variables API page.

Look at the genotypes?#

Call display_genotypes():

In [7]: sg.display_genotypes(ds, max_variants=10)
Out[7]: 
samples    S0   S1   S2   S3   S4  ...  S45  S46  S47  S48  S49
variants                           ...                         
0         0/0  1/0  1/0  0/1  1/0  ...  1/1  0/0  1/0  0/0  1/1
1         1/1  1/0  1/.  ./0  1/0  ...  1/1  0/1  1/0  1/1  1/0
2         0/1  1/1  1/1  1/0  1/1  ...  0/0  0/1  0/0  0/0  1/1
3         1/1  0/0  1/1  ./1  0/1  ...  0/1  1/0  0/1  0/.  0/.
4         1/0  0/1  0/1  0/1  0/0  ...  1/0  1/1  0/0  1/.  1/0
...       ...  ...  ...  ...  ...  ...  ...  ...  ...  ...  ...
6         1/1  1/1  ./0  1/1  0/1  ...  0/0  0/.  1/0  1/0  0/1
7         1/.  1/0  ./0  0/1  1/0  ...  0/1  1/.  0/0  1/0  0/0
8         0/1  0/0  0/0  0/1  0/0  ...  0/1  0/1  1/0  1/0  0/0
9         1/1  0/0  ./1  1/0  0/0  ...  0/0  0/0  1/1  0/1  1/0
10        1/1  0/.  0/0  0/1  1/.  ...  1/0  0/.  0/1  0/1  0/0

[100 rows x 50 columns]

Subset the variables?#

Use Xarray’s pandas-like method for selecting variables:

In [8]: ds[["variant_contig", "variant_position", "variant_allele"]]
Out[8]: 
<xarray.Dataset>
Dimensions:           (variants: 100, alleles: 2)
Dimensions without coordinates: variants, alleles
Data variables:
    variant_contig    (variants) int64 0 0 0 0 0 1 1 1 ... 21 21 21 22 22 22 22
    variant_position  (variants) int64 0 1 2 3 4 0 1 2 3 4 ... 3 0 1 2 3 0 1 2 3
    variant_allele    (variants, alleles) |S1 b'T' b'C' b'C' ... b'T' b'T' b'A'
Attributes:
    contigs:  ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', ...
    source:   sgkit-unknown

Alternatively, you can drop variables that you want to remove:

In [9]: ds.drop_vars(["variant_contig", "variant_position", "variant_allele"])
Out[9]: 
<xarray.Dataset>
Dimensions:             (contigs: 23, samples: 50, variants: 100, ploidy: 2)
Dimensions without coordinates: contigs, samples, variants, ploidy
Data variables:
    contig_id           (contigs) <U2 '0' '1' '2' '3' ... '19' '20' '21' '22'
    sample_id           (samples) <U3 'S0' 'S1' 'S2' 'S3' ... 'S47' 'S48' 'S49'
    call_genotype       (variants, samples, ploidy) int8 0 0 1 0 1 ... 1 0 1 0 0
    call_genotype_mask  (variants, samples, ploidy) bool False False ... False
Attributes:
    contigs:  ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', ...
    source:   sgkit-unknown

Subset to a genomic range?#

Set an index on the dataset, then call xarray.Dataset.sel():

In [10]: ds.set_index(variants=("variant_contig", "variant_position")).sel(variants=(0, slice(2, 4)))
Out[10]: 
<xarray.Dataset>
Dimensions:             (contigs: 23, variants: 3, alleles: 2, samples: 50,
                         ploidy: 2)
Coordinates:
  * variants            (variants) object MultiIndex
  * variant_contig      (variants) int64 0 0 0
  * variant_position    (variants) int64 2 3 4
Dimensions without coordinates: contigs, alleles, samples, ploidy
Data variables:
    contig_id           (contigs) <U2 '0' '1' '2' '3' ... '19' '20' '21' '22'
    variant_allele      (variants, alleles) |S1 b'T' b'G' b'G' b'G' b'C' b'G'
    sample_id           (samples) <U3 'S0' 'S1' 'S2' 'S3' ... 'S47' 'S48' 'S49'
    call_genotype       (variants, samples, ploidy) int8 0 1 1 1 1 ... 1 -1 1 0
    call_genotype_mask  (variants, samples, ploidy) bool False False ... False
Attributes:
    contigs:  ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', ...
    source:   sgkit-unknown

An API to make this easier is under discussion. Please add your requirements to pystatgen/sgkit#658.

Get the list of samples?#

Get the values for the sample_id variable:

In [11]: ds.sample_id.values
Out[11]: 
array(['S0', 'S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10',
       'S11', 'S12', 'S13', 'S14', 'S15', 'S16', 'S17', 'S18', 'S19',
       'S20', 'S21', 'S22', 'S23', 'S24', 'S25', 'S26', 'S27', 'S28',
       'S29', 'S30', 'S31', 'S32', 'S33', 'S34', 'S35', 'S36', 'S37',
       'S38', 'S39', 'S40', 'S41', 'S42', 'S43', 'S44', 'S45', 'S46',
       'S47', 'S48', 'S49'], dtype='<U3')

Subset the samples?#

Call xarray.Dataset.sel() and xarray.DataArray.isin():

In [12]: ds.sel(samples=ds.sample_id.isin(["S30", "S32"]))
Out[12]: 
<xarray.Dataset>
Dimensions:             (contigs: 23, variants: 100, alleles: 2, samples: 2,
                         ploidy: 2)
Dimensions without coordinates: contigs, variants, alleles, samples, ploidy
Data variables:
    contig_id           (contigs) <U2 '0' '1' '2' '3' ... '19' '20' '21' '22'
    variant_contig      (variants) int64 0 0 0 0 0 1 1 ... 21 21 21 22 22 22 22
    variant_position    (variants) int64 0 1 2 3 4 0 1 2 3 ... 3 0 1 2 3 0 1 2 3
    variant_allele      (variants, alleles) |S1 b'T' b'C' b'C' ... b'T' b'A'
    sample_id           (samples) <U3 'S30' 'S32'
    call_genotype       (variants, samples, ploidy) int8 0 -1 0 0 0 ... 1 1 0 0
    call_genotype_mask  (variants, samples, ploidy) bool False True ... False
Attributes:
    contigs:  ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', ...
    source:   sgkit-unknown

Define a new variable based on others?#

Use Xarray’s dictionary like methods, or xarray.Dataset.assign():

In [13]: ds["pos0"] = ds.variant_position - 1

In [14]: ds.assign(pos0 = ds.variant_position - 1) # alternative
Out[14]: 
<xarray.Dataset>
Dimensions:             (contigs: 23, variants: 100, alleles: 2, samples: 50,
                         ploidy: 2)
Dimensions without coordinates: contigs, variants, alleles, samples, ploidy
Data variables:
    contig_id           (contigs) <U2 '0' '1' '2' '3' ... '19' '20' '21' '22'
    variant_contig      (variants) int64 0 0 0 0 0 1 1 ... 21 21 21 22 22 22 22
    variant_position    (variants) int64 0 1 2 3 4 0 1 2 3 ... 3 0 1 2 3 0 1 2 3
    variant_allele      (variants, alleles) |S1 b'T' b'C' b'C' ... b'T' b'A'
    sample_id           (samples) <U3 'S0' 'S1' 'S2' 'S3' ... 'S47' 'S48' 'S49'
    call_genotype       (variants, samples, ploidy) int8 0 0 1 0 1 ... 1 0 1 0 0
    call_genotype_mask  (variants, samples, ploidy) bool False False ... False
    pos0                (variants) int64 -1 0 1 2 3 -1 0 1 ... -1 0 1 2 -1 0 1 2
Attributes:
    contigs:  ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', ...
    source:   sgkit-unknown

Get summary stats?#

Call sample_stats() or variant_stats() as appropriate:

In [15]: sg.sample_stats(ds)
Out[15]: 
<xarray.Dataset>
Dimensions:             (samples: 50, contigs: 23, variants: 100, alleles: 2,
                         ploidy: 2)
Dimensions without coordinates: samples, contigs, variants, alleles, ploidy
Data variables: (12/14)
    sample_n_called     (samples) int64 dask.array<chunksize=(50,), meta=np.ndarray>
    sample_call_rate    (samples) float64 dask.array<chunksize=(50,), meta=np.ndarray>
    sample_n_het        (samples) int64 dask.array<chunksize=(50,), meta=np.ndarray>
    sample_n_hom_ref    (samples) int64 dask.array<chunksize=(50,), meta=np.ndarray>
    sample_n_hom_alt    (samples) int64 dask.array<chunksize=(50,), meta=np.ndarray>
    sample_n_non_ref    (samples) int64 dask.array<chunksize=(50,), meta=np.ndarray>
    ...                  ...
    variant_position    (variants) int64 0 1 2 3 4 0 1 2 3 ... 3 0 1 2 3 0 1 2 3
    variant_allele      (variants, alleles) |S1 b'T' b'C' b'C' ... b'T' b'A'
    sample_id           (samples) <U3 'S0' 'S1' 'S2' 'S3' ... 'S47' 'S48' 'S49'
    call_genotype       (variants, samples, ploidy) int8 0 0 1 0 1 ... 1 0 1 0 0
    call_genotype_mask  (variants, samples, ploidy) bool False False ... False
    pos0                (variants) int64 -1 0 1 2 3 -1 0 1 ... -1 0 1 2 -1 0 1 2
Attributes:
    contigs:  ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', ...
    source:   sgkit-unknown

In [16]: sg.variant_stats(ds)
Out[16]: 
<xarray.Dataset>
Dimensions:                   (variants: 100, alleles: 2, contigs: 23,
                               samples: 50, ploidy: 2)
Dimensions without coordinates: variants, alleles, contigs, samples, ploidy
Data variables: (12/17)
    variant_n_called          (variants) int64 dask.array<chunksize=(100,), meta=np.ndarray>
    variant_call_rate         (variants) float64 dask.array<chunksize=(100,), meta=np.ndarray>
    variant_n_het             (variants) int64 dask.array<chunksize=(100,), meta=np.ndarray>
    variant_n_hom_ref         (variants) int64 dask.array<chunksize=(100,), meta=np.ndarray>
    variant_n_hom_alt         (variants) int64 dask.array<chunksize=(100,), meta=np.ndarray>
    variant_n_non_ref         (variants) int64 dask.array<chunksize=(100,), meta=np.ndarray>
    ...                        ...
    variant_position          (variants) int64 0 1 2 3 4 0 1 2 ... 1 2 3 0 1 2 3
    variant_allele            (variants, alleles) |S1 b'T' b'C' ... b'T' b'A'
    sample_id                 (samples) <U3 'S0' 'S1' 'S2' ... 'S47' 'S48' 'S49'
    call_genotype             (variants, samples, ploidy) int8 0 0 1 0 ... 1 0 0
    call_genotype_mask        (variants, samples, ploidy) bool False ... False
    pos0                      (variants) int64 -1 0 1 2 3 -1 0 ... 1 2 -1 0 1 2
Attributes:
    contigs:  ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', ...
    source:   sgkit-unknown

Filter variants?#

Call xarray.Dataset.sel() on the variants dimension:

In [17]: ds2 = sg.hardy_weinberg_test(ds)

In [18]: ds2.sel(variants=(ds2.variant_hwe_p_value > 1e-2).compute())
Out[18]: 
<xarray.Dataset>
Dimensions:                 (variants: 99, genotypes: 3, contigs: 23,
                             alleles: 2, samples: 50, ploidy: 2)
Coordinates:
  * genotypes               (genotypes) <U3 '0/0' '0/1' '1/1'
Dimensions without coordinates: variants, contigs, alleles, samples, ploidy
Data variables:
    variant_hwe_p_value     (variants) float64 dask.array<chunksize=(99,), meta=np.ndarray>
    variant_genotype_count  (variants, genotypes) uint64 dask.array<chunksize=(99, 3), meta=np.ndarray>
    genotype_id             (genotypes) <U3 dask.array<chunksize=(3,), meta=np.ndarray>
    contig_id               (contigs) <U2 '0' '1' '2' '3' ... '20' '21' '22'
    variant_contig          (variants) int64 0 0 0 0 1 1 1 ... 21 21 22 22 22 22
    variant_position        (variants) int64 0 2 3 4 0 1 2 3 ... 0 1 2 3 0 1 2 3
    variant_allele          (variants, alleles) |S1 b'T' b'C' b'T' ... b'T' b'A'
    sample_id               (samples) <U3 'S0' 'S1' 'S2' ... 'S47' 'S48' 'S49'
    call_genotype           (variants, samples, ploidy) int8 0 0 1 0 ... 0 1 0 0
    call_genotype_mask      (variants, samples, ploidy) bool False ... False
    pos0                    (variants) int64 -1 1 2 3 -1 0 1 ... 0 1 2 -1 0 1 2
Attributes:
    contigs:  ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', ...
    source:   sgkit-unknown

Note

The call to compute is needed to avoid an Xarray error.

Find which new variables were added by a method?#

Use xarray.Dataset.data_vars to compare the new dataset variables to the old:

In [19]: ds2 = sg.sample_stats(ds)

In [20]: set(ds2.data_vars) - set(ds.data_vars)
Out[20]: 
{'sample_call_rate',
 'sample_n_called',
 'sample_n_het',
 'sample_n_hom_alt',
 'sample_n_hom_ref',
 'sample_n_non_ref'}

Save results to a Zarr file?#

Call save_dataset():

In [21]: sg.save_dataset(ds, "ds.zarr")

Note

Zarr datasets must have equal-sized chunks (except for the final chunk, which may be smaller), so you may have to rechunk the dataset first.

Load a dataset from Zarr?#

Call load_dataset():

In [22]: ds = sg.load_dataset("ds.zarr")