31:T40d,consensusSeq = (set)=> { let diff = 0; let consensus = _.map (_.range (set[0].length), (d)=>new Object()); /*_.each (consensus, (d,i)=> { consensus[i] = _.maxBy (_.map (_.countBy (_.map (set, (s)=>s[i])), (n,i)=>[i,n]), (d)=>d[1])[0]; if (consensus[i] != referenceSequence[i]) { diff++; } });*/ for (let s = 0; s < set.length; s++) { for (let n = 0; n < set[0].length; n++) { let l = set[s][n]; if (consensus[n][l]) { consensus[n][l] ++; } else { consensus[n][l] = 1; } } } for (let s = 0; s < consensus.length; s++) { let max = 0, maxi = 0; for (let n = 0; n < alphabet.length; n++) { if ( consensus[s][alphabet[n]] > max) { max = consensus[s][alphabet[n]]; maxi = alphabet[n]; } } consensus[s] = maxi; if (maxi != referenceSequence[s]) diff++; } if (i==1) { mutable counter = 0; } if (diff == 0) { mutable counter ++; } return [consensus.join(""), diff]; }32:[] 2f:{"title":"Resolving sequencing error by consensus","cells":[{"id":0,"value":"md`# Resolving sequencing error by consensus\n\nThis is a simple illustration how increasing the depth of **coverage** and relying on the **random and independent** distrbution of sequencing errors can help recover the correct sequence even when individual reads are quite noisy\n\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":154,"value":"viewof i = Scrubber(replicates, {autoplay: false,delay: 50, loop:false})","pinned":false,"hidden":false,"mode":"ojs"},{"id":216,"value":"viewof coverage = slider({\n min: 1, \n max: 100, \n step: 1, \n value: 10, \n title: \"Sequencing coverage\", \n //description: \"Integers from zero through 100\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":217,"value":"viewof rate = slider({\n min: 0.0,\n max: 0.5,\n value: 0.05,\n precision: 3,\n title: \"Error rate\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":45,"value":"md`\nCorrectly sequenced: ${counter}/${(i)}\n\n

\n\n\nTrue seq.: ${referenceSequence}
\nConsensus: ${highlightDifferences(referenceSequence,consensusSequence[0],\"red\",1)}\n
Seq. read:
Seq. read:\n${_.map (simulatedReads, (r)=>highlightDifferences (referenceSequence,r)).join (\"
Seq. read: \")}\n\n\n\n

\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":6,"value":"md`---- \nLibraries and data from this point on`","pinned":false,"hidden":false,"mode":"ojs"},{"id":25,"value":"N = 75 ","pinned":false,"hidden":false,"mode":"ojs"},{"id":156,"value":"replicates = _.range (1,101)","pinned":false,"hidden":false,"mode":"ojs"},{"id":28,"value":"alphabet = \"ACGT\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":64,"value":"mutants = {\n let r = {};\n _.each (alphabet, (l,i) => {\n r[l] = [];\n _.each (alphabet, (l2, i2)=> {\n if (i2 != i) {\n r[l].push (l2);\n }\n }); \n });\n return r;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":32,"value":"randomNuc = d3.randomInt(alphabet.length)","pinned":false,"hidden":false,"mode":"ojs"},{"id":55,"value":"randomNucM1 = d3.randomInt (alphabet.length-1)","pinned":false,"hidden":false,"mode":"ojs"},{"id":31,"value":"generateSequence = (alph, length) => {\n return _.map (_.range (length), (d)=>alphabet [randomNuc()]).join (\"\")\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":38,"value":"referenceSequence = generateSequence (alphabet, N)","pinned":false,"hidden":false,"mode":"ojs"},{"id":50,"value":"mutateSequence = (seq, rate) => {\n let mutSeq = new Array (seq.length);\n _.each (seq, (n,i)=> {\n if (Math.random() < rate) {\n mutSeq[i] = _.sample(mutants[n]);\n } else {\n mutSeq[i] = n; \n }\n });\n return mutSeq.join (\"\");\n};","pinned":false,"hidden":false,"mode":"ojs"},{"id":95,"value":"simulatedReads = {\n i;\n return _.map (_.range (coverage), (d)=>mutateSequence (referenceSequence, rate));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":124,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":199,"value":"mutable counter = 0","pinned":false,"hidden":false,"mode":"ojs"},{"id":165,"value":"consensusSequence = consensusSeq (simulatedReads)","pinned":false,"hidden":false,"mode":"ojs"},{"id":78,"value":"highlightDifferences = (ref, qry,diff_color,opacity)=> {\n diff_color = diff_color || \"black\";\n opacity = opacity || 0.2;\n let r = new Array (ref.length);\n _.each (ref, (l,i)=> {\n if (l != qry [i]) {\n r[i] = \"\" + qry[i] + \"\"\n } else {\n r[i] = \"\" + qry[i] + \"\"; \n }\n });\n return r.join (\"\");\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":10,"value":"_ = require(\"lodash\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":13,"value":"d3 = require (\"d3\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":16,"value":"import { Scrubber } from \"@mbostock/scrubber\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":19,"value":"import {select, slider, checkbox, text, autoSelect} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W32","stdlib":"1","titling":"auto"}