The free version of GraphDB can be downloaded and installed from the website.

octofludb can be installed through pip:

You may need to replace pip with pip3, depending on your setup.

Documentation for all subcommands in octofludb can be accessed by calling octofludb -h. Subcommand documentation can be accessed with octofludb <subcommand> -h.

Once GraphDB is installed on your system, you can initialize an empty repo with octofludb init.

A local database can be accessed by submitting SPARQL queries. A SPARQL script specifies a pattern over a graph and when given to the database will print a row of data for every occurrence of the pattern.

Below is a very simple script that returns the strain names of every gamma:

The PREFIX line defines an abbreviatin, f, for the long URI https://flu-crew.org/term/. While the URI looks like a web link, it does not need to point to a real web page (although it may, and that is a good place to put documentation). Every edge type defined specifically for OctofluDB has this prefix.

The SELECT term is followed by a list of all the columns we want printed out for us, in this case just one.

?sid f:strain_name ?name . is a graph pattern restricting the returned data to things related to eachother by the f:strain_name predicate. Words that begin with ? are variables. There is nothing special about ?sid and ?name, these variables could be named anything. By convention, I normally use ?sid for a *s*train *id*entifier and ?gid for a *g*enbank *id*enditifer. Every statement must end with a ., otherwise GraphDB will spill out a very ugly error statement.

?sid f:has_segment ?gid . links our strain node to its child segment node. This particular link will appear in many queries.

?gid f:clade "gamma" links the segment to the literal string value "gamma" by the f:clade relation.

The query can be sent to the database with the command:

The extension .rq is convention for SPARQL scripts.

I will describe how to delete data with a few examples.

There are three basic steps to uploading data to the database.

The first is to prepare the raw data in the "proper" way.

The second step is to convert the raw data to the turtle-format GraphDB understands. This is done through the mk_* commands in octofludb. These commands do more than just translate the data, they use specialized parsers to extract flu specific information from the data. For example, if

The third step is to upload the turtle files to the database. If you are working with a local database, this step is easy: octofludb upload <turtle-fiels>.

Recipe functions are data synthesis functions that are built into octofludb. They will be added on an as-needed database.

You have just retrieved H1 pandemic sequence from the last few years and want to do something awesome with it. With smof, you can see how many sequences you have:

Oh crap, one of our sequences has a length of 1.

The HAs were trimmed to the HA1 region, so these sequences are simply truncated sequences. We can either filter them out or change our methods to not trim. I’ll go with the third option of leaving them in and building the tree (since this is the tree section).

This will take awhile to run, since there are nearly 27 thousand strains. The result tree is shown below.

We want to root it on California/07/2009 and then sort it. This can be done with the treebender utility from the phylommand package (https://github.com/RybergGroup/phylommand):

The outgroup name needs to be an exact match against the tip label. The result is much cleaner:

But this tree is far too large for our purposes. We need to subsample. However, since we are interested in human-swine transmission, we need to preserve the overall topology of the tree with respect to the two species. If we just sampled 5% of the tips, we would loose most of the non-persisting spill-over cases. smot can sample for us while preserving the topology.

This yields a much more manageable tree with only 2643 tips. The black and orange represent human and swine sequences, respectively.

If you just want to see the swine, smot can grep them out:

Or if you want to remove the private VDL sequences:

Publicly available HA segment sequences are inconsistent with regards to what part of the gene the include. Some include DNA before the start codon, some begin at the start codon, some begin after the start codon. This variety introduces non-biological noise into phylogenetic studies. Furthermore, the HA protein is processed after translation to remove the initial signal peptide and to cut the protein into two subunits, HA1 and HA2. For phylogenetic work, we usually focus on HA1. For more consistent phylogenetic results, it is good practice to trim the HA DNA sequences to the part that codes for HA1.

flutile trim can extract HA1 regions for H1 and H3 sequences (handling for the other HA subtypes could be easily added given a reference strain and motifs). flutile finds the HA1 by aligning the input sequences against the reference sequence (California/07/2009) and then extracting from the alignment the regions that map to the reference HA1. This will correctly find HA1 regions even if they are truncated. Here are a few examples:

For tree building, the --conversion=dna2dna option is strongly recommended.

To build a gamma tree from the octoFLU reference file, for example, we might use the following command:

The realignment is needed since the alignment flutile makes internally is of translated inputs against the protein reference.

flutile aadiff takes a multiple-sequence alignment as an input and creates a character difference table. This command is designed for preparing amino acid difference tables. Below is an example of a comparison of 4 H1 sequences.

flutile aadiff --subtype=H1 mufile.faa

The --subtype=H1 argument tells flutile to align the inputs against an H1 reference (A/United Kingdom/1/1933). The reference is used to determine relative indices (the sites column). The index reference is used only for indexing and does not appear in the final table. The first three rows (sites -3, -2, -1) align to the three residues at the end of the signal peptide. Site 1 is the first residue in the mature peptide. Any gaps in the reference alignment are indexed as <ref_id>+<offset>, for example 154+1 and 154+2 are positions 1 and 2 residues after the reference position 154. flutile uses the references from [2] (summarized in table below).

The reference strains selected by [2] may be used for indexing HA, but no comparable study has been done for NA or the internal genes.

Indels are rather rare in the internal genes (with the interesting exception of NS), so the genes are covered well with a single reference. I chose the H1N1 A/Puerto_Rico/8/1934 (EPI_ISL_146807) for this role since it 1) has a complete genome sequence with no ambiguous residues and 2) has been the subject of extensive structural research.

I retrieved all available public WGS data including around 70000 to 80000 sequences for each internal gene. For each segment (other than NS) the 1st, 2nd and 3rd quantiles were all equal (indicating at least 50% of the strains in each group are of exactly the same length). The lengths of the translated A/Puerto_Rico/8/1934 segments exactly matched these quantiles. For NS, the translated product matched the 2nd quantile value (the median).

So overall, while A/Puerto_Rico/8/1934 may be very different in sequence from modern strains, it will work fine as an index reference.

The NA’s are a bit more complicated. The NA length varies dramatically even within subtypes, particularly with variation in length of the stalk-motif [5, 6].

flutile can also use these subtype reference to extract HA motifs relative to Burke2014 numbering.

For example, to exact the classical H1 antigenic motifs: