31:T450,viewof table = Inputs.table(filteredSummary, { columns: [ "gene", "sites", "N", "T", "positive_sites", "negative_sites", "diff_sites", "comparison_group", "group_N", "group_T", "group_dN/dS", ], header: { "N": "branches", "T": "tree length", "positive_sites": "+ sites", "negative_sites": "- sites", "diff_sites": "diff sites", "comparison_group": "gtype", "group_N": "branches(gtype)", "group_T": "tree length(gtype)", "group_dN/dS": "dN/dS(gtype)", }, width: { "gene": 240, "sites": 50, "N": 60, "T": 70, "positive_sites": 50, "negative_sites": 50, "diff_sites": 60, "comparison_group": 50, }, format: { "group_T": (x) => d3.format(".4f")(x), "group_dN/dS": (x) => d3.format(".4f")(x), "BUSTED_pval": (x) => d3.format(".4")(x), "BUSTED_omega3": (x) => d3.format(".4")(x), "BUSTED_prop_sites_in_omega3": (x) => d3.format(".4")(x), "RELAX_clade_K": (x) => d3.format(".4")(x), "RELAX_overall_pval": (x) => d3.format(".4")(x), } })32:Taab,vegalite({ "$schema": "https://vega.github.io/schema/vega-lite/v5.json", "data": { "values": combined_summary }, "facet": { "field": "gene", "type": "nominal", "title": null, }, "columns": 5, "spec": { "width": 200, // Change this value to increase/decrease facet width "height": 200, // Optional: Adjust height if needed "layer": [ { "mark": "bar", "transform": [ {"calculate": "log(datum.RELAX_K)", "as": "log_k"}, { "calculate": "max(-2, min(2, log(datum.RELAX_K)))", "as": "log_k_clamped" } ], "encoding": { "x": { "field": "comparison_group", "type": "nominal", "title": "Genotype", "axis": {"labelAngle": 0} }, "y": { "field": "log_k_clamped", "type": "quantitative", "title": "log(RELAX K)", "scale": { "domain": [ -2, 2 ] } }, "color": { "condition": {"test": "datum.RELAX_overall_pval < 0.01", "value": "darkorange"}, "value": "steelblue" } } }, { "mark": { "type": "text", "dy": -5, "fontSize": 12 }, "transform": [ {"calculate": "log(datum.RELAX_K)", "as": "log_k"}, { "calculate": "max(-2, min(2, log(datum.RELAX_K)))", "as": "log_k_clamped" } ], "encoding": { "x": {"field": "comparison_group", "type": "nominal"}, "y": {"field": "log_k_clamped", "type": "quantitative"}, "text": {"field": "RELAX_K", "type": "quantitative", "format": ".3f"} } }, { "mark": { "type": "text", "dy": -100, "dx": -60, "fontSize": 14, //"color": "red", "fontWeight": "bold" }, "transform": [ {"calculate": "'p=' + (datum.RELAX_overall_pval >= 0.01 ? format(datum.RELAX_overall_pval, '.2f') : format(datum.RELAX_overall_pval, '.2e'))", "as": "formatted_pval"} ], "encoding": { "y": {"datum": -0.3}, "text": {"field": "formatted_pval", "type": "nominal", //"color": ("formatted_pval" >= 0.01 ? "black" : "red") } } } ] } })33:Ta49,**FEL** infers nonsynoymous (dN) and synonymous (dS) substitution rates for each site in a gene alignment. FEL assumes that the selection pressure for each site is constant across the entire phylogeny; _i.e._, a site must have evidence of negative selection across _all branches_ to have dN/dS < 1. Here, we report 'positive' if dN/dS > 1 and p<0.05, 'negative' if dN/dS < 1 and p<0.05, and 'neutral' if p>=0.05. **MEME** indicates that a site has been subject to episodic positive selection at some point in its history. Unlike BUSTED, MEME is not a gene-wide test for positive selection and is not subject to multiple testing correction, so positive results may be false positives and should be treated with some caution if BUSTED didn't show a significant result for the gene in question. MEME is also a better choice than FEL when searching for evidence of site-specific positive selection, because MEME does not require that all branches have evidence of positive selection for a site to be significant. **Contrast-FEL** searches for evidence that a site may be associated with adaptation to a different environment. A significant result at a site means that dN/dS (β/α) is different between the two sets of branches, with either an increase or a decrease on the foreground branches relative to the reference branches (Q-value < 0.05). A significant finding does not make any claims about positive (dN/dS > 1) or negative (dN/dS < 1), just that dN/dS differs among sets of branches, _i.e._, a difference need not change the mode of selection. Non-significant results do not necessarily mean that there is no difference in selection pressure between foreground and reference branches; it is also possible that whatever difference there may be does not rise to the level of statistical significance. **PRIME** searches for evidence that a site is experiencing selective pressure to maintain a particular physicochemical property. We have chosen to focus on the set of five biochemical properties described by Atchley et. al. (2005) when assessing property conservation at each site. **Intensified positive selection:** A site with evidence of episodic positive selection (_i.e._, a significant MEME p-value) and a statistically significant increase in the nonsynonymous substitution rate in the foreground relative to the reference (Contrast-FEL significant p-value and Beta_foreground > Beta_reference) may be subject to intensified positive selection. Such sites may also have different majority amino acid residues in the foreground and reference branches. Filters with a * require that Contrast-FEL results are present.34:T476,vegalite({ $schema: 'https://vega.github.io/schema/vega-lite/v5.json', description: 'Treemap of Composition', title: "Site composition", width: 350, // Change this value to increase/decrease facet width height: 350, // Optional: Adjust height if needed data: { values: focal_site_composition }, mark: { type: 'rect'}, encoding: { x: { field: 'category', type: 'ordinal', title: "Amino acid", axis: {"labelAngle": 0, "labelFontSize": 15}, sort: {x: "y"} }, y: { field: 'size', type: 'quantitative', title: "# Branches" }, // size: { field: 'size', type: 'quantitative' }, text: { field: 'size', type: 'quantitative' } }, layer: [ { mark: 'rect', encoding: { color: { field: 'category', type: 'nominal', legend: false }, size: { field: 'size', type: 'quantitative' } } }, { mark: { type: 'text', align: 'center', baseline: 'bottom' }, encoding: { text: { field: 'size', type: 'quantitative' }, } }, ] })35:T49b,vegalite({ $schema: 'https://vega.github.io/schema/vega-lite/v5.json', description: 'most common substitutions at the current site', title: "Most Common Substitutions", width: 350, // Change this value to increase/decrease facet width height: 350, // Optional: Adjust height if needed data: { values: focal_site_substitutions }, mark: { type: 'rect'}, encoding: { x: { field: 'category', type: 'ordinal', title: "Amino acid", axis: {"labelAngle": 0, "labelFontSize": 15}, sort: {x: "y"} }, y: { field: 'size', type: 'quantitative', title: "Substitutions" }, // size: { field: 'size', type: 'quantitative' }, text: { field: 'size', type: 'quantitative' } }, layer: [ { mark: 'rect', encoding: { color: { field: 'category', type: 'nominal', legend: false }, size: { field: 'size', type: 'quantitative' } } }, { mark: { type: 'text', align: 'center', baseline: 'bottom' }, encoding: { text: { field: 'size', type: 'quantitative' }, } }, ] })36:Taf5,serotypes = [ {"serotype": "DENV1", "combined_summary": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV1/combined_summary.csv", d3.autoType), "combined_sites": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV1/combined_sites.csv", d3.autoType), "combined_comparison_summary": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV1/combined_comparison_summary.csv", d3.autoType), "combined_comparison_sites": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV1/combined_comparison_site.csv", d3.autoType)}, {"serotype": "DENV2", "combined_summary": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV2/combined_summary.csv", d3.autoType), "combined_sites": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV2/combined_sites.csv", d3.autoType), "combined_comparison_summary": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV2/combined_comparison_summary.csv", d3.autoType), "combined_comparison_sites": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV2/combined_comparison_site.csv", d3.autoType)}, {"serotype": "DENV3", "combined_summary": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV3/combined_summary.csv", d3.autoType), "combined_sites": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV3/combined_sites.csv", d3.autoType), "combined_comparison_summary": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV3/combined_comparison_summary.csv", d3.autoType), "combined_comparison_sites": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV3/combined_comparison_site.csv", d3.autoType)}, {"serotype": "DENV4", "combined_summary": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV4/combined_summary.csv", d3.autoType), "combined_sites": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV4/combined_sites.csv", d3.autoType), "combined_comparison_summary": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV4/combined_comparison_summary.csv", d3.autoType), "combined_comparison_sites": await d3.csv("https://raw.githubusercontent.com/hverdonk/DENV_genotype_data/refs/heads/main/DENV4/combined_comparison_site.csv", d3.autoType)}, ]37:T438,// // Set of gene+site pairs that (1) have significant MEME results, // // (2) have significant CFEL results, and (3) Foreground cfel_beta > Reference cfel_beta // intensifiedSites = { // const byKey = d3.group(sites, siteKey); // Map key -> rows for that site (both groups) // const out = new Set(); // for (const [key, rows] of byKey) { // // must have a significant MEME result at this site // const hasMeme = rows.some(r => r.meme_qval r.cfel_marker < q_cut); // if (!hasCFEL) continue; // // find Foreground and Reference rows // const fg = rows.find(r => r.comparison_group === "Foreground"); // const bg = rows.find(r => r.comparison_group === "Reference"); // if (!fg || !bg) continue; // // numeric compare (coerce just in case) // const fgBeta = +fg.cfel_beta; // const bgBeta = +bg.cfel_beta; // if (Number.isFinite(fgBeta) && Number.isFinite(bgBeta) && fgBeta > bgBeta) { // out.add(key); // } // } // return out; // }38:T4eb,selectionDataAll = sitesForPlottingAll .map(d => { const site_num = +d.site; if (!Number.isFinite(site_num)) return null; const cfel_raw = +d.cfel_marker; const fel_raw = d.fel_qval != null && d.fel_qval !== "" ? +d.fel_qval : +d.fel_pval; const meme_raw = d.meme_qval != null && d.meme_qval !== "" ? +d.meme_qval : +d.meme_marker; const prime_raw = +d.prime_qval; const fel_beta_raw = +d.fel_beta; const fel_alpha_raw = +d.fel_alpha; const fel_beta_eps = (Number.isFinite(fel_beta_raw) ? fel_beta_raw : 0) + epsilon; const fel_alpha_eps = (Number.isFinite(fel_alpha_raw) ? fel_alpha_raw : 0) + epsilon; const fel_effect = Math.log(fel_beta_eps / fel_alpha_eps); return [ { method: "CFEL", raw_val: cfel_raw }, { method: "FEL", raw_val: fel_raw }, { method: "MEME", raw_val: meme_raw }, { method: "PRIME", raw_val: prime_raw } ].map(m => ({ ...d, site_num, site_start: site_num - 0.5, site_end: site_num + 0.5, method: m.method, raw_val: m.raw_val, significant: Number.isFinite(m.raw_val) && m.raw_val <= q_cut, fel_alpha_raw, fel_beta_raw, fel_effect })); }) .filter(Boolean) .flat()39:T1233,selectionGeneSpec = { return function selectionGeneSpec(gene, values, { showYAxis = true } = {}) { const maxSite = d3.max(values, d => d.site_num); const width = selectionPanelWidth(maxSite); return { title: { text: selectionGeneDisplayName(gene), anchor: "middle", fontSize: 16, fontWeight: "bold", offset: 8 }, vconcat: selectionSerotypeOrder.map((serotype, index) => ({ title: { text: serotype, anchor: "start", fontSize: 12, fontWeight: "bold", offset: 6 }, width, height: 88, data: {values: values.filter(d => d.serotype === serotype)}, layer: [ { transform: [{filter: "datum.significant && datum.method !== 'FEL'"}], mark: {type: "tick", orient: "vertical", thickness: 2.2, size: 22}, encoding: { x: { field: "site_num", type: "quantitative", scale: {domain: [0.5, maxSite + 0.5], nice: false}, axis: index === selectionSerotypeOrder.length - 1 ? {title: "Site", titleAnchor: "end", titlePadding: 10} : null }, y: { field: "method", type: "nominal", sort: selectionMethodOrder, scale: { domain: selectionMethodOrder, paddingInner: 0.3, paddingOuter: 0.08 }, axis: showYAxis ? {title: null} : {title: null, labels: false, ticks: false, domain: true, tickSize: 6} }, color: { field: "method", type: "nominal", scale: { domain: ["FEL", "MEME", "PRIME", "CFEL"], range: ["#56B4E9", "#E69F00", "#009E73", "#D55E00"] }, legend: null }, tooltip: [ {field: "serotype", title: "Serotype"}, {field: "gene", title: "Gene"}, {field: "site_num", title: "Site"}, {field: "method", title: "Method"}, {field: "raw_val", title: "q (or p fallback)"} ] } }, { transform: [{filter: "datum.significant && datum.method === 'FEL' && isFinite(datum.fel_effect)"}], mark: {type: "rect"}, encoding: { x: { field: "site_start", type: "quantitative", scale: {domain: [0.5, maxSite + 0.5], nice: false}, axis: index === selectionSerotypeOrder.length - 1 ? {title: "Site", titleAnchor: "end", titlePadding: 10} : null }, x2: {field: "site_end"}, y: { field: "method", type: "nominal", sort: selectionMethodOrder, scale: { domain: selectionMethodOrder, paddingInner: 0.3, paddingOuter: 0.08 }, axis: showYAxis ? {title: null} : {title: null, labels: false, ticks: false, domain: true, tickSize: 6} }, color: { field: "fel_effect", type: "quantitative", scale: { domain: [-selectionMaxAbsFelEffectAll, selectionMaxAbsFelEffectAll], scheme: "redblue", domainMid: 0, reverse: true }, legend: null }, tooltip: [ {field: "serotype", title: "Serotype"}, {field: "gene", title: "Gene"}, {field: "site_num", title: "Site"}, {field: "method", title: "Method"}, {field: "raw_val", title: "q (or p fallback)"}, {field: "fel_alpha_raw", title: "FEL alpha"}, {field: "fel_beta_raw", title: "FEL beta"}, {field: "fel_effect", title: "log(ω)", format: ".3f"} ] } } ] })), spacing: 8 }; } }3a:T810,/** * Apply per-gene BH to MEME and FEL. * * Adds: * - meme_q_bh * - fel_q_bh * * Options: * - geneKey: field to group by (default "gene") * - methods: array of { pField, qField } definitions * - filter: optional predicate to select rows (e.g. Foreground only) */ function addPerGeneBH(sites, { geneKey = "gene", methods = [ { pField: "meme_marker", qField: "meme_qval" }, { pField: "fel_pval", qField: "fel_qval" } // <-- change if your FEL p-value field is named differently ], filter = null } = {}) { // Work on a shallow copy so we don't mutate caller data const out = sites.map(d => ({ ...d })); // Helper: parse p-values that might be "-", null, etc. const toP = (v) => { if (v === null || v === undefined) return null; if (typeof v === "string") { const t = v.trim(); if (t === "" || t === "-") return null; const n = Number(t); return Number.isFinite(n) ? n : null; } if (typeof v === "number") return Number.isFinite(v) ? v : null; const n = Number(v); return Number.isFinite(n) ? n : null; }; // Group row indices by gene const geneToIndices = new Map(); for (let i = 0; i < out.length; i++) { if (filter && !filter(out[i])) continue; const g = out[i][geneKey]; if (!g) continue; if (!geneToIndices.has(g)) geneToIndices.set(g, []); geneToIndices.get(g).push(i); } // For each gene, compute BH per method for (const [gene, indices] of geneToIndices.entries()) { for (const { pField, qField } of methods) { const indexedPvals = []; for (const idx of indices) { const p = toP(out[idx][pField]); // Accept only p in [0,1] if (p === null) continue; if (p < 0 || p > 1) continue; indexedPvals.push({ idx, p }); } const qMap = benjaminiHochbergForIndices(indexedPvals); // Write q-values back; leave missing as null for (const idx of indices) { out[idx][qField] = qMap.has(idx) ? qMap.get(idx) : null; } } } return out; }3b:[] 2f:{"title":"Dengue genotype-specific evolutionary analysis results","cells":[{"id":0,"value":"# Dengue genotype-specific evolutionary analysis results","pinned":false,"hidden":false,"mode":"md"},{"id":985,"value":"viewof focal_serotype_data = Inputs.select(serotypes, {label: \"Select serotype\", value: d => d, format: x => x.serotype})","pinned":false,"hidden":false,"mode":"ojs"},{"id":117,"value":"## Data summary","pinned":false,"hidden":false,"mode":"md"},{"id":305,"value":"BUSTED, FEL, MEME, and PRIME were all run on the full serotype, and should be identical to [what we've seen previously](https://observablehq.com/d/562fa8ed5fbda59f).\n\n`+ sites` indicates how many sites had evidence of positive selection (significant MEME results)\n\n`- sites` indicates how many sites had evidence of negative selection (significant FEL results with site-level dN/dS < 1)\n\n`diff_sites` indicates how many sites had significant test results for contrast-FEL, indicating that selection may be acting differently for those sites across clades.","pinned":false,"hidden":false,"mode":"md"},{"id":42,"value":"//viewof selectedClade = Inputs.select([\"All\", \"DENV1\", \"DENV2\", \"DENV3\", \"DENV4\"], {label: \"Filter by clade\"})\nviewof sites_selections_0 = Inputs.form({\n option1: Inputs.select([\"All\", ...comparison_groups], {label: \"Filter by group\"}),\n option2: Inputs.select([\"All\", ...new Set(summary.map(d => d.gene))], {label: \"Filter by gene\"}),\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":29,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":107,"value":"## Gene-wide evidence for positive selection","pinned":false,"hidden":false,"mode":"md"},{"id":109,"value":"BUSTED (Murrell et al. 2015) was used to search for gene-wide evidence of positive selection for at least one site on at least one branch of the provided gene tree. BUSTED estimates the proportion of sites belonging to each of three rate classes, constrained as ω1 ≤ ω2 ≤ 1 ≤ ω3. The positive class is the rate class where we assign sites with evidence of positive selection, for which we estimate ω3.\n\nA **p-value** < 0.05 indicates that BUSTED detected statistically significant evidence of positive selection. **omega3** is the evolutionary rate ratio estimated from all the sites in the positive class. **% of sites in positive class** indicates how many sites in the original alignment contributed to the evidence for positive selection. BUSTED itself isn't designed to detect *which* sites may be under positive selection; see **[Per-site evidence of selection](https://observablehq.com/d/562fa8ed5fbda59f#cell-146)** instead.","pinned":false,"hidden":false,"mode":"md"},{"id":120,"value":"viewof keep_significant = Inputs.toggle({label: html`Show only significant results`, value: true, default: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":63,"value":"viewof table2 = Inputs.table(BUSTED_results, {\n columns: [\n \"gene\",\n \"BUSTED_pval\",\n \"BUSTED_omega3\",\n \"BUSTED_prop_sites_in_omega3\",\n ],\n header: {\n \"N\": \"internal branches\",\n \"T\": \"tree length\",\n \"BUSTED_pval\": \"P-value\",\n \"BUSTED_omega3\": \"omega3\",\n \"BUSTED_prop_sites_in_omega3\": \"% of sites in positive class\",\n \n },\n format: {\n \"BUSTED_pval\": (x) => d3.format(\".3\")(x),\n \"BUSTED_omega3\": (x) => d3.format(\".4\")(x),\n //\"BUSTED_prop_sites_in_omega3\": (x) => d3.format(\".4\"), //(x * 100),\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":92,"value":"## Relaxed / intensified selection across genotypes","pinned":false,"hidden":false,"mode":"md"},{"id":94,"value":"The K-value is the relaxation/intensification parameter from RELAX, plotted here on a log scale. Lower K-values (log(K) < 0) indicate relaxed selection compared to the Reference branches, while higher K-values (log(K) > 0) indicate intensified - but not *necessarily* positive - selection. p < 0.01 indicates that the differences are significant (orange bars).","pinned":false,"hidden":false,"mode":"md"},{"id":1040,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":146,"value":"## Per-site evidence of selection","pinned":false,"hidden":false,"mode":"md"},{"id":1274,"value":"Results of site-level selection tests are shown below, faceted by gene and by serotype. A tick mark indicates that a site had significant evidence of selection under the particular method. **FEL:** Sites are colored by the direction and strength of selection: blue for negative selection and red for positive selection. **MEME:** Colored sites indicate episodic positive selection at that position. **PRIME:** Colored sites indicate selection to maintain one or more biochemical properties (bipolar, structure, volume, composition, charge, or overall (the shared effect of all five properties)). **Contrast-FEL:** Colored sites indicate a difference in selection pressure for at least one genotype among all tested genotypes for each serotype. The number of sequences that mapped to the NS3 gene in the DENV 2 serotype was insufficient to build a gene tree, so no testing could be performed.","pinned":false,"hidden":false,"mode":"md"},{"id":1155,"value":"viewof q_cut = Inputs.range([0, 1], {value: 0.2, step: 0.05, label: \"Significance threshold (q-value):\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1161,"value":"vegalite({\n \"$schema\": \"https://vega.github.io/schema/vega-lite/v5.json\",\n \"hconcat\": [\n felLegendSpecAll,\n {\n \"vconcat\": selectionRowsAll,\n \"spacing\": 14\n }\n ],\n \"spacing\": 20,\n \"config\": {\n \"font\": \"Helvetica\",\n \"axis\": {\n \"labelFontSize\": 14,\n \"titleFontSize\": 13,\n \"tickSize\": 3,\n \"domain\": true,\n \"grid\": false\n },\n \"view\": { \"stroke\": null }\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":148,"value":"### Selection across all branches","pinned":false,"hidden":false,"mode":"md"},{"id":295,"value":"$33","pinned":false,"hidden":false,"mode":"md"},{"id":168,"value":"viewof sites_selections_1 = Inputs.form({\n option1: Inputs.select([\"All\", ...comparison_groups], {label: \"Filter by group\"}),\n option2: Inputs.select([\"All\", ...new Set(serotypeSitesWithBH.map(d => d.gene))], {label: \"Filter by gene\"}),\n option3: Inputs.toggle({label: \"Show only significant FEL results\", value: false, default: false}),\n option4: Inputs.toggle({label: \"Show only significant MEME results\", value: false, default: false}),\n option5: Inputs.toggle({label: \"*Show only significant Contrast-FEL results\", value: false, default: false}),\n option6: Inputs.toggle({label: \"Show only significant PRIME results\", value: false, default: false}),\n //option7: Inputs.toggle({label: \"*Show intensified positive selection\", value: false, default: false}),\n option8: Inputs.toggle({label: \"*Keep only sites with different majority residues\", value: false, default: false}),\n // option3: Inputs.select([\"All\", ...new Set(sites.map(d => d.site))], {label: \"Filter by site\"}),\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":281,"value":"viewof table3 = Inputs.table(filteredSites, {\n columns: [\n \"gene\",\n \"comparison_group\",\n \"site\",\n \"majority_residue\",\n \"fel_selection\",\n \"fel_qval\",\n \"meme_qval\",\n \"cfel_marker\",\n \"cfel_beta\",\n \"prime_marker\"\n ],\n header: {\n \"comparison_group\": \"branch set\",\n \"majority_residue\": \"majority residue\",\n \"fel_selection\": \"FEL results\",\n \"fel_qval\": \"FEL qval\",\n \"meme_qval\": \"MEME qval\",\n \"cfel_marker\": \"CFEL qval\",\n \"cfel_beta\": \"CFEL Beta\",\n \"prime_marker\": \"PRIME properties under selection\",\n },\n width: {\n \"gene\": 200,\n \"site\": 30,\n // \"cfel_marker\": 120,\n // \"cfel_beta\": 120,\n \"prime_marker\": 250\n },\n format: {\n \"meme_marker\": (x) => formatMarker(x),\n \"cfel_marker\": (x) => formatMarker(x),\n \"cfel_beta\": (x) => d3.format(\".4\")(x),\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":256,"value":"Site composition indicates the prevalence of each amino acid at a given site. It does not include how many sites had missing data (Codon \"NNN\"), as branches with missing site data could not contribute to the evolutionary analyses and were ignored for the given site.","pinned":false,"hidden":false,"mode":"md"},{"id":195,"value":"//viewof choose_site = Inputs.select([...new Set(filteredSites.map(d => d.site))], {label: \"Sites in current serotype + gene with significant results\"})\nviewof choose_gene = Inputs.select([...new Set(serotypeSitesWithBH.map(d => d.gene))], {label: \"Filter by gene\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":468,"value":"// viewof choose_group = Inputs.select(\n// [...comparison_groups], \n// {label: \"Filter by group\"}\n// )\nviewof choose_group = comparison_groups.size > 0\n ? Inputs.select([...comparison_groups], {value: \"Foreground\", label: \"Filter by group\"})\n : Inputs.select([], {label: \"Filter by group\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":459,"value":"viewof choose_site = Inputs.select(\n [...new Set(serotypeSitesWithBH\n .filter(d => (!choose_group || d.comparison_group === choose_group) &&\n (!choose_gene || d.gene === choose_gene))\n .map(d => d.site))], \n {value: 2, label: \"Filter by site\"}\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":155,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":279,"value":"The following chart indicates the most common substitutions of one amino acid for another at the selected genotype and site. Hyphy analyses use _time reversible_ codon evolutionary models, so these substitutions should NOT be viewed as directional. For example, a substitution listed as Q:R could equally be written as R:Q. The total number of substitutions would reflect the sum of both Q:R and R:Q substitutions. ","pinned":false,"hidden":false,"mode":"md"},{"id":273,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":126,"value":"## Data","pinned":false,"hidden":false,"mode":"md"},{"id":983,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":837,"value":"combined_summary = focal_serotype_data.combined_summary\n//combined_summary = combined_summary_file.csv({typed: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":839,"value":"combined_sites = focal_serotype_data.combined_sites\n// combined_sites = combined_sites_file.csv({typed: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":841,"value":"combined_comparison_summary = focal_serotype_data.combined_comparison_summary\n// combined_comparison_summary = (\n// combined_comparison_summary_file && combined_comparison_summary_file.name.endsWith(\".csv\")\n// )\n// ? await combined_comparison_summary_file.csv({ typed: true })\n// : []","pinned":false,"hidden":false,"mode":"ojs"},{"id":843,"value":"combined_comparison_sites = focal_serotype_data.combined_comparison_sites\n// combined_comparison_sites = (\n// combined_comparison_sites_file && combined_comparison_sites_file.name.endsWith(\".csv\")\n// )\n// ? await combined_comparison_sites_file.csv({ typed: true })\n// : []","pinned":false,"hidden":false,"mode":"ojs"},{"id":621,"value":"comparison_groups = {\n const groups = new Set(summary.map(d => d.comparison_group));\n groups.delete(\"background\");\n return groups\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":20,"value":"//summary = FileAttachment(\"capsid_protein_C_summary@2.csv\").csv()\n// match gene AND comparison_group\nsummary = mergePrimaryIntoSecondary(combined_summary, combined_comparison_summary, { key: d => `${d.gene}::${d.comparison_group}`}) // combineJsonData(summary_urls)\n .filter(d => d.comparison_group != 'Nuisance') // remove nuisance background (leaf) branches\n .map(r => (r.comparison_group && r.comparison_group.toLowerCase() === \"reference\") \n ? { ...r, RELAX_K: 1 } // reset RELAX K of 'Reference' to 1 (RELAX default for reference branches)\n : r\n )","pinned":false,"hidden":false,"mode":"ojs"},{"id":54,"value":"filteredSummary = summary\n .filter(d => sites_selections_0.option1 === \"All\" || d.comparison_group === sites_selections_0.option1)\n .filter(d => sites_selections_0.option2 === \"All\" || d.gene === sites_selections_0.option2);","pinned":false,"hidden":false,"mode":"ojs"},{"id":130,"value":"BUSTED_results = keep_significant\n ? combined_summary.filter(d => d.BUSTED_pval < 0.05).filter(d => d.comparison_group === [...comparison_groups][0])\n : combined_summary.filter(d => d.comparison_group === [...comparison_groups][0]);","pinned":false,"hidden":false,"mode":"ojs"},{"id":17,"value":"// sites = combineJsonData(sites_urls)\n// merge on gene AND site number\nserotypeSites = mergePrimaryIntoSecondary(combined_sites, combined_comparison_sites, {key: siteKey})\n .filter(d => d.comparison_group != 'Nuisance')","pinned":false,"hidden":false,"mode":"ojs"},{"id":1150,"value":"// serotypeSitesWithBH = addPerGeneBH(serotypeSites, {\n// geneKey: \"gene\",\n// methods: [\n// { pField: \"meme_marker\", qField: \"meme_qval\" },\n// { pField: \"fel_pval\", qField: \"fel_qval\" }\n// ],\n// });\nserotypeSitesWithBH = serotypeSites","pinned":false,"hidden":false,"mode":"ojs"},{"id":182,"value":"filteredSites = serotypeSitesWithBH\n .filter(d => sites_selections_1.option1 === \"All\" || d.comparison_group === sites_selections_1.option1)\n .filter(d => sites_selections_1.option2 === \"All\" || d.gene === sites_selections_1.option2)\n .filter(d => !sites_selections_1.option3 || d.fel_qval < q_cut)\n .filter(d => !sites_selections_1.option4 || d.meme_qval !sites_selections_1.option5 || d.cfel_marker < q_cut)\n .filter(d => !sites_selections_1.option6 || d.prime_qval < q_cut)\n //.filter(d => !sites_selections_1.option7 || intensifiedSites.has(siteKey(d)))\n .filter(d => !sites_selections_1.option8 || majorityDiffPairs.has(siteKey(d)));","pinned":false,"hidden":false,"mode":"ojs"},{"id":723,"value":"$37","pinned":true,"hidden":false,"mode":"ojs"},{"id":787,"value":"// Set of (gene,site) pairs whose majority_residue differs among the set of tested genotypes\nmajorityDiffPairs = {\n const byKey = d3.group(serotypeSitesWithBH, siteKey); // Map\n const out = new Set();\n\n for (const [key, rows] of byKey) {\n // collect distinct residues (after normalizing via residue())\n const distinct = new Set(\n rows\n .map(r => residue(r.majority_residue))\n .filter(v => v != null && v !== \"\")\n );\n\n // if any differs from the others, there will be >1 distinct residue\n if (distinct.size > 1) out.add(key);\n }\n \n return out;\n};","pinned":false,"hidden":false,"mode":"ojs"},{"id":207,"value":"//focal_site = sites.filter(d => d.site === choose_site && d.comparison_group === choose_group && d.gene === choose_gene)\nfocal_site = serotypeSitesWithBH.filter(d => \n d.site === choose_site &&\n d.gene === choose_gene &&\n (!choose_group || d.comparison_group === choose_group)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":213,"value":"focal_site_composition = {\n const raw = focal_site?.[0]?.composition;\n if (!raw || raw === \"-\") return [];\n\n return raw\n .split(/\\s*,\\s*/) // e.g., [\"E:9\", \"T:1\"]\n .filter(Boolean)\n .map(tok => {\n const i = tok.lastIndexOf(\":\");\n if (i <= 0) return null; // malformed\n const category = tok.slice(0, i).trim();\n const size = Number(tok.slice(i + 1));\n return Number.isFinite(size) ? { category, size } : null;\n })\n .filter(Boolean);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":260,"value":"focal_site_substitutions = {\n const s = focal_site?.[0]?.substitutions;\n if (!s || s === \"-\") return [];\n\n return s\n .split(/\\s*,\\s*/) // [\"E:E:9\", \"E:T:1\"]\n .filter(Boolean)\n .map(entry => {\n const i = entry.lastIndexOf(\":\");\n const category = entry.slice(0, i).trim(); // \"E:E\"\n const size = Number(entry.slice(i + 1)); // 9\n return { category, size };\n })\n .filter(d => Number.isFinite(d.size));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1190,"value":"### Summary plot data","pinned":false,"hidden":false,"mode":"md"},{"id":1249,"value":"#### Load cross-serotype sites","pinned":false,"hidden":false,"mode":"md"},{"id":1163,"value":"combined_sites_all = serotypes.flatMap(({serotype, combined_sites}) => combined_sites.map(d => ({...d, serotype})))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1187,"value":"combined_comparison_sites_all = serotypes.flatMap(({serotype, combined_comparison_sites}) => combined_comparison_sites.map(d => ({...d, serotype})))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1251,"value":"#### Merge cross-serotype sites and apply per-serotype, per-gene Benjamini-Hotchberg multiple testing correction to MEME and FEL","pinned":false,"hidden":false,"mode":"md"},{"id":1192,"value":"siteKeyAll = d => `${d.serotype}::${d.gene}::${d.site}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1196,"value":"sitesAll = mergePrimaryIntoSecondary(combined_sites_all, combined_comparison_sites_all, {key: siteKeyAll})\n .filter(d => d.comparison_group != 'Nuisance')","pinned":false,"hidden":false,"mode":"ojs"},{"id":1199,"value":"sitesWithBHAll = addPerGeneBH(sitesAll, {\n geneKey: d => `${d.serotype}::${d.gene}`,\n methods: [\n { pField: \"meme_marker\", qField: \"meme_qval\" },\n { pField: \"fel_pval\", qField: \"fel_qval\" }\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1254,"value":"#### Prep data for plotting","pinned":false,"hidden":false,"mode":"md"},{"id":1203,"value":"plotFieldsAll = [\n \"serotype\",\n \"gene\",\n \"site\",\n \"cfel_marker\",\n \"fel_qval\",\n \"fel_pval\",\n \"meme_qval\",\n \"meme_marker\",\n \"prime_qval\",\n \"fel_beta\",\n \"fel_alpha\"\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1206,"value":"// remove genotype-level data & drop duplicate rows for easier plotting of summary chart\nsitesForPlottingAll = {\n const seen = new Set();\n\n return sitesWithBHAll\n .map(row => Object.fromEntries(plotFieldsAll.map(field => [field, row[field] ?? null])))\n .filter(row => {\n const key = plotFieldsAll.map(field => JSON.stringify(row[field] ?? null)).join(\"::\");\n if (seen.has(key)) return false;\n seen.add(key);\n return true;\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1210,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"},{"id":1216,"value":"selectionByGeneAll = d3.group(selectionDataAll, d => d.gene)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1219,"value":"// cap for FEL color range\nselectionMaxAbsFelEffectAll = d3.max(\n selectionDataAll\n .filter(d => d.method === \"FEL\" && d.significant && Number.isFinite(d.fel_effect))\n .map(d => Math.abs(d.fel_effect))\n) ?? 1","pinned":false,"hidden":false,"mode":"ojs"},{"id":1257,"value":"#### Configure plot","pinned":false,"hidden":false,"mode":"md"},{"id":1223,"value":"selectionMethodOrder = [\"FEL\", \"MEME\", \"PRIME\", \"CFEL\"]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1226,"value":"selectionSerotypeOrder = [\"DENV1\", \"DENV2\", \"DENV3\", \"DENV4\"]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1229,"value":"selectionDefaultRows = [\n [\"capsid_protein_C\", \n \"membrane_glycoprotein_precursor_prM\", \n \"envelope_protein_E\",\n \"nonstructural_protein_NS1\", \n \"nonstructural_protein_NS2A\",\n \"nonstructural_protein_NS2B\",\n \"nonstructural_protein_NS3\", \n \"nonstructural_protein_NS4A\", \n \"nonstructural_protein_NS4B\",\n \"RNA-dependent_RNA_polymerase_NS5\"]\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1232,"value":"selectionMinWidth = 160","pinned":false,"hidden":false,"mode":"ojs"},{"id":1235,"value":"selectionMaxWidth = 520","pinned":false,"hidden":false,"mode":"ojs"},{"id":1238,"value":"selectionPanelWidth = {\n const geneSizes = Array.from(selectionByGeneAll, ([gene, values]) => ({\n gene,\n nSites: d3.max(values, d => d.site_num)\n }));\n\n const minSites = d3.min(geneSizes, d => d.nSites);\n const maxSites = d3.max(geneSizes, d => d.nSites);\n const widthScale = d3.scaleSqrt()\n .domain([minSites, maxSites])\n .range([selectionMinWidth, selectionMaxWidth]);\n\n return nSites => Math.round(widthScale(nSites));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1241,"value":"selectionGeneDisplayName = gene => ({\n \"capsid_protein_C\": \"C\",\n \"membrane_glycoprotein_precursor_prM\": \"prM\",\n \"envelope_protein_E\": \"E\",\n \"nonstructural_protein_NS1\": \"NS1\",\n \"nonstructural_protein_NS2A\": \"NS2A\",\n \"nonstructural_protein_NS2B\": \"NS2B\",\n \"nonstructural_protein_NS3\": \"NS3\",\n \"nonstructural_protein_NS4A\": \"NS4A\",\n \"nonstructural_protein_NS4B\": \"NS4B\",\n \"RNA-dependent_RNA_polymerase_NS5\": \"NS5\"\n}[gene] ?? gene)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1246,"value":"$39","pinned":false,"hidden":false,"mode":"ojs"},{"id":1264,"value":"selectionRowsAll = selectionDefaultRows.map(rowGenes => ({\n hconcat: rowGenes\n .filter(gene => selectionByGeneAll.has(gene))\n .map((gene, colIndex) => selectionGeneSpec(gene, selectionByGeneAll.get(gene), { showYAxis: colIndex === 0 })),\n spacing: 24\n}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1267,"value":"felLegendSpecAll = {\n return {\n width: 1,\n height: 1,\n data: {\n values: [\n {fel_effect: -selectionMaxAbsFelEffectAll},\n {fel_effect: selectionMaxAbsFelEffectAll}\n ]\n },\n mark: {type: \"point\", opacity: 0, filled: true},\n encoding: {\n color: {\n field: \"fel_effect\",\n type: \"quantitative\",\n scale: {\n domain: [-selectionMaxAbsFelEffectAll, selectionMaxAbsFelEffectAll],\n scheme: \"redblue\",\n domainMid: 0,\n reverse: true\n },\n legend: {\n title: \"FEL log(ω)\",\n orient: \"left\",\n direction: \"vertical\",\n gradientLength: 180,\n gradientThickness: 16\n }\n }\n }\n };\n }","pinned":false,"hidden":false,"mode":"ojs"},{"id":7,"value":"### Functions & Imports","pinned":false,"hidden":false,"mode":"md"},{"id":9,"value":"vegalite = require(\"vega-embed\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":707,"value":"siteKey = d => `${d.gene}::${d.site}`;","pinned":false,"hidden":false,"mode":"ojs"},{"id":789,"value":"// Helper to normalize residue letters (and ignore missing/unknown)\nresidue = x => {\n if (x == null) return null;\n const r = String(x).trim().toUpperCase();\n return (r === \"-\" || r === \"?\") ? null : r;\n};","pinned":false,"hidden":false,"mode":"ojs"},{"id":739,"value":"formatMarker = {\n const fmt = d3.format(\".4~g\");\n return x => (x == null || x === \"-\" ? x : fmt(+x));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":584,"value":"// Merge \"primary\" fields into each row of \"secondary\" (if it exists).\n// Matching key defaults to `gene`, but you can pass a custom key fn.\nfunction mergePrimaryIntoSecondary(primary, secondary, { key = d => d.gene } = {}) {\n // If there's no secondary, return primary as-is\n if (!Array.isArray(secondary) || secondary.length === 0) return primary;\n\n // Index primary by the chosen key\n const idx = new Map(primary.map(p => [key(p), p]));\n\n // For each secondary row, overlay its fields on top of the matching primary row\n return secondary.map(s => {\n const p = idx.get(key(s));\n // If no match, keep the secondary row; otherwise merge (secondary wins on conflicts)\n return p ? { ...p, ...s } : s;\n });\n}\n\n// --- Usage ---\n// Basic (match by gene)\n// const merged = mergePrimaryIntoSecondary(firstArray, secondArray);\n\n// // Stricter match (gene + fasta_id)\n// const mergedStrict = mergePrimaryIntoSecondary(\n// firstArray,\n// secondArray,\n// { key: d => `${d.gene}::${d.fasta_id}` }\n// );","pinned":false,"hidden":false,"mode":"ojs"},{"id":403,"value":"async function combineJsonData(urls) {\n const combinedData = [];\n\n // Fetch and combine data from each URL\n for (const url of urls) {\n try {\n const data = await d3.csv(url, d3.autoType);\n combinedData.push(...data); // Combine arrays\n } catch (error) {\n console.error(`Failed to load data from ${url}:`, error);\n }\n }\n\n return combinedData;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":322,"value":"function pivotCfelMarker(data) {\n return data.map(item => {\n const newItem = { ...item };\n if (item.cfel_marker && item.cfel_marker !== '-') {\n const markers = item.cfel_marker.split(',').map(m => m.trim());\n markers.forEach(marker => {\n const [key, value] = marker.split(':').map(s => s.trim());\n newItem[key] = value;\n });\n }\n delete newItem.cfel_marker;\n return newItem;\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1058,"value":"first_genotype = Array.from(comparison_groups).sort()[0]","pinned":false,"hidden":false,"mode":"ojs"},{"id":956,"value":"// determines the extent of the log(RELAX_K values) on the RELAX chart, so we can add a little padding to the extent. This way, K and p values aren't clipped by the x-axis\nextent = d3.extent(\n combined_summary\n .map(d => Math.log(d.RELAX_K))\n .filter(v => Number.isFinite(v))\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1142,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs"},{"id":1144,"value":"/**\n * Benjamini–Hochberg (BH) correction.\n * Input: array of objects that include a numeric p-value field.\n * Output: Map from row index -> q-value (BH-adjusted p).\n */\nfunction benjaminiHochbergForIndices(indexedPvals) {\n // indexedPvals: [{ idx, p }]\n const m = indexedPvals.length;\n if (m === 0) return new Map();\n\n // Sort ascending by p\n const sorted = [...indexedPvals].sort((a, b) => a.p - b.p);\n\n // Compute raw BH q = p * m / rank\n const qRaw = new Array(m);\n for (let i = 0; i < m; i++) {\n const rank = i + 1; // 1-based\n qRaw[i] = (sorted[i].p * m) / rank;\n }\n\n // Enforce monotonicity from the end: q_i = min(q_i, q_{i+1})\n const qMon = new Array(m);\n qMon[m - 1] = Math.min(qRaw[m - 1], 1);\n for (let i = m - 2; i >= 0; i--) {\n qMon[i] = Math.min(qRaw[i], qMon[i + 1], 1);\n }\n\n // Map back to original indices\n const out = new Map();\n for (let i = 0; i < m; i++) {\n out.set(sorted[i].idx, qMon[i]);\n }\n return out;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1212,"value":"// tiny offset value so we don't have to worry about alpha or beta being 0 when we calculate log(beta/alpha)\nepsilon = 1e-12","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W3b","stdlib":"1","titling":"auto"}