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-0.6.1.dev2+gcc728043
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-0.6.1.dev2+gcc728043
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-0.6.1.dev2+gcc728043
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-0.6.1.dev2+gcc728043
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-0.6.1.dev2+gcc728043
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-0.6.1.dev2+gcc728043
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 81 78 80 77 88 82 ... 76 80 77 90 81 81
sample_call_rate (samples) float64 0.81 0.78 0.8 0.77 ... 0.9 0.81 0.81
sample_n_het (samples) int64 40 37 40 39 41 33 ... 47 37 40 49 50 37
sample_n_hom_ref (samples) int64 20 22 19 21 21 23 ... 9 19 20 17 19 24
sample_n_hom_alt (samples) int64 21 19 21 17 26 26 ... 20 24 17 24 12 20
sample_n_non_ref (samples) int64 61 56 61 56 67 59 ... 67 61 57 73 62 57
... ...
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-0.6.1.dev2+gcc728043
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 40 42 43 40 41 ... 42 38 43 43 40
variant_call_rate (variants) float64 0.8 0.84 0.86 ... 0.86 0.86 0.8
variant_n_het (variants) int64 22 30 18 26 29 ... 23 22 26 21 19
variant_n_hom_ref (variants) int64 11 3 13 3 6 8 ... 6 10 9 12 8 13
variant_n_hom_alt (variants) int64 7 9 12 11 6 17 ... 15 9 7 5 14 8
variant_n_non_ref (variants) int64 29 39 30 37 35 ... 32 29 31 35 27
... ...
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-0.6.1.dev2+gcc728043
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))
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-0.6.1.dev2+gcc728043
Note
Filtering causes an eager Xarray computation.
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")