You can run this tutorial in two ways:
However, since TileDB Cloud has a free tier, we strongly recommend that you sign up and run everything there, as that requires no installations or deployment.
The TileDB-VCF API ingests variant data for rapid access into a 3D array. One of these three dimensions corresponds to the sample_name of the variant data. This enables consumers to query their dataset by one or more samples, extracting attributes for that sample.
Sometimes, for more complex experiments, it may be necessary to generate sample cohorts to use for querying a tiledbvcf.Dataset as a prerequisite step. For example, a researcher may desire to pull all variant data from a tiledbvcf.Dataset related to a particular patient ID, where a patient may have a 1-to-1 or 1-to-many relationship with samples. Because that information isn’t stored in VCF files and, as a result, is not included in a tiledbvcf.Dataset, you can use an accompanying array relating patient ID to sample_name. With this setup, you would first query the array for all samples belonging to a patient and then those samples used to query the tiledbvcf.Dataset.
In this tutorial, you’ll learn how to extend TileDB-VCF with the help of TileDB arrays to perform complex queries.
import os
import shutil
import tiledb
import tiledbvcf
import pandas as pd
print("TileDB core version: {}".format(tiledb.libtiledb.version()))
print("TileDB-Py version: {}".format(tiledb.version()))
print("TileDB-VCF version: {}".format(tiledbvcf.version))TileDB core version: (2, 24, 2)
TileDB-Py version: (0, 30, 2)
TileDB-VCF version: 0.30.0rc0First, pull a public TileDB-VCF dataset and choose some random samples for this tutorial. The TileDB team has made available several public TileDB-VCF datasets.
This tutorial uses the Chr15 canine dataset to explore canine variant data.
Now that you have a tiledbvcf.Dataset and some example sample names, simulate some corresponding metadata. This metadata will include other observations related to those samples such as parent, weight, and whether the sample is a pure breed. First create a pandas.DataFrame and then convert that to an array.
metadata_df = pd.DataFrame(
{
"Sample": query_samples,
"Parent": ["X400", "X400", "X301", "X203"],
"Weight": [45, 23, 11, 57],
"PureBreed": [True, True, False, True],
},
)
metadata_uri = "local_metadata_array"
def convert():
"""Convert pandas.DataFrame to TileDB array.
Delete array if it exists.
"""
try:
tiledb.from_pandas(
metadata_uri,
metadata_df,
index_dims=["Sample"],
)
except tiledb.TileDBError:
print("Array already exists, refreshing.")
if os.path.exists(metadata_uri):
shutil.rmtree(metadata_uri)
convert()
convert()Verifying the schema looks as expected.
Note that during the array transformation, TileDB established the data types in the metadata array schema because the original types were not explicitly defined by the caller. Most data types should look familiar, with the exception of the Parent attribute, which is of type <U0. This represents a variable-length Unicode string in numpy, as described here.
ArraySchema(
domain=Domain(*[
Dim(name='Sample', domain=('', ''), tile=None, dtype='|S0', var=True, filters=FilterList([ZstdFilter(level=-1), ])),
]),
attrs=[
Attr(name='Parent', dtype='<U0', var=True, nullable=False, enum_label=None, filters=FilterList([ZstdFilter(level=-1), ])),
Attr(name='Weight', dtype='int64', var=False, nullable=False, enum_label=None, filters=FilterList([ZstdFilter(level=-1), ])),
Attr(name='PureBreed', dtype='bool', var=False, nullable=False, enum_label=None, filters=FilterList([ZstdFilter(level=-1), ])),
],
cell_order='row-major',
tile_order='row-major',
capacity=10000,
sparse=True,
allows_duplicates=True,
)
Starting with a given Parent, query the metadata array to extract sample names to subsequently query the tiledbvcf.Dataset.
# example parent we'd like to get samples for
query_parent = "X400"
with tiledb.open(metadata_uri, "r") as Ar:
# query metadata array for
sample_result = Ar.query(cond=f"Parent == '{query_parent}'")[:]
samples = [s.decode() for s in sample_result["Sample"]]
print(samples)['YorkshireTerrier73', 'YorkshireTerrier75']variant_results = ds.read(
regions=["chr15:1-5000"],
samples=samples,
attrs=[
"sample_name",
"contig",
"pos_start",
"pos_end",
"alleles",
"fmt_GT",
],
)
print(variant_results) sample_name contig pos_start pos_end alleles fmt_GT
0 YorkshireTerrier73 chr15 177 177 [T, G] [-1, -1]
1 YorkshireTerrier75 chr15 177 177 [T, G] [-1, -1]
2 YorkshireTerrier73 chr15 213 213 [C, A] [0, 0]
3 YorkshireTerrier75 chr15 213 213 [C, A] [0, 0]
4 YorkshireTerrier73 chr15 231 231 [C, T] [0, 0]
.. ... ... ... ... ... ...
983 YorkshireTerrier75 chr15 4980 4980 [C, T] [0, 0]
984 YorkshireTerrier73 chr15 4985 4986 [CG, C] [0, 0]
985 YorkshireTerrier75 chr15 4985 4986 [CG, C] [0, 0]
986 YorkshireTerrier73 chr15 4994 4995 [CT, C] [0, 0]
987 YorkshireTerrier75 chr15 4994 4995 [CT, C] [-1, -1]
[988 rows x 6 columns]As demonstrated here, storing metadata in a separate array is a useful strategy for querying a tiledbvcf.Dataset using supporting data to the TileDB-VCF API.