From 0627c139b88598c172d16bf2b495d062c9fe11a6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20=C3=81lvarez=20Herrera?= Date: Mon, 9 Feb 2026 19:55:51 +0100 Subject: [PATCH 1/5] feat: use all-sites quality filter for the evolutionary rate calculation --- workflow/rules/distances.smk | 49 +++++++++++++++++++++++++++++++++++- workflow/rules/vaf.smk | 16 +++++++++--- 2 files changed, 61 insertions(+), 4 deletions(-) diff --git a/workflow/rules/distances.smk b/workflow/rules/distances.smk index 11fae38..65fe3d4 100644 --- a/workflow/rules/distances.smk +++ b/workflow/rules/distances.smk @@ -1,3 +1,50 @@ +rule compile_fail_sites_vcf: + params: + header = ("#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"), + filter_text = "mask", + sub_text = "NA", + exc_text = "site_qual" + input: + sites = OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv", + output: + sites = temp(OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf"), + run: + import pandas as pd + sites = ( + pd.read_table(input.sites, sep="\t") + .drop_duplicates(subset=("CHROM", "POS", "REF")) + .rename(columns={"CHROM": "#CHROM"}) + ) + sites["ID"] = "." + sites["ALT"] = "." + sites["QUAL"] = "." + sites["FILTER"] = params.filter_text + sites["INFO"] = f"SUB={params.sub_text};EXC={params.exc_text}" + sites[list(params.header)].to_csv(output.sites, sep="\t", index=False) + + +rule merge_sites: + params: + header = ("#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO") + input: + lambda wildcards: select_problematic_vcf(), + OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf", + output: + sites = temp(OUTDIR / "all_mask_sites.vcf"), + run: + import pandas as pd + ( + pd.concat( + [pd.read_table(path, sep="\t", comment="#", names=params.header) for path in input], + axis="rows", + ignore_index=True + ) + .drop_duplicates(subset=("#CHROM", "POS", "FILTER"), keep="first") + .sort_values(list(params.header)) + .to_csv(output.sites, sep="\t", index=False) + ) + + rule extract_afwdist_variants: conda: "../envs/biopython.yaml" params: @@ -8,7 +55,7 @@ rule extract_afwdist_variants: mask_class = ["mask"], input: variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv", - mask_vcf = lambda wildcards: select_problematic_vcf(), + mask_vcf = OUTDIR / "all_mask_sites.vcf", ancestor = OUTDIR/f"{OUTPUT_NAME}.ancestor.fasta", reference = OUTDIR/"reference.fasta", output: diff --git a/workflow/rules/vaf.smk b/workflow/rules/vaf.smk index c42d580..6cccf20 100644 --- a/workflow/rules/vaf.smk +++ b/workflow/rules/vaf.smk @@ -237,7 +237,8 @@ rule filter_mpileup_all_sites: input: OUTDIR / "all_sites" / "{sample}.query.tsv", output: - temp(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv"), + sites_pass = temp(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv"), + sites_fail = temp(OUTDIR / "all_sites" / "{sample}.fail_sites.tsv"), log: LOGDIR / "filter_mpileup_all_sites" / "{sample}.txt" run: @@ -248,11 +249,13 @@ rule filter_mpileup_all_sites: df["TOTAL_AD"] = df.AD.str.split(",").apply(lambda values: sum(int(n) for n in values)) df["TOTAL_ADF"] = df.ADF.str.split(",").apply(lambda values: sum(int(n) for n in values)) df["TOTAL_ADR"] = df.ADR.str.split(",").apply(lambda values: sum(int(n) for n in values)) - df[ + mask = ( (df.TOTAL_AD >= params.min_total_AD) & (df.TOTAL_ADF >= params.min_total_ADF) & (df.TOTAL_ADR >= params.min_total_ADR) - ].to_csv(output[0], sep="\t", index=False) + ) + df[mask].to_csv(output.sites_pass, sep="\t", index=False) + df[~mask].to_csv(output.sites_fail, sep="\t", index=False) use rule concat_vcf_fields as merge_filtered_mpileup_all_sites with: @@ -262,6 +265,13 @@ use rule concat_vcf_fields as merge_filtered_mpileup_all_sites with: OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv", +use rule concat_vcf_fields as merge_fail_mpileup_all_sites with: + input: + expand(OUTDIR / "all_sites" / "{sample}.fail_sites.tsv", sample=iter_samples()), + output: + OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv", + + rule fill_all_sites: conda: "../envs/renv.yaml" input: From c14600f9f9bd4198cd9da321c23b655dfe9d21a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20=C3=81lvarez=20Herrera?= Date: Tue, 10 Feb 2026 11:49:18 +0100 Subject: [PATCH 2/5] refactor: mark failed sites in merged VCF filter Remove VCF header from params as well --- workflow/rules/distances.smk | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/workflow/rules/distances.smk b/workflow/rules/distances.smk index 65fe3d4..db1e798 100644 --- a/workflow/rules/distances.smk +++ b/workflow/rules/distances.smk @@ -1,15 +1,15 @@ rule compile_fail_sites_vcf: params: - header = ("#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"), - filter_text = "mask", + filter_text = "fail_site", sub_text = "NA", - exc_text = "site_qual" + exc_text = "site_qual", input: sites = OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv", output: sites = temp(OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf"), run: import pandas as pd + HEADER = ["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"] sites = ( pd.read_table(input.sites, sep="\t") .drop_duplicates(subset=("CHROM", "POS", "REF")) @@ -20,12 +20,10 @@ rule compile_fail_sites_vcf: sites["QUAL"] = "." sites["FILTER"] = params.filter_text sites["INFO"] = f"SUB={params.sub_text};EXC={params.exc_text}" - sites[list(params.header)].to_csv(output.sites, sep="\t", index=False) + sites[HEADER].to_csv(output.sites, sep="\t", index=False) -rule merge_sites: - params: - header = ("#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO") +rule merge_mask_sites_vcf: input: lambda wildcards: select_problematic_vcf(), OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf", @@ -33,15 +31,16 @@ rule merge_sites: sites = temp(OUTDIR / "all_mask_sites.vcf"), run: import pandas as pd + HEADER = ["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"] ( pd.concat( - [pd.read_table(path, sep="\t", comment="#", names=params.header) for path in input], + [pd.read_table(path, sep="\t", comment="#", names=HEADER, dtype={"POS": "int64"}) for path in input], axis="rows", ignore_index=True ) - .drop_duplicates(subset=("#CHROM", "POS", "FILTER"), keep="first") - .sort_values(list(params.header)) - .to_csv(output.sites, sep="\t", index=False) + .drop_duplicates(subset=("#CHROM", "POS", "FILTER"), keep="first") + .sort_values(by=["#CHROM", "POS"]) + .to_csv(output.sites, sep="\t", index=False) ) @@ -52,7 +51,7 @@ rule extract_afwdist_variants: position_col = "POS", sequence_col = "ALT", frequency_col = "ALT_FREQ", - mask_class = ["mask"], + mask_class = ["mask", "fail_site"], input: variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv", mask_vcf = OUTDIR / "all_mask_sites.vcf", From 28d55d2bed8ab347786def2955202866126fdc99 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20=C3=81lvarez=20Herrera?= Date: Tue, 10 Feb 2026 12:11:35 +0100 Subject: [PATCH 3/5] refactor: organize site quality filter in a new snakefile Reverts marking failed sites in filter with a different string (back to just "mask") --- workflow/core.smk | 1 + workflow/rules/distances.smk | 48 +------------- workflow/rules/sites.smk | 123 +++++++++++++++++++++++++++++++++++ workflow/rules/vaf.smk | 79 ---------------------- 4 files changed, 125 insertions(+), 126 deletions(-) create mode 100644 workflow/rules/sites.smk diff --git a/workflow/core.smk b/workflow/core.smk index 2bfe396..38ff64e 100644 --- a/workflow/core.smk +++ b/workflow/core.smk @@ -17,5 +17,6 @@ include: "rules/fetch.smk" include: "rules/fasta.smk" include: "rules/asr.smk" include: "rules/vaf.smk" +include: "rules/sites.smk" include: "rules/distances.smk" include: "rules/evolution.smk" diff --git a/workflow/rules/distances.smk b/workflow/rules/distances.smk index db1e798..77a4ae9 100644 --- a/workflow/rules/distances.smk +++ b/workflow/rules/distances.smk @@ -1,49 +1,3 @@ -rule compile_fail_sites_vcf: - params: - filter_text = "fail_site", - sub_text = "NA", - exc_text = "site_qual", - input: - sites = OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv", - output: - sites = temp(OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf"), - run: - import pandas as pd - HEADER = ["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"] - sites = ( - pd.read_table(input.sites, sep="\t") - .drop_duplicates(subset=("CHROM", "POS", "REF")) - .rename(columns={"CHROM": "#CHROM"}) - ) - sites["ID"] = "." - sites["ALT"] = "." - sites["QUAL"] = "." - sites["FILTER"] = params.filter_text - sites["INFO"] = f"SUB={params.sub_text};EXC={params.exc_text}" - sites[HEADER].to_csv(output.sites, sep="\t", index=False) - - -rule merge_mask_sites_vcf: - input: - lambda wildcards: select_problematic_vcf(), - OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf", - output: - sites = temp(OUTDIR / "all_mask_sites.vcf"), - run: - import pandas as pd - HEADER = ["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"] - ( - pd.concat( - [pd.read_table(path, sep="\t", comment="#", names=HEADER, dtype={"POS": "int64"}) for path in input], - axis="rows", - ignore_index=True - ) - .drop_duplicates(subset=("#CHROM", "POS", "FILTER"), keep="first") - .sort_values(by=["#CHROM", "POS"]) - .to_csv(output.sites, sep="\t", index=False) - ) - - rule extract_afwdist_variants: conda: "../envs/biopython.yaml" params: @@ -51,7 +5,7 @@ rule extract_afwdist_variants: position_col = "POS", sequence_col = "ALT", frequency_col = "ALT_FREQ", - mask_class = ["mask", "fail_site"], + mask_class = ["mask"], input: variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv", mask_vcf = OUTDIR / "all_mask_sites.vcf", diff --git a/workflow/rules/sites.smk b/workflow/rules/sites.smk new file mode 100644 index 0000000..493d06c --- /dev/null +++ b/workflow/rules/sites.smk @@ -0,0 +1,123 @@ +rule bcftools_mpileup_all_sites: + threads: 1 + conda: "../envs/var_calling.yaml" + params: + min_mq = 0, + min_bq = config["VC"]["MIN_QUALITY"], + mpileup_extra = "--no-BAQ" + input: + bam = get_input_bam, + reference = OUTDIR/"vaf"/"{sample}.reference.fasta", + output: + mpileup = temp(OUTDIR / "all_sites" / "{sample}.mpileup.vcf"), + query = temp(OUTDIR / "all_sites" / "{sample}.query.tsv"), + log: + mpileup = LOGDIR / "bcftools_mpileup_all_sites" / "{sample}.mpileup.txt", + query = LOGDIR / "bcftools_mpileup_all_sites" / "{sample}.query.txt", + shell: + "bcftools mpileup {params.mpileup_extra} -a AD,ADF,ADR --fasta-ref {input.reference:q} --threads {threads} -Q {params.min_bq} -q {params.min_mq} -Ov -o {output.mpileup:q} {input.bam:q} >{log.mpileup:q} 2>&1 && " + "echo 'CHROM\tPOS\tREF\tALT\tDP\tAD\tADF\tADR' >{output.query:q} && " + "bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%DP\t[ %AD]\t[ %ADF]\t[ %ADR]\n' {output.mpileup:q} >>{output.query:q} 2>{log.query:q}" + + +rule filter_mpileup_all_sites: + threads: 1 + params: + min_total_AD = config["VC"]["MIN_DEPTH"], + min_total_ADF = 0, + min_total_ADR = 0, + input: + OUTDIR / "all_sites" / "{sample}.query.tsv", + output: + sites_pass = temp(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv"), + sites_fail = temp(OUTDIR / "all_sites" / "{sample}.fail_sites.tsv"), + log: + LOGDIR / "filter_mpileup_all_sites" / "{sample}.txt" + run: + import pandas as pd + df = pd.read_csv(input[0], sep="\t") + df["SAMPLE"] = wildcards.sample + df["REF_AD"] = df.AD.str.split(",").apply(lambda values: int(values[0])) + df["TOTAL_AD"] = df.AD.str.split(",").apply(lambda values: sum(int(n) for n in values)) + df["TOTAL_ADF"] = df.ADF.str.split(",").apply(lambda values: sum(int(n) for n in values)) + df["TOTAL_ADR"] = df.ADR.str.split(",").apply(lambda values: sum(int(n) for n in values)) + mask = ( + (df.TOTAL_AD >= params.min_total_AD) & + (df.TOTAL_ADF >= params.min_total_ADF) & + (df.TOTAL_ADR >= params.min_total_ADR) + ) + df[mask].to_csv(output.sites_pass, sep="\t", index=False) + df[~mask].to_csv(output.sites_fail, sep="\t", index=False) + + +use rule concat_vcf_fields as merge_filtered_mpileup_all_sites with: + input: + expand(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv", sample=iter_samples()), + output: + OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv", + + +use rule concat_vcf_fields as merge_fail_mpileup_all_sites with: + input: + expand(OUTDIR / "all_sites" / "{sample}.fail_sites.tsv", sample=iter_samples()), + output: + OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv", + + +rule fill_all_sites: + conda: "../envs/renv.yaml" + input: + variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv", + sites = OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv", + output: + variants = OUTDIR/f"{OUTPUT_NAME}.variants.all_sites.tsv", + log: + LOGDIR / "fill_all_sites" / "log.txt" + script: + "../scripts/fill_all_sites.R" + + +rule compile_fail_sites_vcf: + params: + filter_text = "mask", + sub_text = "NA", + exc_text = "site_qual", + input: + sites = OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv", + output: + sites = temp(OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf"), + run: + import pandas as pd + HEADER = ["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"] + sites = ( + pd.read_table(input.sites, sep="\t") + .drop_duplicates(subset=("CHROM", "POS", "REF")) + .rename(columns={"CHROM": "#CHROM"}) + ) + sites["ID"] = "." + sites["ALT"] = "." + sites["QUAL"] = "." + sites["FILTER"] = params.filter_text + sites["INFO"] = f"SUB={params.sub_text};EXC={params.exc_text}" + sites[HEADER].to_csv(output.sites, sep="\t", index=False) + + +rule merge_mask_sites_vcf: + input: + lambda wildcards: select_problematic_vcf(), + OUTDIR / f"{OUTPUT_NAME}.fail_sites.vcf", + output: + sites = temp(OUTDIR / "all_mask_sites.vcf"), + run: + import pandas as pd + HEADER = ["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"] + ( + pd.concat( + [pd.read_table(path, sep="\t", comment="#", names=HEADER, dtype={"POS": "int64"}) for path in input], + axis="rows", + ignore_index=True + ) + .drop_duplicates(subset=("#CHROM", "POS", "FILTER"), keep="first") + .sort_values(by=["#CHROM", "POS"]) + .to_csv(output.sites, sep="\t", index=False) + ) diff --git a/workflow/rules/vaf.smk b/workflow/rules/vaf.smk index 6cccf20..b517ad9 100644 --- a/workflow/rules/vaf.smk +++ b/workflow/rules/vaf.smk @@ -204,82 +204,3 @@ use rule concat_vcf_fields as concat_variants with: expand(OUTDIR/"vaf"/"{sample}.variants.tsv", sample=iter_samples()), output: OUTDIR/f"{OUTPUT_NAME}.variants.tsv", - - -rule bcftools_mpileup_all_sites: - threads: 1 - conda: "../envs/var_calling.yaml" - params: - min_mq = 0, - min_bq = config["VC"]["MIN_QUALITY"], - mpileup_extra = "--no-BAQ" - input: - bam = get_input_bam, - reference = OUTDIR/"vaf"/"{sample}.reference.fasta", - output: - mpileup = temp(OUTDIR / "all_sites" / "{sample}.mpileup.vcf"), - query = temp(OUTDIR / "all_sites" / "{sample}.query.tsv"), - log: - mpileup = LOGDIR / "bcftools_mpileup_all_sites" / "{sample}.mpileup.txt", - query = LOGDIR / "bcftools_mpileup_all_sites" / "{sample}.query.txt", - shell: - "bcftools mpileup {params.mpileup_extra} -a AD,ADF,ADR --fasta-ref {input.reference:q} --threads {threads} -Q {params.min_bq} -q {params.min_mq} -Ov -o {output.mpileup:q} {input.bam:q} >{log.mpileup:q} 2>&1 && " - "echo 'CHROM\tPOS\tREF\tALT\tDP\tAD\tADF\tADR' >{output.query:q} && " - "bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%DP\t[ %AD]\t[ %ADF]\t[ %ADR]\n' {output.mpileup:q} >>{output.query:q} 2>{log.query:q}" - - -rule filter_mpileup_all_sites: - threads: 1 - params: - min_total_AD = config["VC"]["MIN_DEPTH"], - min_total_ADF = 0, - min_total_ADR = 0, - input: - OUTDIR / "all_sites" / "{sample}.query.tsv", - output: - sites_pass = temp(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv"), - sites_fail = temp(OUTDIR / "all_sites" / "{sample}.fail_sites.tsv"), - log: - LOGDIR / "filter_mpileup_all_sites" / "{sample}.txt" - run: - import pandas as pd - df = pd.read_csv(input[0], sep="\t") - df["SAMPLE"] = wildcards.sample - df["REF_AD"] = df.AD.str.split(",").apply(lambda values: int(values[0])) - df["TOTAL_AD"] = df.AD.str.split(",").apply(lambda values: sum(int(n) for n in values)) - df["TOTAL_ADF"] = df.ADF.str.split(",").apply(lambda values: sum(int(n) for n in values)) - df["TOTAL_ADR"] = df.ADR.str.split(",").apply(lambda values: sum(int(n) for n in values)) - mask = ( - (df.TOTAL_AD >= params.min_total_AD) & - (df.TOTAL_ADF >= params.min_total_ADF) & - (df.TOTAL_ADR >= params.min_total_ADR) - ) - df[mask].to_csv(output.sites_pass, sep="\t", index=False) - df[~mask].to_csv(output.sites_fail, sep="\t", index=False) - - -use rule concat_vcf_fields as merge_filtered_mpileup_all_sites with: - input: - expand(OUTDIR / "all_sites" / "{sample}.filtered_sites.tsv", sample=iter_samples()), - output: - OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv", - - -use rule concat_vcf_fields as merge_fail_mpileup_all_sites with: - input: - expand(OUTDIR / "all_sites" / "{sample}.fail_sites.tsv", sample=iter_samples()), - output: - OUTDIR / f"{OUTPUT_NAME}.fail_sites.tsv", - - -rule fill_all_sites: - conda: "../envs/renv.yaml" - input: - variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv", - sites = OUTDIR / f"{OUTPUT_NAME}.filtered_sites.tsv", - output: - variants = OUTDIR/f"{OUTPUT_NAME}.variants.all_sites.tsv", - log: - LOGDIR / "fill_all_sites" / "log.txt" - script: - "../scripts/fill_all_sites.R" From 960920818b92f9e30625f93c64646b115c4ee5ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20=C3=81lvarez=20Herrera?= Date: Tue, 10 Feb 2026 12:43:03 +0100 Subject: [PATCH 4/5] refactor: calculate days after the first sample consistently --- workflow/scripts/calculate_dnds.R | 8 +++++--- workflow/scripts/report/time_signal_data.R | 8 +++++--- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/workflow/scripts/calculate_dnds.R b/workflow/scripts/calculate_dnds.R index f545c27..a0495cb 100644 --- a/workflow/scripts/calculate_dnds.R +++ b/workflow/scripts/calculate_dnds.R @@ -29,9 +29,11 @@ variants <- read_delim( log_info("Reading metadata table") metadata <- read_delim(snakemake@input[["metadata"]]) %>% mutate( - interval = as.numeric( - as.Date(CollectionDate) - min(as.Date(CollectionDate)) - ) + interval = difftime( + as.Date(CollectionDate), + min(as.Date(CollectionDate)), + units = "days" + ) |> as.numeric() ) %>% select(ID, interval) %>% rename(SAMPLE = ID) diff --git a/workflow/scripts/report/time_signal_data.R b/workflow/scripts/report/time_signal_data.R index d209528..04448b2 100644 --- a/workflow/scripts/report/time_signal_data.R +++ b/workflow/scripts/report/time_signal_data.R @@ -44,9 +44,11 @@ time.signal <- distRoot( by = "ID" ) %>% mutate( - date_interval = as.numeric( - as.Date(CollectionDate) - min(as.Date(CollectionDate)) - ) + date_interval = difftime( + as.Date(CollectionDate), + min(as.Date(CollectionDate)), + units = "days" + ) |> as.numeric() ) # Save table From fb8be5e15ab6d09351673b68c2c11153e211b8d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miguel=20=C3=81lvarez=20Herrera?= Date: Tue, 10 Feb 2026 12:43:30 +0100 Subject: [PATCH 5/5] style: update phrasing about evolutionary rate --- template.qmd | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/template.qmd b/template.qmd index 8e3b6c3..c2e2fdf 100644 --- a/template.qmd +++ b/template.qmd @@ -77,7 +77,6 @@ stats <- append( ) correlation <- stats[["r2"]] sub_rate <- stats[["sub_rate"]] -sub_rate <- display.num(sub_rate, 2) p_value_lm <- stats[["pvalue"]] # NV counts @@ -262,7 +261,7 @@ frequency-weighted distances.](`r params$tree`){#fig-tree} To estimate the evolutionary rate, root-to-tip distances measured on the previous tree (@fig-tree) have been correlated with time, obtaining a $R^2$ of $`r display.num(correlation, 4)`$ and a p-value of $`r p_value_lm`$. The estimated evolutionary -rate is $`r sub_rate`$ number of changes per year (@fig-tempest). +rate is $`r display.num(sub_rate, 2)`$ changes per year (@fig-tempest). ![Scatterplot depicting the relationship between root-to-tip distances and the number of days passed since the first sample. The solid