Write basic queries
TODO
In this tutorial, we will write basic queries to work with metadata and alignments in RhyDB. It is recommended to first read this short introduction into the query language but it is not strictly necessary.
Loading the data
We will use a small dataset consisting a single table with 1,000 SARS-CoV-2 sequences. We will use the WebAssembly version of RhyDB and run everything within your browser. Click on the following button to index and load the sequences into the browser. You can also download the data file here if you would like to inspect it (note: the file is compressed with Zstadard and you may need to install a software to decompress it).
The dataset has a single table called “default”. Let’s inspect it:
Load the whole table by running the following query:
defaultSee the schema of the table:
default.schema()The table describes where and when each sample was taken: region, country and division hold the sampling location and date holds the sampling date. How a sequence was classified is recorded in pango_lineage, nextstrain_clade and who_clade. The sequences themselves are in unaligned_main, which holds the nucleotide sequence as submitted, in main, which holds the same sequence aligned against the reference genome, and in one field per gene, such as S, ORF1a and E, holding the aligned amino acid sequences.
Important: The query language is case-sensitive both in terms of the names of the tables and fields and in terms of the operators.
Filter and sort metadata
To get familiar with the basic operators of RhyDB, we will first run a few queries on the metadata.
Use project() to define the columns that should be returned. The operator accepts a set of field names.
default.project({ strain, country, pango_lineage })Use filter() to filter the rows — here, all sequences from Switzerland with a coverage of at least 95%:
default.filter(country = 'Switzerland' && coverage >= 0.95)Use orderBy() to sort the data and limit() to restrict the number of returned rows. orderBy() accepts a set of
fields, as it is possible to sort by several of them. Fields are sorted ascending by default; add .desc() to a
field name for a descending order. This returns the 20 most recent entries:
default.orderBy({ date.desc() }).limit(20)Now combine the operators yourself: get the 10 earliest sequences from Europe outside of Switzerland and look at the
fields date, country and pango_lineage. From which country is the earliest sequence and which lineage did it
have? (Hint: the operator for unequal is <>, and .isNotNull() filters out missing data.)
The United Kingdom, with lineage B.1.
Grouping and aggregations
Instead of returning individual sequences, groupBy() aggregates them. Its first argument is a record of aggregates, currently only count(), and its optional second argument is the set of fields to group by.
Without grouping fields, the whole table is aggregated into a single row.
default.groupBy({ count := count() })Add a set of grouping fields to get one row per value.
default.groupBy({ count := count() }, { region })Group by several fields and sort the result to get the five most frequent combinations of country and lineage.
default.groupBy({ count := count() }, { country, pango_lineage }).orderBy({ count.desc() }).limit(5)Combine the operators yourself: which country contributed the most Delta sequences, and how many? (Hint: filter on
who_clade first, then group by country.)
The United Kingdom, with 123 sequences.
Filtering and grouping mutations (TODO: refinement needed after the new S[123] syntax has been implemented)
main is aligned against the reference genome and the gene fields against the reference proteins, so a position means the same thing in every sequence. nucleotideEquals(position := n, symbol := s, sequenceName := 'main') keeps the sequences carrying symbol s at the 1-based reference position n, and aminoAcidEquals(position := n, symbol := s, sequenceName := 'S') does the same within a gene. hasMutation() and hasAAMutation() take only a position and a sequence name, and keep the sequences that differ from the reference there, whatever the symbol.
Count the sequences carrying the S:D614G substitution — nearly the whole dataset.
default.filter(aminoAcidEquals(position := 614, symbol := 'G', sequenceName := 'S')).groupBy({ count := count() })The same substitution can be queried on the nucleotide level, where it is A23403G. It matches a few more sequences than the amino acid query, because a nucleotide can be called even when the surrounding codon cannot.
default.filter(nucleotideEquals(position := 23403, symbol := 'G', sequenceName := 'main')).groupBy({ count := count() })Position 484 of the spike protein carries different substitutions in different clades. Count every sequence that
differs from the reference there, then count only those carrying alanine with aminoAcidEquals(position := 484, symbol := 'A', sequenceName := 'S').
default.filter(hasAAMutation(position := 484, sequenceName := 'S')).groupBy({ count := count() })532 sequences differ from the reference, 489 of them carry A.
Combine mutations with && and group the matches. Which clades combine both substitutions?
default.filter(aminoAcidEquals(position := 501, symbol := 'Y', sequenceName := 'S') && aminoAcidEquals(position := 484, symbol := 'K', sequenceName := 'S')).groupBy({ count := count() }, { who_clade }).orderBy({ count.desc() })Omicron, Gamma, Mu and Beta, plus a few sequences without a clade assignment.
Rather than asking about one symbol at a time, you can also read the symbol itself: at() returns the symbol a sequence carries at a reference position, and map() turns that expression into a named column. groupBy() groups by existing columns only, so map the positions first and then group by them — the result is a co-occurrence table of symbol combinations and how often each occurs. Positions that could not be called appear as X in the genes and as N in main, and deletions appear as -.
Look at a single position: how many sequences carry each amino acid at position 501 of the spike protein?
default.map({ s501 := S.at(501) }).groupBy({ count := count() }, { s501 }).orderBy({ count.desc() })Map a second position to get a co-occurrence table. The most frequent combination is 501Y with 484A.
default.map({ s501 := S.at(501), s484 := S.at(484) }).groupBy({ count := count() }, { s501, s484 }).orderBy({ count.desc() }).limit(5)Restrict a co-occurrence table to a clade: build one over positions 339 and 371 of S for the sequences with
who_clade = 'Omicron'. Then group the sequences carrying L at position 371 by nextstrain_clade. Which
Nextstrain clades are they?
21K only — the earliest Omicron clade. The later ones carry F.