Skip to content

Commit e7b07da

Browse files
anna-parkeranna-parker
committed
feat(silo): read in phylogenetic tree files and create tree structure (#806)
* feat(silo): parse auspice.json * feat(lapis): add newick parser and tests * feat(silo): add tests for edge cases and error handling * feat(silo): switch to using pointers * feat(silo): try to create TreeNodeId type * feat(silo): add to e2e tests * feat(silo): add link between metadata field and tree in initializer * feat(silo): refactor reading from file code * feat(silo): add docs and improve variable names * feat(silo): update names * feat(silo): add newick file to e2e tests * feat(silo): move phylo_tree class to common * feat(silo): add better error type * feat(silo): improve function structure
1 parent 3f8f2ac commit e7b07da

28 files changed

Lines changed: 635 additions & 18 deletions

‎documentation/phylogenetic_queries.md‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
# Documentation: Phylogenetic Queries
2+
3+
By default SILO does not support phylogeny-based queries. In order to enable these queries a file containing the phylogeny of all (or a subset) of sequences as leaf nodes needs to be imported during pre-processing. The tree file can be specified with the `lineageDefinitionsFilename` field of the preprocessing config. Additionally, to enable querying users must specify which metadata column corresponds to the nodes in the tree file. This should be done by adding the `phyloTreeNodeIdentifier` to the respective metadata field in the database config.
4+
5+
> **Warning**
6+
>
7+
> Like the `generateLineageIndex` option the `phyloTreeNodeIdentifier` option can only be applied to metadata fields that are strings. However, unlike the `generateLineageIndex` option the `generateIndex` should NOT be set. (Internally lineage indexes are queried using a bitmap structure, however for phylogenetic queries we use the tree structure directly).
8+
9+
## Tree File Format
10+
11+
We allow users to specify a tree phylogeny using two standard formats: [newick](https://en.wikipedia.org/wiki/Newick_format) and [auspice.json (v2)](https://docs.nextstrain.org/projects/auspice/en/stable/releases/v2.html#new-dataset-json-format). We do not currently support ancestral reassortment graphs. Unlike the standard format we additionally require that **all nodes (internal and leaves) of the tree should be uniquely labelled**.
12+
13+
// TODO: Add details of the phylogenetic queries that SILO now supports
14+

‎src/config/source/yaml_file.test.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -58,6 +58,7 @@ inputDirectory: "./testBaseData/exampleDataset/"
5858
outputDirectory: "./output/"
5959
ndjsonInputFilename: "input_file.ndjson"
6060
lineageDefinitionsFilename: "lineage_definitions.yaml"
61+
phyloTreeFilename: "phylogenetic_tree.yaml"
6162
referenceGenomeFilename: "reference_genomes.json"
6263
)"
6364
)
@@ -70,6 +71,7 @@ referenceGenomeFilename: "reference_genomes.json"
7071
{YamlFile::stringToConfigKeyPath("ndjsonInputFilename"), YAML::Node{"input_file.ndjson"}},
7172
{YamlFile::stringToConfigKeyPath("lineageDefinitionsFilename"),
7273
YAML::Node{"lineage_definitions.yaml"}},
74+
{YamlFile::stringToConfigKeyPath("phyloTreeFilename"), YAML::Node{"phylogenetic_tree.yaml"}},
7375
{YamlFile::stringToConfigKeyPath("referenceGenomeFilename"),
7476
YAML::Node{"reference_genomes.json"}},
7577
};

‎src/silo/common/phylo_tree.cpp‎

Lines changed: 230 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,230 @@
1+
#include "silo/common/phylo_tree.h"
2+
3+
#include <fstream>
4+
#include <sstream>
5+
6+
#include <nlohmann/json.hpp>
7+
8+
#include "silo/preprocessing/preprocessing_exception.h"
9+
#include "silo/query_engine/batched_bitmap_reader.h"
10+
11+
namespace silo::common {
12+
using silo::common::TreeNodeId;
13+
14+
std::shared_ptr<TreeNode> parse_auspice_tree(
15+
const nlohmann::json& j,
16+
std::optional<std::shared_ptr<TreeNode>> parent,
17+
std::unordered_map<TreeNodeId, std::shared_ptr<TreeNode>>& node_map,
18+
int depth = 0
19+
) {
20+
auto node = std::make_shared<TreeNode>();
21+
if (!j.contains("name")) {
22+
throw silo::preprocessing::PreprocessingException(
23+
"Invalid File: Auspice JSON node does not contain a 'name' entry."
24+
);
25+
}
26+
node->node_id = TreeNodeId{j.at("name").get<std::string>()};
27+
node->parent = parent;
28+
node->depth = depth;
29+
30+
const auto& children = j.contains("children") ? j["children"] : nlohmann::json::array();
31+
32+
for (const auto& child : children) {
33+
auto child_node = parse_auspice_tree(child, node, node_map, depth + 1);
34+
node->children.push_back(child_node);
35+
}
36+
37+
// Insert node into the map *after* children so it's fully constructed
38+
node_map[node->node_id] = node;
39+
return node;
40+
}
41+
42+
PhyloTree PhyloTree::fromAuspiceJSONString(const std::string& json_string) {
43+
nlohmann::json j = nlohmann::json::parse(json_string);
44+
45+
if (!j.contains("tree")) {
46+
throw silo::preprocessing::PreprocessingException(
47+
"Invalid File: Auspice JSON does not contain a 'tree' entry."
48+
);
49+
}
50+
51+
PhyloTree file;
52+
auto root = parse_auspice_tree(j["tree"], std::nullopt, file.nodes);
53+
return file;
54+
}
55+
56+
PhyloTree PhyloTree::fromAuspiceJSONFile(const std::filesystem::path& json_path) {
57+
std::ifstream file(json_path, std::ios::in | std::ios::binary);
58+
if (!file) {
59+
throw silo::preprocessing::PreprocessingException(
60+
fmt::format("Could not open the JSON file: '{}'", json_path.string())
61+
);
62+
}
63+
64+
std::ostringstream contents;
65+
if (file.peek() != std::ifstream::traits_type::eof()) {
66+
contents << file.rdbuf();
67+
if (contents.fail()) {
68+
throw silo::preprocessing::PreprocessingException(
69+
fmt::format("Error when reading the JSON file: '{}'", json_path.string())
70+
);
71+
}
72+
}
73+
try {
74+
return fromAuspiceJSONString(contents.str());
75+
} catch (const nlohmann::json::parse_error& parse_exception) {
76+
throw silo::preprocessing::PreprocessingException(
77+
fmt::format("The JSON file '{}' does not contain valid JSON.", json_path.string())
78+
);
79+
}
80+
}
81+
82+
bool isValidLabelChar(char c) {
83+
return isalnum(c) || c == '_' || c == '.' || c == '-';
84+
}
85+
86+
double parseBranchLength(std::string_view& sv) {
87+
std::string number;
88+
while (!sv.empty() && (isdigit(sv.front()) || sv.front() == '.' || sv.front() == '-' ||
89+
sv.front() == '+' || sv.front() == 'e')) {
90+
number += sv.front();
91+
sv.remove_prefix(1);
92+
}
93+
return number.empty() ? 0.0 : std::stod(number);
94+
}
95+
96+
TreeNodeId parseLabel(std::string_view& sv) {
97+
std::string label;
98+
while (!sv.empty() && isValidLabelChar(sv.front())) {
99+
label += sv.front();
100+
sv.remove_prefix(1);
101+
}
102+
if (!sv.empty() && sv.front() == ':') {
103+
sv.remove_prefix(1);
104+
parseBranchLength(sv);
105+
}
106+
return TreeNodeId{label};
107+
}
108+
109+
void skipWhitespace(std::string_view& sv) {
110+
while (!sv.empty() && std::isspace(sv.front())) {
111+
sv.remove_prefix(1);
112+
}
113+
}
114+
115+
std::shared_ptr<TreeNode> parseSubtree(
116+
std::string_view& sv,
117+
std::unordered_map<TreeNodeId, std::shared_ptr<TreeNode>>& node_map,
118+
int depth = 0,
119+
std::optional<std::shared_ptr<TreeNode>> parent = std::nullopt
120+
) {
121+
auto node = std::make_shared<TreeNode>();
122+
node->depth = depth;
123+
node->parent = parent;
124+
125+
skipWhitespace(sv);
126+
if (!sv.empty() && sv.front() == '(') {
127+
sv.remove_prefix(1);
128+
depth++;
129+
do {
130+
auto child_node = parseSubtree(sv, node_map, depth, node);
131+
node->children.push_back(child_node);
132+
skipWhitespace(sv);
133+
if (!sv.empty() && sv.front() == ',') {
134+
sv.remove_prefix(1);
135+
}
136+
} while (!sv.empty() && sv.front() != ')');
137+
if (!sv.empty() && sv.front() == ')') {
138+
sv.remove_prefix(1);
139+
depth--;
140+
}
141+
}
142+
143+
if (depth != node->depth) {
144+
throw silo::preprocessing::PreprocessingException(
145+
"Parenthesis mismatch in Newick string - depth does not match"
146+
);
147+
}
148+
149+
skipWhitespace(sv);
150+
node->node_id = parseLabel(sv);
151+
skipWhitespace(sv);
152+
153+
node_map[node->node_id] = node;
154+
155+
return node;
156+
}
157+
158+
PhyloTree PhyloTree::fromNewickString(const std::string& newick_string) {
159+
PhyloTree file;
160+
161+
std::string_view sv(newick_string);
162+
if (sv.empty()) {
163+
throw silo::preprocessing::PreprocessingException(
164+
"Error when parsing the Newick string - The string is empty"
165+
);
166+
}
167+
if (sv.back() != ';') {
168+
throw silo::preprocessing::PreprocessingException(fmt::format(
169+
"Error when parsing the Newick string: '{}' - string does not end in ';'", newick_string
170+
));
171+
}
172+
sv.remove_suffix(1);
173+
try {
174+
auto root = parseSubtree(sv, file.nodes, 0);
175+
if (!sv.empty()) {
176+
throw silo::preprocessing::PreprocessingException(fmt::format(
177+
"Error when parsing the Newick string: '{}' - extra characters found", newick_string
178+
));
179+
}
180+
} catch (const std::exception& e) {
181+
throw silo::preprocessing::PreprocessingException(
182+
fmt::format("Error when parsing the Newick string: '{}'", newick_string)
183+
);
184+
}
185+
186+
return file;
187+
}
188+
189+
PhyloTree PhyloTree::fromNewickFile(const std::filesystem::path& newick_path) {
190+
std::ifstream file(newick_path, std::ios::in | std::ios::binary);
191+
if (!file) {
192+
throw silo::preprocessing::PreprocessingException(
193+
fmt::format("Could not open the Newick file: '{}'", newick_path.string())
194+
);
195+
}
196+
197+
std::ostringstream contents;
198+
if (file.peek() != std::ifstream::traits_type::eof()) {
199+
contents << file.rdbuf();
200+
if (contents.fail()) {
201+
throw silo::preprocessing::PreprocessingException(
202+
fmt::format("Error when reading the Newick file: '{}'", newick_path.string())
203+
);
204+
}
205+
}
206+
try {
207+
return fromNewickString(contents.str());
208+
} catch (const std::exception& e) {
209+
throw silo::preprocessing::PreprocessingException(
210+
fmt::format("Error when parsing the Newick file: '{}'", newick_path.string())
211+
);
212+
}
213+
}
214+
215+
PhyloTree PhyloTree::fromFile(const std::filesystem::path& path) {
216+
auto ext = path.extension().string();
217+
218+
std::transform(ext.begin(), ext.end(), ext.begin(), ::tolower);
219+
220+
if (ext == ".nwk") {
221+
return common::PhyloTree::fromNewickFile(path);
222+
} else if (ext == ".json") {
223+
return common::PhyloTree::fromAuspiceJSONFile(path);
224+
}
225+
throw silo::preprocessing::PreprocessingException(fmt::format(
226+
"Error when parsing tree file: '{}'. Path must end with .nwk or .json", path.string()
227+
));
228+
}
229+
230+
} // namespace silo::common

‎src/silo/common/phylo_tree.h‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
#pragma once
2+
3+
#include <filesystem>
4+
#include <vector>
5+
6+
#include <nlohmann/json.hpp>
7+
8+
#include "silo/common/tree_node_id.h"
9+
10+
namespace silo::common {
11+
using silo::common::TreeNodeId;
12+
13+
class TreeNode {
14+
public:
15+
TreeNodeId node_id;
16+
std::vector<std::shared_ptr<TreeNode>> children;
17+
std::optional<std::shared_ptr<TreeNode>> parent;
18+
int depth;
19+
};
20+
21+
class PhyloTree {
22+
public:
23+
std::unordered_map<TreeNodeId, std::shared_ptr<TreeNode>> nodes;
24+
25+
static PhyloTree fromAuspiceJSONFile(const std::filesystem::path& json_path);
26+
27+
static PhyloTree fromAuspiceJSONString(const std::string& json_string);
28+
29+
static PhyloTree fromNewickFile(const std::filesystem::path& newick_path);
30+
31+
static PhyloTree fromNewickString(const std::string& newick_string);
32+
33+
static PhyloTree fromFile(const std::filesystem::path& path);
34+
};
35+
36+
} // namespace silo::common

‎src/silo/common/phylo_tree.test.cpp‎

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,97 @@
1+
#include "silo/common/phylo_tree.h"
2+
3+
#include <gmock/gmock.h>
4+
#include <gtest/gtest.h>
5+
#include <yaml-cpp/exceptions.h>
6+
7+
#include "silo/preprocessing/preprocessing_exception.h"
8+
9+
using silo::common::PhyloTree;
10+
using silo::common::TreeNodeId;
11+
12+
TEST(PhyloTree, correctlyParsesFromJSON) {
13+
auto phylo_tree_file = PhyloTree::fromAuspiceJSONString(
14+
R"({
15+
"version": "schema version",
16+
"meta": {},
17+
"tree": {
18+
"name": "ROOT",
19+
"children": [
20+
{
21+
"name": "CHILD",
22+
"children": [
23+
{
24+
"name": "CHILD2"
25+
}
26+
]
27+
}
28+
]
29+
}
30+
})"
31+
);
32+
ASSERT_EQ(phylo_tree_file.nodes.size(), 3);
33+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"ROOT"})->parent, std::nullopt);
34+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"ROOT"})->depth, 0);
35+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"ROOT"})->children.size(), 1);
36+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"CHILD"})->depth, 1);
37+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"CHILD"})->children.size(), 1);
38+
ASSERT_EQ(
39+
phylo_tree_file.nodes.at(TreeNodeId{"CHILD"})->children.at(0)->node_id, TreeNodeId{"CHILD2"}
40+
);
41+
ASSERT_EQ(
42+
phylo_tree_file.nodes.at(TreeNodeId{"CHILD2"})->parent->get()->node_id, TreeNodeId{"CHILD"}
43+
);
44+
}
45+
46+
TEST(PhyloTree, throwsOnInvalidJSON) {
47+
EXPECT_THROW(
48+
PhyloTree::fromAuspiceJSONString("{\"invalid\": \"json\"}"),
49+
silo::preprocessing::PreprocessingException
50+
);
51+
}
52+
53+
TEST(PhyloTree, correctlyParsesFromNewick) {
54+
auto phylo_tree_file = PhyloTree::fromNewickString("((CHILD2)CHILD)ROOT;");
55+
ASSERT_EQ(phylo_tree_file.nodes.size(), 3);
56+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"ROOT"})->parent, std::nullopt);
57+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"ROOT"})->depth, 0);
58+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"ROOT"})->children.size(), 1);
59+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"CHILD"})->depth, 1);
60+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"CHILD"})->children.size(), 1);
61+
ASSERT_EQ(
62+
phylo_tree_file.nodes.at(TreeNodeId{"CHILD"})->children.at(0)->node_id, TreeNodeId{"CHILD2"}
63+
);
64+
ASSERT_EQ(
65+
phylo_tree_file.nodes.at(TreeNodeId{"CHILD2"})->parent->get()->node_id, TreeNodeId{"CHILD"}
66+
);
67+
}
68+
69+
TEST(PhyloTree, correctlyParsesFromNewickWithBranchLengths) {
70+
auto phylo_tree_file =
71+
PhyloTree::fromNewickString("((CHILD2:0.5, CHILD3:1)CHILD:0.1, CHILD4:1.5)ROOT;");
72+
ASSERT_EQ(phylo_tree_file.nodes.size(), 5);
73+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"ROOT"})->parent, std::nullopt);
74+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"ROOT"})->depth, 0);
75+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"ROOT"})->children.size(), 2);
76+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"CHILD"})->depth, 1);
77+
ASSERT_EQ(phylo_tree_file.nodes.at(TreeNodeId{"CHILD"})->children.size(), 2);
78+
ASSERT_EQ(
79+
phylo_tree_file.nodes.at(TreeNodeId{"CHILD"})->children.at(0)->node_id, TreeNodeId{"CHILD2"}
80+
);
81+
ASSERT_EQ(
82+
phylo_tree_file.nodes.at(TreeNodeId{"CHILD2"})->parent->get()->node_id, TreeNodeId{"CHILD"}
83+
);
84+
}
85+
86+
TEST(PhyloTree, throwsOnInvalidNewick) {
87+
EXPECT_THROW(
88+
PhyloTree::fromNewickString("((CHILD2)CHILD;"), silo::preprocessing::PreprocessingException
89+
);
90+
}
91+
92+
TEST(PhyloTree, throwsOnInvalidNewickNoSemicolon) {
93+
EXPECT_THROW(
94+
PhyloTree::fromNewickString("((CHILD2)CHILD)ROOT"),
95+
silo::preprocessing::PreprocessingException
96+
);
97+
}

0 commit comments

Comments
 (0)