31:Tc67,plot_power_curve = (data, x_title, title) => { title = title || "Power"; return vegalite ({"data": {values : data}, "layer" : [ { "mark": { "type" : "area", "opacity" : 0.2, "interpolate" : "monotone" }, "width": 500, "height":300, "encoding": { "x": { "field": "x", "type": "quantitative", "axis":{ "title" : x_title, grid : false, titleFontSize: 20, labelFontSize: 14, format : "r" } }, "y2": { "field": "p10", "type": "quantitative" }, "y": { "axis":{ "title" : title, grid : false, titleFontSize: 20, labelFontSize: 14, format : "r" }, "field": "p1", "type": "quantitative", "scale" : { "type" : "linear", } } , "axis": null, } }, { "data" : {"values" : {"reference" : 0.05}}, "mark": { "type": "rule", "color" : "grey", "strokeDash" : [5,1] }, "encoding": { "y": { "field": "reference", "type": "quantitative" }, "size": {"value": 3} } }, { "mark": { "type" : "line", "strokeWidth" : 5, "interpolate" : "monotone" }, "encoding": { "x": { "field": "x", "type": "quantitative", }, "y": { "field": "p5", "type": "quantitative", "scale" : { "type" : "linear", } }, } }, { "mark": { "type" : "line", "strokeWidth" : 5, "color" : "red", "interpolate" : "monotone" }, "encoding": { "x": { "field": "x", "type": "quantitative" }, "y": { "field": "pp5", "type": "quantitative", "scale" : { "type" : "linear", } }, } }, { "mark": { "type" : "line", "strokeWidth" : 5, "color" : "grey", "interpolate" : "monotone" }, "encoding": { "x": { "field": "x", "type": "quantitative", }, "y": { "field": "q20", "type": "quantitative", "scale" : { "type" : "linear", } }, } }, { "mark": { "type" : "circle", "color" : "#f58518" }, "encoding": { "x": { "field": "x", "type": "quantitative", }, "y": { "field": "p5", "type": "quantitative", }, "size": { "field": "count", "type": "quantitative", "legend" : {"title" : "Samples", "titleFontSize": 20, "orient" : "right", "labelFontSize" : 14}, } } } ], "config": { "view": { "stroke": "transparent" } }}) }32:T40e,sim_data = { var data = raw_data; _.each (raw_data, (d)=>{ d.kind = (d.sim_alpha == d.sim_beta ? "NULL" : (d.sim_alpha > d.sim_beta_internals ? "-" : "+")); d.variable = (d.mle_alpha + d.mle_beta) > 0; d.rate_diff = Math.abs (d.sim_beta_internals-d.sim_alpha); d.logalpha = Math.log10 (d.sim_alpha); d.logbeta = Math.log10 (d.sim_beta_internals); d.omega = d.sim_beta_internals/d.sim_alpha; d.branches = (d.FileID.includes("all") ? "all" : "internal"); d.geneID = (d.SimID.includes("S") ? "S" : (d.SimID.includes("M") ? "M" : (d.SimID.includes("N") ? "N" : d.SimID.includes("ORF1a") ? "ORF1a": (d.SimID.includes("ORF1b") ? "ORF1b" : (d.SimID.includes("ORF3a") ? "ORF3a" : (d.SimID.includes("ORF6") ? "ORF6": d.SimID.includes("ORF7a") ? "ORF7a": "ORF8" )))))); if (d.p <= 0.1) { d.call10 = true; if (d.p <= 0.05) { d.call05 = true; if (d.p <= 0.01) { d.call01 = true; } } } }); return data; }33:T402,compute_metrics = (data, accessor, p_crit, label) => { var false_positives = 0, false_negatives = 0, true_positives = 0, true_negatives = 0; _.each (data, (d) => { var p = accessor (d); var l = label (d); if (l) { if (p <= p_crit) { true_positives++; } else { false_negatives++; } } else { if (p <= p_crit) { false_positives++; } else { true_negatives++; } } }); var N = true_positives + false_positives + true_negatives + false_negatives; return { 'p' : p_crit, 'N' : N, 'tp' : true_positives, 'fp' : false_positives, 'tn' : true_negatives, 'fn' : false_negatives, 'fpr' : false_positives / (true_negatives + false_positives), 'power' : true_positives / (true_positives + false_negatives), 'ppv' : true_positives / (true_positives + false_positives) }; };34:T47a,compute_roc = (data, accessor, step, label) => { var aoc_by_step = [[0,0,0,0]]; var sorted_list = []; var total_positives = 0, total_negatives = 0; _.each (data, (d) => { var p = accessor (d); var l = label (d); if (l) { total_positives ++; } else { total_negatives ++; } sorted_list.push ([p,l]); }); sorted_list.sort ((a,b) => a[0]-b[0]); var index = 0; var true_calls = 0, false_calls = 0; for (var p = step; p <= 1 + (step/2); p+=step) { while (index < sorted_list.length && sorted_list[index][0] <= p) { if (sorted_list[index][1]) { true_calls++; } else { false_calls ++; } index ++; } aoc_by_step.push ([p, true_calls+false_calls,true_calls/total_positives,false_calls/total_negatives, true_calls/(true_calls+false_calls)]); } return _.map (aoc_by_step, (d) => {return {'fp' : d[3], 'power' : d[2], 'ppv' : d[4], "p" : d[0], "b" : total_positives/(total_positives + total_negatives)}}); };35:Te84,plot_prc = (plot_data, labels, cv) => { cv = cv || [[0.01,true],[0.05,true],[0.1,true]]; var merged_data = []; var merged_baseline = []; const fmt = d3.format (".2f"); _.each (plot_data, (d,i) => { merged_baseline.push ({'reference' : d[0]['b'], 'label' : labels[i]}); _.each (d, (e) => { e.label = labels[i]; e.pair = "R=" + fmt(e.power) + ", P=" + fmt (e.ppv); merged_data.push (e); }); }); var crticial_values = _.filter (merged_data, (d)=>d3.min (_.map (cv, (c)=>Math.abs (c[0]-d.p))) < 1e-5); var cv_shown = _.filter (cv, c=>c[1]); var text_to_show = _.filter (merged_data, (d)=>d3.min (_.map (cv_shown, (c)=>Math.abs (c[0]-d.p))) < 1e-5); var baseline = plot_data[0][0]['b']; return vegalite ({ "layer" : [{ "data": {values : merged_data}, "mark": { "type": "line", "interpolate": "monotone", "strokeWidth" : 5 }, "width" : 500, "height" : 500, "encoding": { "x": { "field": "power", "axis":{ "title" : "Recall", grid : false, titleFontSize: 18, labelFontSize: 14, format : "g" }, "type": "quantitative" }, "y": { "field": "ppv", "type": "quantitative", "axis":{ "title" : "Precision", grid : false, titleFontSize: 18, labelFontSize: 14, format : "g" } }, "color" : { "field" : "label", "type" : "nominal", "legend" : {"title" : "Simulation", "titleFontSize": 1, "orient" : "right", "labelFontSize" : 14}, } } }, { "data": {values : crticial_values}, "mark": { "type": "point", "size" : 400, "filled" : true }, "encoding": { "x": { "field": "power", "axis":{ "title" : "Recall", grid : false, titleFontSize: 14, labelFontSize: 12, format : "g" }, "type": "quantitative" }, "y": { "field": "ppv", "type": "quantitative", "axis":{ "title" : "Precision", grid : false, titleFontSize: 16, labelFontSize: 12, format : "g" } }, "color" : { "field" : "label", "type" : "nominal" } } }, { "data" : {"values" : merged_baseline}, "mark": { "type": "rule", "strokeDash" : [5,1] }, "encoding": { "y": { "field": "reference", "type": "quantitative" }, "size": {"value": 3}, "color" : { "field" : "label", "type" : "nominal" } } }, { "data": {values : text_to_show}, "mark": { "type": "text", "dx" : 9, "dy" : 0, "fontSize" : 10, "align" : "left", //"angle" : 90. }, "encoding": { "x": { "field": "power", "axis":{ "title" : "Recall", grid : false, titleFontSize: 14, labelFontSize: 12, format : "g" }, "type": "quantitative" }, "y": { "field": "ppv", "type": "quantitative", "axis":{ "title" : "Precision", grid : false, titleFontSize: 16, labelFontSize: 12, format : "g" } }, /*"color" : { "field" : "label", "type" : "nominal" },*/ "text": { "field" : "pair", } } } ], "config": { "view": { "stroke": "transparent" } }}); }36:T576,plot_qq = (data, label) => { var color = d3.scaleOrdinal(d3.schemeCategory10); return vegalite ({"layer": _.map (data, (plot_data,i) => {return { "data": {values : plot_data}, "layer" : [{ "mark": { "type": "line", "interpolate": "monotone"//, //"color" : color (i) }, "width" : 300, "height" : 300, "encoding": { "x": { "field": "expected", "axis":{ "title" : "Nominal p-value", grid : false, titleFontSize: 14, labelFontSize: 12, format : "g" }, "type": "quantitative" }, "y": { "field": "observed", "type": "quantitative", "axis":{ "title" : "Observed rejection rate", grid : false, titleFontSize: 16, labelFontSize: 12, format : "g" } }, "color" : { "field" : "geneID" } } }, { "mark": { "type": "line", "interpolate": "monotone", "color" : "grey", "strokeDash" : [3,3] }, "width" : 800, "height" : 800, "encoding": { "x": { "field": "expected", "type": "quantitative" }, "y": { "field": "expected", "type": "quantitative", } } } ], "config": { "view": { "stroke": "transparent" } }}})}); }37:T44e,by_N = { const p = 0.05; var bin_map = function (n) { if (n <= 8) { return "" + n; } if (n <= 16) { let n2 = 2*Math.floor (n/2.); return "" + n2 + "-" + (n2 + 1); } if (n <= 32) { let n4 = 4*Math.floor (n/4.); return "" + n4 + "-" + (n4 + 3); } if (n <= 64) { let n8 = 8*Math.floor (n/8.); return "" + n8 + "-" + (n8 + 7); } return ">64"; }; var binned = {}; // 0 - TN // 1 - FP // 2 - FN // 3 - TP _.each (power_sims_detectable, (d) => { let k = bin_map (d.subs_test); if (! (k in binned)) { binned [k] = [0,0,0,0,0,0]; } binned [k][2] ++; if (d.p < p) { binned [k][3] ++; } }); _.each (null_sims_detectable, (d) => { let k = bin_map (d.subs_test); if (! (k in binned)) { binned [k] = [0,0,0,0,0,0]; } binned [k][0] ++; if (d.p < p) { binned [k][1] ++; } }); _.each (binned, (d) => {d[4] = d[3]/(d[1]+d[3]); d[5] = d[0]/(d[2]+d[0]);}); return binned; };38:T4f8,bin_detection = (data, field, step, from, to) => { from = from || 0; to = to || d3.max (data, (d)=>d[field]); var values = []; var current_floor = from; for (var upto = step + from; upto < to + step * 0.5; upto += step) { values.push ({'count' : 0, 'x' : upto - step * 0.5, 'p5' : 0, 'p1' : 0, 'p10' : 0, 'pp5' : 0, 'q20' : 0 }); } _.each (data, (record)=> { var index = Math.floor ((record[field]-from)/step); if (index < 0) index = 0; if (index >= values.length) index = values.length - 1; values[index].count++; if (record.call05) { values[index].p5++ }; if (record.call01) { values[index].p1++ }; if (record.call10) { values[index].p10++ }; if (record.pcall) { values[index].pp5++ }; if (record.qcall20) { values[index].q20++ }; }); _.each (values, (record) => { _.each (['p1','p5','p10','pp5','q20'], (l) => { record[l] /= record.count; }); }); return values; }39:[] 2f:{"title":"SARS-CoV-2 simulation data","cells":[{"id":0,"value":"md`# SARS-CoV-2 simulation data\n\nExploration of simulation results for the FEL method.\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1444,"value":"viewof branchType= select({\n title: \"branchType\",\n options:[\"all\", \"internal\"]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1447,"value":"viewof geneType= select({\n title: \"geneID\",\n options:[\"N\", \"M\", \"S\", \"ORF1a\", \"ORF1b\", \"ORF3a\", \"ORF6\", \"ORF7a\", \"ORF8\"]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1449,"value":"plot_qq ([_.filter (_.map (compute_roc (_.filter (sim_data, (d) => d.variable && d.kind != \"+\" && d.branches == branchType && d.geneID == geneType), (d)=>d.p, 0.005, (d)=>d.sim_alpha < d.sim_beta), (d)=>{return {\"expected\" : d.p, \"observed\" : d.fp}}), (d) => d.expected <= 0.5),_.filter (_.map (compute_roc (_.filter (sim_data, (d) => d.variable && d.kind == \"SAME\" && d.max_rate > 2.5 && d.max_rate < 3.5), (d)=>d.p, 0.005), (d)=>{return {\"expected\" : d.p, \"observed\" : d.fp}}), (d) => d.expected <= 0.5)])","pinned":false,"hidden":false,"mode":"ojs"},{"id":1414,"value":"md`qqplot for entire data set`","pinned":false,"hidden":false,"mode":"ojs"},{"id":384,"value":"plot_qq ([_.filter (_.map (compute_roc (_.filter (sim_data, (d) => d.variable && d.kind != \"+\"), (d)=>d.p, 0.005, (d)=>d.sim_alpha < d.sim_beta), (d)=>{return {\"expected\" : d.p, \"observed\" : d.fp}}), (d) => d.expected <= 0.5),_.filter (_.map (compute_roc (_.filter (sim_data, (d) => d.variable && d.kind == \"SAME\" && d.max_rate > 2.5 && d.max_rate < 3.5), (d)=>d.p, 0.005), (d)=>{return {\"expected\" : d.p, \"observed\" : d.fp}}), (d) => d.expected <= 0.5)])","pinned":false,"hidden":false,"mode":"ojs"},{"id":1417,"value":"md` qqplot for Internal flag only`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1410,"value":"plot_qq ([_.filter (_.map (compute_roc (_.filter (sim_data, (d) => d.variable && d.kind != \"+\" && d.branches == \"internal\"), (d)=>d.p, 0.005, (d)=>d.sim_alpha < d.sim_beta), (d)=>{return {\"expected\" : d.p, \"observed\" : d.fp}}), (d) => d.expected <= 0.5),_.filter (_.map (compute_roc (_.filter (sim_data, (d) => d.variable && d.kind == \"SAME\" && d.max_rate > 2.5 && d.max_rate < 3.5), (d)=>d.p, 0.005), (d)=>{return {\"expected\" : d.p, \"observed\" : d.fp}}), (d) => d.expected <= 0.5)])","pinned":false,"hidden":false,"mode":"ojs"},{"id":1422,"value":"md` qqplot for All flag only`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1424,"value":"plot_qq ([_.filter (_.map (compute_roc (_.filter (sim_data, (d) => d.variable && d.kind != \"+\" && d.branches == \"all\"), (d)=>d.p, 0.005, (d)=>d.sim_alpha < d.sim_beta), (d)=>{return {\"expected\" : d.p, \"observed\" : d.fp}}), (d) => d.expected <= 0.5),_.filter (_.map (compute_roc (_.filter (sim_data, (d) => d.variable && d.kind == \"SAME\" && d.max_rate > 2.5 && d.max_rate < 3.5), (d)=>d.p, 0.005), (d)=>{return {\"expected\" : d.p, \"observed\" : d.fp}}), (d) => d.expected <= 0.5)])","pinned":true,"hidden":false,"mode":"ojs"},{"id":1351,"value":"_.filter (sim_data, (d) => d.variable)","pinned":true,"hidden":false,"mode":"ojs"},{"id":239,"value":"plot_power_curve (bin_detection (_.filter (power_sims_detectable, (d)=>true), \"rate_diff\",1, 0, 10), \"Number of substitutions in the test set\", \"Detection rate\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":555,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":15,"value":"md`------`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1406,"value":"_.filter (sim_data, (d) => d.branches == \"internal\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":19,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":287,"value":"raw_data = await d3.csv (\"https://dl.dropbox.com/s/m2iw53bj22ieqtp/temp_res_5_22_20-2.csv\", d3.autoType)//await d3.csv (\"https://dl.dropbox.com/s/pstq3h9s9o468un/results-with-fel.csv\",d3.autoType);","pinned":false,"hidden":false,"mode":"ojs"},{"id":278,"value":"null_sims_detectable = _.filter (sim_data, d=>d.kind == \"SAME\" && (d.mle_alpha + d.mle_beta))","pinned":false,"hidden":false,"mode":"ojs"},{"id":280,"value":"power_sims_detectable = _.filter (sim_data, d=>d.kind == \"+\" && (d.mle_alpha + d.mle_beta))","pinned":false,"hidden":false,"mode":"ojs"},{"id":137,"value":"p_exp = _.filter (sim_data, d => d.kind == 'SAME' && (d.sim_alpha > 0 || d.sim_beta_leaves > 0 || d.sim_beta_internal>0)).length * binning_step","pinned":false,"hidden":false,"mode":"ojs"},{"id":520,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":300,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":330,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":381,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":6,"value":"_=require ('lodash')","pinned":false,"hidden":false,"mode":"ojs"},{"id":9,"value":"d3 = require(\"d3\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":140,"value":"binning_step = 0.01","pinned":false,"hidden":false,"mode":"ojs"},{"id":251,"value":"binning_step_power_curve = 0.5","pinned":false,"hidden":false,"mode":"ojs"},{"id":10,"value":"vegalite = require(\"vega-embed@3\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":201,"value":"by_sim = {\n const p = 0.05;\n \n var binned = {};\n _.each (sim_data, (d) => {\n if (! (d.SimID in binned)) {\n binned [d.SimID] = [0,0,0,d.SimID];\n }\n if (d.variable && d.sim_beta_internals == d.sim_alpha) {\n binned [d.SimID][0] ++;\n if (d.p < p) {\n binned [d.SimID][1] ++;\n }\n }\n });\n \n _.each (binned, (d) => d[2] = d[1]/d[0]);\n \n return binned;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":266,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":514,"value":"import { table } from \"@gampleman/table\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1430,"value":"import {select, slider, checkbox} from \"@jashkenas/inputs\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":756,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"},{"id":1043,"value":"float_f = d3.format (\".3g\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1250,"value":"NP = _.filter (sim_data, (d) => d.variable && d.beta1 != d.beta2).length","pinned":false,"hidden":false,"mode":"ojs"},{"id":1255,"value":"NN=_.filter (sim_data, (d) => d.variable && d.beta1 == d.beta2).length","pinned":false,"hidden":false,"mode":"ojs"},{"id":1259,"value":"cutoff = 0.10","pinned":false,"hidden":false,"mode":"ojs"},{"id":1329,"value":"has_fel = _.filter (sim_data, (d)=>(\"FEL_TEST_p\" in d))","pinned":true,"hidden":false,"mode":"ojs"},{"id":1266,"value":"CM_CF = [[_.filter (sim_data, (d) => d.variable && d.sim_beta_internals == d.sim_beta_leaves && d.p > cutoff).length, _.filter (sim_data, (d) => d.variable && d.sim_beta_internals == d.sim_beta_leaves && d.p < cutoff).length],[_.filter (sim_data, (d) => d.variable && d.sim_beta_internals != d.sim_beta_leaves && d.p > cutoff).length, _.filter (sim_data, (d) => d.variable && d.sim_beta_internals != d.beta2 && d.p <= cutoff).length]]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1281,"value":"same_class = (d)=> {return d.FEL_TEST_p <= cutoff && d.FEL_REF_p <= cutoff && d.FEL_TEST_diff * d.FEL_REF_diff > 0 || d.FEL_TEST_p > cutoff && d.FEL_REF_p > cutoff}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1277,"value":"CM_SF = [[_.filter (sim_data, (d) => d.variable && d.sim_beta_internals == d.sim_beta_leaves && same_class(d)).length, _.filter (sim_data, (d) => d.variable && d.beta1 == d.beta2 && ! same_class(d)).length],[_.filter (sim_data, (d) => d.variable && d.sim_beta_internals != d.sim_beta_leaves && same_class(d)).length, _.filter (sim_data, (d) => d.variable && d.sim_beta_internals != d.sim_beta_leaves && !same_class(d)).length]]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1301,"value":"fp = _.map ([CM_CF,CM_SF],(m)=> m[1][1]/(m[0][1]+m[1][1]))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1305,"value":"ppv = _.map ([CM_CF,CM_SF],(m)=> m[0][0]/(m[0][0]+m[1][0]))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1297,"value":"_.filter (sim_data, (d) => d.variable && d.sim_beta_internals == d.sim_beta_leaves && d.p > cutoff && !same_class (d))","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W39","stdlib":"1","titling":"auto"}