@@ -8,7 +8,6 @@ sink(log, type = "output")
88library(tidyverse )
99library(jsonlite )
1010library(ape )
11- library(ggtree )
1211library(logger )
1312
1413log_threshold(INFO )
@@ -23,40 +22,44 @@ TREE_LEGEND_NAMES["boot_alrt_pass"] <- sprintf(
2322 snakemake @ params [[" alrt_th" ]], " %"
2423)
2524
25+ log_info(" Reading tree" )
2626tree_ml <- read.tree(snakemake @ input [[" tree" ]]) %> %
27- root(
28- snakemake @ params [[" ref_name" ]],
29- resolve.root = TRUE
30- )
27+ root(
28+ snakemake @ params [[" ref_name" ]],
29+ resolve.root = TRUE
30+ )
31+ log_debug(" Read tree with {length(tree_ml$node.label)} labels" )
3132
33+ log_info(" Reading target names from FASTA" )
3234target_names <- read.dna(
33- snakemake @ input [[" target_fasta" ]],
34- format = " fasta" ,
35- as.matrix = FALSE ,
35+ snakemake @ input [[" target_fasta" ]],
36+ format = " fasta" ,
37+ as.matrix = FALSE ,
3638)
39+ log_debug(" Read {length(target_names)} records" )
3740
41+ log_info(" Processing target names" )
3842target_names <- target_names [
39- ! startsWith(names(target_names ), snakemake @ config [[" ALIGNMENT_REFERENCE" ]])
43+ ! startsWith(names(target_names ), snakemake @ config [[" ALIGNMENT_REFERENCE" ]])
4044]
4145target_names <- names(target_names )
46+ log_debug(" {length(target_names)} records remaining after processing" )
4247
4348# ML tree with context data
4449# Internal nodes color
4550# Node labels contain SH-aLRT/UFboot values
51+ log_info(" Reading support values from labels" )
52+ labels <- strsplit(tree_ml $ node.label , " /" )
4653aLRT.values <- sapply(
47- strsplit(tree_ml $ node.label , " /" ),
48- function (x ) {
49- as.numeric(x [1 ])
50- }
54+ labels ,
55+ function (x ) as.numeric(x [1 ])
5156)
52-
5357bootstrap.values <- sapply(
54- strsplit(tree_ml $ node.label , " /" ),
55- function (x ) {
56- as.numeric(x [2 ])
57- }
58+ labels ,
59+ function (x ) as.numeric(x [2 ])
5860)
5961
62+ log_info(" Calculating support mask for the given thresholds" )
6063aLRT.mask <- aLRT.values > = snakemake @ params [[" alrt_th" ]]
6164boot.mask <- bootstrap.values > = snakemake @ params [[" boot_th" ]]
6265
@@ -73,15 +76,15 @@ tree_ml.labels <- tree_ml$tip.label[1:tree_ml.nodes]
7376tree_ml.node.pass <- c(rep(FALSE , tree_ml.ntips ), aLRT.mask & boot.mask )
7477
7578ml.tree.annot <- tibble(
76- node = 1 : tree_ml.nodes ,
79+ node = 1 : tree_ml.nodes ,
7780) %> %
78- mutate(
79- Class = case_when(
80- tree_ml.labels %in% target_names ~ TREE_LEGEND_NAMES [" tip_label" ],
81- tree_ml.node.pass ~ TREE_LEGEND_NAMES [" boot_alrt_pass" ],
82- TRUE ~ NA
83- )
81+ mutate(
82+ Class = case_when(
83+ tree_ml.labels %in% target_names ~ TREE_LEGEND_NAMES [" tip_label" ],
84+ tree_ml.node.pass ~ TREE_LEGEND_NAMES [" boot_alrt_pass" ],
85+ TRUE ~ NA
8486 )
87+ )
8588
8689# Write output files
8790log_info(" Writing tree annotation" )
@@ -90,22 +93,22 @@ write_csv(ml.tree.annot, snakemake@output$annotation)
9093log_info(" Writing JSON data" )
9194target.node <- tree_ml $ node.label [target.mrca - length(tree_ml $ tip.label )]
9295list (
93- " boot" = strsplit(target.node , " /" )[[1 ]][2 ] %> % as.numeric(),
94- " alrt" = strsplit(target.node , " /" )[[1 ]][1 ] %> % as.numeric(),
95- " monophyly" = ifelse(
96- is.monophyletic(tree_ml , target_names ),
97- " are" ,
98- " are not"
99- ),
100- " target_mrca" = target.mrca ,
101- " clade_tips" = target.mrca.clade.ntips ,
102- " max_tip_length" = max(node.depth.edgelength(tree_ml )[
103- 1 : length(tree_ml $ tip.label )
104- ]),
105- " root" = snakemake @ params [[" ref_name" ]]
96+ " boot" = strsplit(target.node , " /" )[[1 ]][2 ] %> % as.numeric(),
97+ " alrt" = strsplit(target.node , " /" )[[1 ]][1 ] %> % as.numeric(),
98+ " monophyly" = ifelse(
99+ is.monophyletic(tree_ml , target_names ),
100+ " are" ,
101+ " are not"
102+ ),
103+ " target_mrca" = target.mrca ,
104+ " clade_tips" = target.mrca.clade.ntips ,
105+ " max_tip_length" = max(node.depth.edgelength(tree_ml )[
106+ 1 : length(tree_ml $ tip.label )
107+ ]),
108+ " root" = snakemake @ params [[" ref_name" ]]
106109) %> %
107- write_json(
108- snakemake @ output $ json ,
109- auto_unbox = TRUE ,
110- digits = NA
111- )
110+ write_json(
111+ snakemake @ output $ json ,
112+ auto_unbox = TRUE ,
113+ digits = NA
114+ )
0 commit comments