2f:T56f, ** Table 1. ** Column definitions are as follows: 1. **ID** the candidate gene/alignment under study 2. **Sites** the total number of codons in the alignment 3. **Total Sequences** the number of sequences in the alignment 4. **BUSTED-E filtered sites** total number of sites BUSTED-E masked as error-prone (where sites filtered in multiple species are counted more than once). 5. **Fraction of sites filtered** = (BUSTED-E filtered sites) ÷ (Sites * Total Sequences) 6. **Error P-value** P-value for BUSTED-E vs BUSTED LRT 7. **Filter?** interpretation of BUSTED-E vs BUSTED LRT; is filtered or unfiltered data recommended? 8. **aBSREL selected branches** number of branches in the test set (homo sapiens) under positive selection 9. **BUSTED(E?) selection p-val** P-value for the test for positive selection in genes. Either BUSTED or BUSTED-E: if Filter?(7) is "filtered", then BUSTED-E is used. If Filter?(7) is "unfiltered", then BUSTED-E is used 10. **BUSTED(E?) ⍵1** first ⍵ value designated as negative selection 11. **BUSTED(E?) ⍵2** second ⍵ value designated as neutral or negative selection 12. **BUSTED(E?) ⍵3** third ⍵ value designated as neutral or positive selection 13. **BUSTED(E?) p1** fraction of sites partitioned as ⍵1 14. **BUSTED(E?) p2** fraction of sites partitioned as ⍵2 15. **BUSTED(E?) p3** fraction of sites partitioned as ⍵3 30:T119c,table_data = { let count = 0; let for_table = [] for (const key in mod_results_json ) { let absrel_u_branches = _.filter(_.map(mod_results_json[key]['aBSREL_UNFILTERED']['branch attributes'], d=> d['Corrected P-value']), d => d != null && d < pval) let absrel_f_branches = _.filter(_.map(mod_results_json[key]['aBSREL_FILTERED']['branch attributes'], d=> d['Corrected P-value']), d => d != null && d < pval) // Avery LRT calc let avery_lrt_val = -2 * (mod_results_json[key]['BUSTED']['fits']['Unconstrained model']['Log Likelihood'] - mod_results_json[key]['BUSTED_E']['fits']['Unconstrained model']['Log Likelihood']) let avery_e_p_val = chisqrprob(2, avery_lrt_val) / 2 // Filter based on the pvals let avery_filt = avery_e_p_val < pval ? ['filtered', 'BUSTED-E'] : ['unfiltered', 'BUSTED'] // end lrt calc // AGL LRT calc // Based on https://sites.warnercnr.colostate.edu/gwhite/wp-content/uploads/sites/73/2017/04/LikelihoodRatioTests.pdf // Get Log likelihoods of the models let simple_model_ll = mod_results_json[key]['BUSTED']['fits']['Unconstrained model']['Log Likelihood']; let alternative_model_ll = mod_results_json[key]['BUSTED_E']['fits']['Unconstrained model']['Log Likelihood']; // you should always see LR ≥ 0 let lrt_val = -2 * simple_model_ll + 2 * alternative_model_ll; console.log("Simple model LL: " + simple_model_ll + " Complex model LL " + alternative_model_ll); console.log("Avery LRT: " + avery_lrt_val); console.log("AGL LRT: " + lrt_val); let e_p_val = chisqrprob(2, lrt_val)/2 // Filter based on the pvals let filt = e_p_val < pval ? ['filtered', 'BUSTED-E'] : ['unfiltered', 'BUSTED'] // end lrt calc if (lrt_val < 0) { count += 1; console.log ("Potential convergence issue! Avery! " + key); } const obj = {} //overall obj["ID"] = key obj["Total sites"] = mod_results_json[key]['BUSTED_E']['input']['number of sites'] obj["Total sequences"] = mod_results_json[key]['BUSTED_E']['input']['number of sequences'] //BUSTED-E obj["BUSTED-E filtered sites"] = _.reduce(mod_results_json[key]['FILTER']['filter'], (count, filter) => count + _.keys(filter).length, 0) obj["Fraction of sites filtered"] = obj["BUSTED-E filtered sites"]/(obj["Total sites"] * obj["Total sequences"]) obj["Error p-val"] = e_p_val obj["Filtered?"] = filt[0] //aBSREL obj["aBSREL selected branches"] = filt[0] === 'filtered' ? absrel_f_branches.length : absrel_u_branches.length //BUSTED obj["BUSTED selection p-val"] = filt[0] === 'filtered' ? mod_results_json[key]['BUSTED_E']['test results']['p-value'] : mod_results_json[key]['BUSTED']['test results']['p-value'] obj["BUSTED(E?) ⍵1"] = filt[0] === 'filtered' ? mod_results_json[key]['BUSTED_E']['fits']['Unconstrained model']['Rate Distributions']['Test'][1]['omega'] : mod_results_json[key]['BUSTED']['fits']['Unconstrained model']['Rate Distributions']['Test'][0]['omega'] obj["BUSTED(E?) ⍵2"] = filt[0] === 'filtered' ? mod_results_json[key]['BUSTED_E']['fits']['Unconstrained model']['Rate Distributions']['Test'][2]['omega'] : mod_results_json[key]['BUSTED']['fits']['Unconstrained model']['Rate Distributions']['Test'][1]['omega'] obj["BUSTED(E?) ⍵3"] = filt[0] === 'filtered' ? mod_results_json[key]['BUSTED_E']['fits']['Unconstrained model']['Rate Distributions']['Test'][3]['omega'] : mod_results_json[key]['BUSTED']['fits']['Unconstrained model']['Rate Distributions']['Test'][2]['omega'] obj["BUSTED(E?) p1"] = filt[0] === 'filtered' ? mod_results_json[key]['BUSTED_E']['fits']['Unconstrained model']['Rate Distributions']['Test'][1]['proportion'] : mod_results_json[key]['BUSTED']['fits']['Unconstrained model']['Rate Distributions']['Test'][0]['proportion'] obj["BUSTED(E?) p2"] = filt[0] === 'filtered' ? mod_results_json[key]['BUSTED_E']['fits']['Unconstrained model']['Rate Distributions']['Test'][2]['proportion'] : mod_results_json[key]['BUSTED']['fits']['Unconstrained model']['Rate Distributions']['Test'][1]['proportion'] obj["BUSTED(E?) p3"] = filt[0] === 'filtered' ? mod_results_json[key]['BUSTED_E']['fits']['Unconstrained model']['Rate Distributions']['Test'][3]['proportion'] : mod_results_json[key]['BUSTED']['fits']['Unconstrained model']['Rate Distributions']['Test'][2]['proportion'] for_table.push(obj) } console.log("Issues in " + count + " files"); return for_table }31:[] 2d:{"title":"Using BUSTED-E filtering script on Alyssa's hominin data","cells":[{"id":3,"value":"# Using BUSTED-E filtering script on Alyssa's hominin data ","pinned":false,"hidden":false,"mode":"md"},{"id":16,"value":"viewof search = Inputs.search(table_data)","pinned":false,"hidden":false,"mode":"ojs"},{"id":13,"value":"table1 = Inputs.table(search, {\n rows : 15\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":911,"value":"$2f","pinned":false,"hidden":false,"mode":"md"},{"id":631,"value":"viewof gene_id = {\n const dd3 = autoSelect ({\n title: \"Gene ID\",\n options: OnlySearchResults,\n placeholder: \"Select a gene to analyze . . .\",\n //value: 'NM_007118_5_23_species_23_seqs'\n \n });\n dd3.input.style.fontSize = \"12px\";\n dd3.input.style.marginTop = \"8px\";\n return dd3;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":402,"value":"viewof pval = Inputs.range([0, 1], {label: \"P-value threshold for error filtering LRT\", step: .01, value: 0.05})","pinned":false,"hidden":false,"mode":"ojs"},{"id":464,"value":"## Chosen gene: ${gene_id}","pinned":false,"hidden":false,"mode":"md"},{"id":773,"value":"### BUSTED-E RESULTS:","pinned":false,"hidden":false,"mode":"md"},{"id":697,"value":"BUSTED-E filtered out **${total_count}**, where **${unique_sites.length}** sites are unique (without repetitions across species) (**p-value ${Math.round(error_p_value*100)/100}**) sites as error-prone(sites with an unusually high omega). These sites were grouped together in their own omega error category to distinguish from the non-erroneous sites. ","pinned":false,"hidden":false,"mode":"md"},{"id":672,"value":"### SELECTION RESULTS","pinned":false,"hidden":false,"mode":"md"},{"id":930,"value":"#### *${which_busted}* results for the *gene-wide test for positive selection*:","pinned":false,"hidden":false,"mode":"md"},{"id":791,"value":"We performed a test for positive selection on the gene ${gene_id} (selected above) \n\n**p = ${matchingEntry['BUSTED selection p-val']}**: there is ${selection_evidence}.","pinned":false,"hidden":false,"mode":"md"},{"id":730,"value":"omega_plot.legend(\"color\", {\n label: \"Omega values →\",\n tickSize: 3,\n ticks: 10,\n //tickLabelRotate: 45,\n //tickFormat: r.tickFormat,\n //round: r.round,\n width: 600,\n height: 45\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":708,"value":"omega_plot = Plot.plot({\n color: {\n type: \"sqrt\",\n //pivot: 4,\n scheme: \"brbg\",\n reverse: true,\n //legend: true\n }, \n marks: [\n Plot.barX(stacked_omegas, {y: \"dist\",x: \"proportion\", fill: 'omega', title: (d) =>\n `${d.dist} omega bin #${d.omega_bin}\\n Omega value: ${d.omega},\\nProportion: ${d.proportion}` }),\n Plot.ruleX([0])\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":823,"value":"----","pinned":false,"hidden":false,"mode":"md"},{"id":933,"value":"#### *aBSREL* results for the *test for positive selection on Homo sapiens branches*:","pinned":false,"hidden":false,"mode":"md"},{"id":955,"value":"p-value: ${absrel_p_val} (test for positive selection on Homo sapiens branch)","pinned":false,"hidden":false,"mode":"md"},{"id":943,"value":"----","pinned":false,"hidden":false,"mode":"md"},{"id":964,"value":"absrel_p_val = matchingEntry['Filtered?'] === \"unfiltered\" ? results_json['NM_182961_6_22_species_22_seqs']['aBSREL_UNFILTERED']['branch attributes']['homo_sapiens']['Corrected P-value'] : results_json['NM_182961_6_22_species_22_seqs']['aBSREL_FILTERED']['branch attributes']['homo_sapiens']['Corrected P-value'] ","pinned":false,"hidden":false,"mode":"ojs"},{"id":825,"value":"#### Data parsing","pinned":false,"hidden":false,"mode":"md"},{"id":687,"value":"which_busted = error_p_value < 0.05 ? \"BUSTED_E\" : \"BUSTED\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":491,"value":"total_count = _.reduce(results_json[gene_id]['FILTER']['filter'], (count, filter) => count + _.keys(filter).length, 0);","pinned":false,"hidden":false,"mode":"ojs"},{"id":497,"value":"unique_sites = {\n let list = []\n let filter_list = mod_results_json[gene_id]['FILTER']['filter']\n for (let key in filter_list){\n for (let i =0; i< Object.keys(filter_list[key]).length; i++){\n list.push(filter_list[key][i])\n }\n }\n return [...new Set(list)]\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":680,"value":"error_p_value = chisqrprob(2, error_LRT)/2","pinned":false,"hidden":false,"mode":"ojs"},{"id":676,"value":"error_LRT = -2 * (mod_results_json[gene_id]['BUSTED']['fits']['Unconstrained model']['Log Likelihood'] - mod_results_json[gene_id]['BUSTED_E']['fits']['Unconstrained model']['Log Likelihood'])","pinned":false,"hidden":false,"mode":"ojs"},{"id":10,"value":"$30","pinned":true,"hidden":false,"mode":"ojs"},{"id":602,"value":"chisqrprob = function chisqrprob ($n,$x) {\n if (($n <= 0) || ((Math.abs($n) - (Math.abs(Math.floor($n)))) != 0)) {\n throw(\"Invalid n: $n\\n\"); /* degree of freedom */\n }\n return _subchisqrprob($n-0, $x-0);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":607,"value":"function _subuprob ($x) {\n var $p = 0; /* if ($absx > 100) */\n var $absx = Math.abs($x);\n\n if ($absx < 1.9) {\n $p = Math.pow((1 +\n $absx * (.049867347\n + $absx * (.0211410061\n + $absx * (.0032776263\n + $absx * (.0000380036\n + $absx * (.0000488906\n + $absx * .000005383)))))), -16)/2;\n } else if ($absx <= 100) {\n for (var $i = 18; $i >= 1; $i--) {\n $p = $i / ($absx + $p);\n }\n $p = Math.exp(-.5 * $absx * $absx) \n / Math.sqrt(2 * Math.PI) / ($absx + $p);\n }\n\n if ($x<0)\n $p = 1 - $p;\n return $p;\n}\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":604,"value":"function _subchisqrprob ($n,$x) {\n var $p;\n\n if ($x <= 0) {\n $p = 1;\n } else if ($n > 100) {\n $p = _subuprob((Math.pow(($x / $n), 1/3)\n - (1 - 2/9/$n)) / Math.sqrt(2/9/$n));\n } else if ($x > 400) {\n $p = 0;\n } else { \n var $a;\n var $i;\n var $i1;\n if (($n % 2) != 0) {\n $p = 2 * _subuprob(Math.sqrt($x));\n $a = Math.sqrt(2/Math.PI) * Math.exp(-$x/2) / Math.sqrt($x);\n $i1 = 1;\n } else {\n $p = $a = Math.exp(-$x/2);\n $i1 = 2;\n }\n\n for ($i = $i1; $i <= ($n-2); $i += 2) {\n $a *= $x / $i;\n $p += $a;\n }\n }\n return $p;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":56,"value":"OnlySearchResults = {\n var data = []\n for(var i = 0;i< search.length;i++){\n data.push(search[i][\"ID\"])\n }\n return data\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":63,"value":"matchingEntry = Object.values(table_data).find(entry => entry.ID === gene_id)","pinned":false,"hidden":false,"mode":"ojs"},{"id":812,"value":"selection_evidence = matchingEntry['Selection p-val'] <0.05 ? \"significant evidence for positive selection\" : \"no significant evidence for positive selection\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":711,"value":"stacked_omegas = {\n\n let hom_full = mod_results_json[gene_id][which_busted]['fits']['Unconstrained model']['Rate Distributions']['Test']\n \n const stacked = [\n ...Object.values(hom_full)\n ]\n let count = 1\n let hom = which_busted === 'BUSTED' ? stacked : stacked.slice(1)\n for (let i = 0; i < hom.length; i++) {\n hom[i]['dist'] = 'Homs';\n hom[i]['omega_bin'] = count;\n count++\n}\n return stacked\n //let omegas = which_busted === \"BUSTED\" ? stacked : stacked.slice(1)\n //return omegas \n\n//return test\n\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":583,"value":"mod_results_json = {\n let mod = {}\n for (let key in results_json){\n if (results_json[key]['FILTER'] && results_json[key]['FILTER'].hasOwnProperty('filter')){\n if (results_json[key]['BUSTED'] && results_json[key]['BUSTED'].hasOwnProperty('fits')){\n if (results_json[key]['BUSTED_E'] && results_json[key]['BUSTED_E'].hasOwnProperty('fits')){\n if (results_json[key]['aBSREL_UNFILTERED'] && results_json[key]['aBSREL_UNFILTERED'].hasOwnProperty('branch attributes')){\n mod[key] = results_json[key]\n }\n }\n }\n }\n }\n \nreturn mod\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":837,"value":"mod_results_json['NM_007118_5_23_species_23_seqs']['aBSREL_UNFILTERED']['branch attributes']['homo_sapiens']","pinned":false,"hidden":false,"mode":"ojs"},{"id":633,"value":"import {select,checkbox, autoSelect} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":718,"value":"import {Plot} from \"@mkfreeman/plot-tooltip\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":7,"value":"results_json = await d3.json (\"https://data.hyphy.org/web/busted-e/alycat_results.json\");","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W31","stdlib":"1","titling":"auto"}