31:T9b4,md`# Testing recombination in the pandemic SARS-CoV-2 strains > *By Hongru Wang, Sergei L Kosakovsky Pond, Anton Nekrutenko, and Rasmus Nielsen* **Updated Jan 1st, 2021 to allow filtering by date range and country** To test whether there is recombination among the current pandemic SARS-CoV-2 strains, we look for linkage disequilibrium (LD) decay signal which is a hallmark of recombination. Over time recombination, when present, will disrupt linkage between SNPs, and this will happen faster for more distant SNPs. We collect full-length SARS-CoV-2 genomic sequences from gisaid-logo

and remove all genomes that contain more missing bases (N) per gene/peptide than a fixed threshold. For this analysis we retained **${N}** whole length sequences sampled between **${time_fmt(dateRange.range[0])} and ${time_fmt(dateRange.range[1])}** collected in **${geoLimiter}**. We performed codon alignment of all the genes on the virus genome except for several short regions that are not in functonal ORFs. We then called SNPs based on alignments. At minimim frequency of **${filter_threshold}**, there are **${variants_above_threshold.length}** SNPs. We next calculate the linkage disequilibrium (both r² and D\') between all pairs of SNPs. Plotting the pairwise linkage disequilibrium versus the distance (**Figure 1**), we do **not** observe a clear pattern of LD decay. A linear regression analysis to test if there is a negative relationship between linkage disequilibrium (r²) and distance between the SNPs, yields a regression coefficient of **${r2_regression.m.toPrecision (3)}** and the proportion of variance explained (R², distinct from r²) of **${r2_estimate.toPrecision (3)}**. To test if the fitted coefficient and R² represent a significant finding, we perform bootstrap tests: we randomly shuffle the coordinates of SNPs ${n_bootstrap} times, and then carry out the same regression test for each shuffling. We construct null distribution of fitted coefficient and R² based on the ${n_bootstrap} replicates. Comparing the fitted values with the null distribution (Figure 2), we conclude that there is no evidence of correlation between distance and LD for the pairwise comparison of SNPs across the **${N}** virus genome, thus there is evidence of no recombination in the current pandemic virus genomes. `32:Tb3e,vegalite({ "title": null, data: {values: pair_data}, hconcat: [{ layer: [{ width : 300, height: 300, "mark": { "type": "point", "filled" : true, clip : true, "tooltip": {"content": "data"} }, "encoding": { "size" : {"value" : 64}, "x": { field: "dist", axis : { "title" : "Distance, bp", "titleFontSize": 14, }, type: "quantitative", }, "y": {"field": "r2", "type": "quantitative", axis: { "title" : "r2", "titleFontSize": 14, "grid" : false, }, "scale": { "type": "linear", }}, } }, { "mark": { "type": "line", "color": "black" }, "transform": [ { "regression": "r2", "on": "dist", //"method" : "poly" } ], "encoding": { "x": { "field": "dist", "type": "quantitative" }, "y": { "field": "r2", "type": "quantitative" } } }, { "transform": [ { "regression": "r2", "on": "dist", //"method" : "poly", "params": true }, {"calculate": "'R²: '+format(datum.rSquared, '.5f')", "as": "R2"} ], "mark": { "type": "text", "color": "firebrick", "x": "width", "align": "right", "y": -5 }, "encoding": { "text": {"type": "nominal", "field": "R2"} } } ] }, { width : 300, height: 300, "mark": { "type": "point", "filled" : true, clip : true, "tooltip": {"content": "data"} }, "encoding": { "size" : {"value" : 64}, "x": { field: "dist", axis : { "title" : "Distance, bp", "titleFontSize": 14, }, type: "quantitative", }, "y": {"field": "Dprime", "type": "quantitative", axis: { "title" : "D'", "titleFontSize": 14, "grid" : false, }, "scale": { "type": "linear", }}, } } ] })33:T8ce,variant_histograms = vegalite({ "vconcat" : [{ "title": "Distribution of r2 under permutations vs the estimate", width : 800, height: 100, data: {values: bootstraps}, layer: [ { "mark": {"type": "bar"}, encoding: { "x": {field: "r2", type: "quantitative", axis: {"title": "Regression variance explained", titleFontSize : 14, grid : false, format: "2f"}, "bin": { "bin" :true, "maxbins" : 100 } }, "y": { axis: {"title": "Count", titleFontSize : 14, grid : false}, "aggregate": "count", "type": "quantitative" } } }, { "mark": "rule", "data" : { "values" : [{"r2" : r2_estimate}] }, "encoding": { "x": { "field": "r2", "type": "quantitative" }, "color": {"value": "red"}, "size": {"value": 3} } } ]}, { "title": "Distribution of the regression coefficient under permutations vs the estimate", width : 800, height: 100, data: {values: bootstraps}, layer: [ { "mark": {"type": "bar"}, encoding: { "x": {field: "slope", type: "quantitative", axis: {"title": "Regression coefficient", titleFontSize : 14, grid : false, format: "2f"}, "bin": { "bin" :true, "maxbins" : 100 } }, "y": { axis: {"title": "Count", titleFontSize : 14, grid : false}, "aggregate": "count", "type": "quantitative" } } }, { "mark": "rule", "data" : { "values" : [{"slope" : r2_regression.m}] }, "encoding": { "x": { "field": "slope", "type": "quantitative" }, "color": {"value": "red"}, "size": {"value": 3} } } ]} ] })34:T41e,table (_.filter (pair_data, (d)=>d.r2 >= 0.7), { limit: 5000, enableCSVDownload: true, sort: true, enableFilter: false, columns: [ { key: 'SNP1', name: 'SNP1' }, { key: 'SNP2', name: 'SNP2' }, { key: 'f1', name: 'Freq (SNP1)', render: (val)=>val.toFixed (3) }, { key: 'f2', name: 'Freq (SNP2)', render: (val)=>val.toFixed (3) }, { key: 'f12', name: 'Freq (Pair)', render: (val)=>val.toFixed (3) }, { key: 'dist', name: 'Distance, bp' }, { key: 'r2', name: 'LD r2', render: (val)=>val.toFixed (3) }, { key: 'Dprime', name: 'LD D\'', render: (val)=>val.toFixed (3) } ] });35:[] 2f:{"title":"Testing recombination in the pandemic SARS-CoV-2 strains","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":260,"value":"viewof filter_threshold = text({placeholder: \"0.03\", description: \"Minimim variant frequency to include in LD calculations; limit 100 variants\", submit: \"Update\", value : 0.03})","pinned":false,"hidden":false,"mode":"ojs"},{"id":452,"value":"viewof geoLimiter = autoSelect({\n options: _.keys(countries_and_codes),\n value: \"Entire World\",\n description: \"Only evaluate sequences from this country\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":542,"value":"md`Limit to sequences collected over a particular period of time`","pinned":false,"hidden":false,"mode":"ojs"},{"id":520,"value":"viewof dateRange = rangeSlider(timed_sequences, d=>d[0], value=>d3.sum(value.values.map(d=>d[1])), {yTicks:2,\n yScale:d3.scalePow().exponent(1)})","pinned":false,"hidden":false,"mode":"ojs"},{"id":295,"value":"md`**Figure 1** Linkage disequlibirum decay patterns`","pinned":false,"hidden":false,"mode":"ojs"},{"id":65,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":313,"value":"md`**Figure 2** Fitted regression coefficient and R² compared to random coordinate permutations`","pinned":false,"hidden":false,"mode":"ojs"},{"id":159,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":320,"value":"md`**Table 1** Pairs of SNPs in linkage disequilibrium (r² ≥ 0.7)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":227,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":29,"value":"md`---- `","pinned":false,"hidden":false,"mode":"ojs"},{"id":367,"value":"md`### Data and libararies`","pinned":false,"hidden":false,"mode":"ojs"},{"id":370,"value":"md`#### Data`","pinned":false,"hidden":false,"mode":"ojs"},{"id":20,"value":"sequence_variants = {\n let rd = await d3.json (\"https://raw.githubusercontent.com/veg/SARS-CoV-2/compact/data/fasta/2021-01-07/sequences.variants.json\");\n let seqDict = {};\n _.each (rd.sequences, (d, i)=> {\n seqDict[i] = {};\n for (let j = 0; j < d.length; j+=2) {\n seqDict[i][d[j]] = d[j+1];\n }\n });\n rd.sequences = seqDict;\n return rd;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":432,"value":"meta_data = await d3.json (\"https://raw.githubusercontent.com/veg/SARS-CoV-2/compact/data/db/sequence-info.json\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":377,"value":"md`#### Computed quantities`","pinned":false,"hidden":false,"mode":"ojs"},{"id":398,"value":"variants_included = {\n let geoCode = +countries_and_codes[geoLimiter];\n let dr = _.map (dateRange.range, (d)=>(d.getTime()-reference_date)/24/3600/1000);\n \n let seqs = _.chain (sequence_variants.sequences).toPairs().\n filter ((d)=>geoCode < 0 || meta_data.meta[d[0]][2] == geoCode).\n filter ((d)=> meta_data.values[0][meta_data.meta[d[0]][0]] >= dr[0] && meta_data.values[0][meta_data.meta[d[0]][0]] <= dr[1]).fromPairs().value();\n \n \n \n let variants = {};\n _.each (seqs, (d,i)=> {\n _.each (d, (n,p)=> {\n let tag = \"\" + p + n;\n if (tag in variants) {\n variants[tag] += 1;\n } else {\n variants[tag] = 1;\n }\n })\n });\n variants = _.map (variants, (c,d)=> [+d.substr (0, d.length - 1), d.substr (d.length - 1, d.length), c]);\n return {sequences: seqs, variants: variants};\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":406,"value":"variants_above_threshold = _.slice(_.sortBy (_.filter (variants_included.variants, (v)=>v[2]/N >= filter_threshold), (d)=>-d[2]), 0, 100)","pinned":false,"hidden":false,"mode":"ojs"},{"id":52,"value":"N = _.size (variants_included.sequences)","pinned":false,"hidden":false,"mode":"ojs"},{"id":51,"value":"pair_data = {\n //return common;\n let pairs = [];\n _.each (variants_above_threshold, (d, i)=>{\n let seqsIn = d[2] < N/2 ? _.filter (variants_included.sequences, (vt)=>vt[d[0]] == d[1]) : null;\n for (let j = i+1; j < variants_above_threshold.length; j++) {\n if (variants_above_threshold [j][0] != d[0]) {\n pairs.push (pair_frequency (d, variants_above_threshold [j], seqsIn));\n }\n }\n });\n \n return pairs;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":444,"value":"countries_and_codes = {\n return _.fromPairs(_.flatten ([[[\"Entire World\",\"-1\"]],_.chain (meta_data.values[2]).toPairs().filter ((d)=>d[1]).map ((d)=>[d[1],d[0]]).value()], true));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":189,"value":"r2_regression = ss.linearRegression(pair_data.map(d => ([d.dist, d.r2])))","pinned":false,"hidden":false,"mode":"ojs"},{"id":132,"value":"r2_estimate = ss.rSquared(pair_data.map(d => ([d.dist, d.r2])), model)","pinned":false,"hidden":false,"mode":"ojs"},{"id":155,"value":"n_bootstrap = 100","pinned":false,"hidden":false,"mode":"ojs"},{"id":129,"value":"model = ss.linearRegressionLine(r2_regression)","pinned":false,"hidden":false,"mode":"ojs"},{"id":490,"value":"reference_date = d3.timeParse (\"%Y%m%d\")(meta_data.reference).getTime()","pinned":false,"hidden":false,"mode":"ojs"},{"id":494,"value":"timed_sequences = _.chain (sequence_variants.sequences).map ((d,i)=>i).countBy ((d)=>meta_data.meta[+d][0]).toPairs().filter ((d)=>meta_data.values[0][+d[0]]).map ((d)=>[new Date (reference_date + meta_data.values[0][+d[0]]*24*3600*1000), d[1]]).value()","pinned":false,"hidden":false,"mode":"ojs"},{"id":374,"value":"md`#### Libraries, Functions, and Widgets`","pinned":false,"hidden":false,"mode":"ojs"},{"id":529,"value":"time_fmt = d3.timeFormat (\"%d %B %Y\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":9,"value":"d3 = require (\"d3\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":11,"value":"vegalite = require (\"vega-embed\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":16,"value":"_=require(\"lodash\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":32,"value":"pair_frequency = (snp1, snp2, seqs) => {\n \n let f1 = snp1[2]/N,\n f2 = snp2[2]/N;\n \n let f12 = _.filter (seqs || variants_included.sequences, (ls)=> {\n return ls[snp1[0]] == snp1[1] && ls[snp2[0]] == snp2[1];\n }).length/N;\n \n let D = f12-f1*f2;\n \n return {'f1' : f1,\n 'f2' :f2,\n 'SNP1' : \"\" + (snp1[0] + 1) + snp1[1],\n 'SNP2' : \"\" + (snp2[0] + 1) + snp2[1],\n 'f12': f12,\n 'r2' : D*D/(f1*(1-f1)*f2*(1-f2)),\n 'Dprime' : D / (D > 0 ? Math.min (f1*(1-f2),f2*(1-f1)) : Math.max (-f1*f2, -(1-f1)*(1-f2))),\n 'dist' : Math.abs (snp1[0]-snp2[0]),\n 'p1' : snp1[0] + 1,\n 'p2' : snp2[0] + 1,\n 'n1' : snp1[1],\n 'n2' : snp2[1]};\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":126,"value":"ss = require(\"simple-statistics@6.0.0/dist/simple-statistics.min.js\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":137,"value":"shuffle_coordinates = (input)=> {\n let coords = _.union (_.map (_.uniqBy (input, \"p1\"), (d)=> d[\"p1\"]), _.map (_.uniqBy (input, \"p2\"), (d)=> d[\"p2\"]));\n let map = {};\n _.each (_.shuffle (coords), (d,i)=> {\n map[coords[i]] = d;\n });\n let data = _.map (input, (d)=> {\n let dc = _.clone (d);\n dc.p1 = map[d.p1];\n dc.p2 = map[d.p2];\n dc.dist = Math.abs (dc.p1-dc.p2);\n return dc;\n });\n let lr = ss.linearRegression(data.map(d => ([d.dist, d.r2])));\n let lrm = ss.linearRegressionLine(lr);\n return [lr.m, ss.rSquared(data.map(d => ([d.dist, d.r2])), lrm)];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":140,"value":"bootstraps = _.map (_.range (n_bootstrap), (i)=>{let s = shuffle_coordinates (pair_data); return {'r2':s[1], 'slope' : s[0] }})","pinned":false,"hidden":false,"mode":"ojs"},{"id":224,"value":"import { table } from \"@gampleman/table\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":256,"value":"import {select, slider, checkbox, text, autoSelect} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":481,"value":"import {rangeSlider} from \"@bumbeishvili/utils\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}