![\"Imgur](\"https://i.imgur.com/1PCleoW.png\")
The difficulty with sequencing nucleic acids is nicely summarized by Hutchinson 2007:
\nIt is therefore not surprising that protein-sequencing was developed before DNA sequencing by Sanger and Tuppy 1951.
\ntRNA was the first complete nucleic acid sequenced (see pioneering work of Robert Holley and colleagues and also Holley’s Nobel Lecture). Conceptually, Holley’s approach was similar to Sanger’s protein sequencing: break molecule into small pieces with RNases, determine sequences of small fragments, and use overlaps between fragments to reconstruct (assemble) the final nucleotide sequence.
\nThe work on finding approaches to sequencing DNA molecules began in the late 60s and early 70s. One of the earliest contributions has been made by Ray Wu (Cornell) and Dave Kaiser (Stanford), who used E. coli DNA polymerase to incorporate radioactively labelled nucleotides into protruding ends of bacteriphage lambda. It took several more years for the development of more “high throughput” technologies by Sanger and Maxam/Gilbert. The Sanger technique has ultimately won over Maxam/Gilbert’s protocol due to its relative simplicity (once dideoxynucleotides has become commercially available) and the fact that it required a smaller amount of starting material as the polymerase was used to generate fragments necessary for sequence determination.
\n\n\nComment: Two Nobel prizes\nFred Sanger is one of only four people, who received two Nobel Prizes in their original form (for scientific, not societal, breakthroughs).
\n
This method builds on the idea of Wu and Kaiser (for minus part) and on the special property of DNA polymerase isolated from phage T4 (for plus part). The schematics of the method is given in the following figure (from Sanger & Coulson: 1975):
\n![\"Plus/minus](\"./images/sanger_plus_minus.png\")
Open image in new tabIn this method, a primer and DNA polymerase is used to synthesize DNA in the presence of P32-labeled nucleotides (only one of four is labeled). This generates P32-labeled copies of DNA being sequenced. These are then purified and (without denaturing) separated into two groups: minus and plus. Each group is further divided into four equal parts.
\nIn the case of minus polymerase a mix of nucleotides minus one is added to each of the four aliquotes: ACG (-T), ACT (-G), CGT (-A), AGT (-C). As a result in each case DNA strand is extended up to a missing nucleotide.
\nIn the case of plus only one nucleotide is added to each of the four aliquotes (+A, +C, +G, and +T) and T4 DNA polymerase is used. T4 DNA polymerase acts as an exonuclease that would degrade DNA\nfrom 3’-end up to a nucleotide that is supplied in the reaction.
\nThe products of these are loaded into a denaturing polyacrylamide gel as eight tracks (four for minus and four for plus; from Sanger & Coulson: 1975):
\n![\"Plus/minus](\"./images/sanger_plus_minus_gel.png\")
Open image in new tabIn this method DNA is terminally labeled with P32, separated into four equal aliquotes. Two of these are treated with Dimethyl sulfate (DMSO) and the remaining two are treated with hydrazine.
\nDMSO methylates G and A residues. Treatment of DMSO-incubated DNA with alkali at high temperature will break DNA chains at G and A with Gs being preferentially broken, while treatment of DMSO-incubated DNA with acid will preferentially break DNA at As. Likewise treating hydrazine-incubated DNA with piperidine breaks DNA at C and T, while DNA treated with hydrazine in the presence of NaCl preferentially breaks at Cs. The four reactions are then loaded on a gel generating the following picture (from Maxam & Gilbert: 1977):
\n![\"Maxam/Gilbert](\"./images/maxam_gilbert.jpeg\")
Open image in new tabThe original Sanger +/- method was not popular and had a number of technical limitations. In a new approach, Sanger took advantage of inhibitors that stop the extension of a DNA strand at particular nucleotides.\nThese inhibitors are dideoxy analogs of normal nucleotide triphosphates (from Sanger et al. 1977):
\n![\"Sanger](\"./images/sanger_ddN.png\")
Open image in new tabIn the original Sanger paper the authors sequenced bacteriophage phiX174 by using its own restriction fragments as primers. This was an ideal set up to show the proof of principle for the new method. This is because phiX174 DNA is homogeneous and can be isolated in large quantities. Now suppose that you would like to sequence a larger genome (say E. coli). Remember that the original version of Sanger method can only sequence fragments up to 200 nucleotides at a time. So to sequence the entire E. coli genome (which by-the-way was not sequenced until 1997) you would need to split the genome into multiple pieces and sequence each of them individually. This is hard because to produce a readable Sanger sequencing gel each sequence must be amplified to a suitable amount (around 1 nanogram) and be homogeneous (you cannot mix multiple DNA fragments in a single reaction as it will be impossible to interpret the gel). Molecular cloning enabled by the availability of commercially available restriction enzymes and cloning vectors simplified this process. Until the onset of next generation sequencing in 2005 the process for sequencing looked something like this:
\n![\"pGEM3z.](\"./images/pgem3z.png\")
Open image in new tabUntil the invention of NGS the above protocol was followed with some degree of automation. But you can see that it was quite laborious if the large number of fragments needed to be sequenced. This is because each of them needed to be subcloned and handled separately. This is in part why the Human Genome Project, which will be discussed in future lectures in detail, took so much time to complete.
\nThe simplest possible sequencing machine is a gel rig with polyacrylamide gel. Sanger used it is his protocol obtaining the following results (from Sanger et al. 1977):
\n![\"Sanger](\"./images/sanger_gel.png\")
Open image in new tabHere for sequencing each fragment four separate reactions are performed (with ddA, ddT, ggC, and ddG) and four lanes on the gel are used. One simplification of this process that\ncame in the 90s was to use fluorescently labeled dideoxy nucleotides. This is easier because everything can be performed in a single tube and uses a single lane on a gel\n(from Applied Biosystems support site):
\n![\"PA](\"./images/gel_versus_chromotogram.png\")
Open image in new tabHowever, there is still substantial labor involved in pouring the gels, loading them, running machines, and cleaning everything post-run. A significant improvement was offered by the development of capillary electrophoresis allowing automation of liquid handling and sample loading. Although several manufacturers have been developing and selling such machines a de facto standard in this area was (and still is) the Applied Biosystems (ABI) Genetics and DNA Analyzer systems. The highest throughput ABI system, 3730xl, had 96 capillaries and could automatically process 384 samples.
\n384 samples may sound like a lot, but it is nothing if we are sequencing an entire genome. The beauty of NGS is that these technologies are not bound by sample handling logistics. They still require the preparation of libraries,\nbut once a library is made (which can be automated) it is processed more or less automatically to generate multiple copies of each fragment\n(in the case of 454, Illumina, Ion Torrent, PacBio, Oxford Nanopore, Element, Complete Genomics etc…) and loaded onto the machine, where millions of individual fragments are sequenced simultaneously.\nWe will learn about these technologies later in this course.
\nIn a series of now classical papers (Paper 1, Paper2) Philip Green and co-workers have developed a quantitative framework for the analysis of data generated by automated DNA sequencers:
\n![\"phred1.](\"./images/phred1.png\")
Open image in new tab![\"phred2.](\"./images/phred2.png\")
Open image in new tabIn particular, they developed a standard metric for describing the reliability of base calls:
\nAn important technical aspect of our work is the use of log-transformed error probabilities rather than untransformed ones, which facilitates working with error rates in the range of most importance (very close to 0).\nSpecifically, we define the quality value \\(q\\) assigned to a base-call to be:
\n\n\\[q = -10\\times log\\_{10}(p)\\]\nwhere \\(p\\) is the estimated error probability for that base-call. Thus a base-call having a probability of 1/1000 of being incorrect is assigned a quality value of 30. Note that high-quality values correspond to low error probabilities, and conversely.
\nWe will be using the concept of “quality score” or “phred-scaled quality score” repeatedly in this course.
\nCan we sequence a genome using the shotgun approach?
\n![\"Myers/Green](\"./images/myers_versus_green.png\")
Open image in new tabIn this lesson we will cover some of the fundamental Python basics including variables, expressions, statements, and functions.
\n\n\nWarning: Indentation or bust!\nPython is an indented language: code blocks are defined using indentation with spaces!
\n
In Python, indentation is used to indicate the scope of control structures such as for loops, if statements, and function and class definitions. The amount of indentation is not fixed, but it must be consistent within a block of code.\nThe recommended amount of indentation is 4 spaces, although some developers prefer to use 2 spaces. Indenting is important in Python because it is used to indicate the level of nesting and structure of the code,\nwhich makes it easier to read and understand. Additionally, indentation is also used to indicate which lines of code are executed together as a block.
In this lecture, we will re-implement our Sanger sequencing simulator from the previous lecture and generate realistic gel images.
\nFirst, we import a module called random which contains a number of functions for generating and working with random numbers
Next, we will write a simple loop that would generate a sequence of pre-set lengths:
\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-3", "source": [ "seq = ''\n", "for _ in range(100):\n", " seq += random.choice('ATCG')" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-4", "source": "\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-5", "source": [ "seq" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-6", "source": "![\"xkcd1306.](\"https://imgs.xkcd.com/comics/sigil_cycle.png\")
CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTCGTCGCGCCCTCACAACTGCAAGGTCGTGGCACC\nThe code below iterates through each element of a sequence seq (assumed to be a string containing nucleotides) and it checks if the current nucleotide is equal to ‘A’. If it is, it generates a random number between 0 and 1 using the random.random() function.
It then checks if the random number is greater than 0.5. If it is, the code does nothing and proceeds to the next iteration. If the random number is less than or equal to 0.5, the code adds the\nlowercase version of the nucleotide (‘a’) to a string called synthesized_strand and then breaks out of the loop.
In every iteration of the loop, regardless of whether the nucleotide is ‘A’ or not, the code then adds the current nucleotide to the synthesized_strand string.
This means that when the current nucleotide is ‘A’, then the generated random number will decide whether the code will add the nucleotide ‘A’ or ‘a’ to the synthesized_strand, and it will break out of the loop\nafter adding the nucleotide to the synthesized_strand. To get a good idea of what is going on let’s visualize the code execution in
This can be simplified by first removing d_or_dd variable:
CTTGCGGCTATa\nand removing an unnecessary group of if ... else statements:
CTTGCGGCTATAGGAATa\nfinally let’s make synthesized_strand += 'a' a bit more generic:
CTTGCGGCTATAGGa
\nTo simulate 10 polymerase molecules we simply wrap the code from above into a for loop:
CTTGCGGCTa\n CTTGCGGCTa\n CTTGCGGCTa\n CTTGCGGCTa\n CTTGCGGCTa\n CTTGCGGCTa\n CTTGCGGCTATa\n CTTGCGGCTATAGGa\n CTTGCGGCTa\n CTTGCGGCTa
\nOne problem with this code is that does not actually save the newly synthesized strand: it simply prints it. To fix this we will create a list\n(or an array) called new_strands and initialize it by assigning an empty array to it:
[‘CTTGCGGCTa’,\n ‘CTTGCGGCTa’,\n ‘CTTGCGGCTATAGGAATAa’,\n ‘CTTGCGGCTATAGGAATa’,\n ‘CTTGCGGCTATa’,\n ‘CTTGCGGCTATAGGa’,\n ‘CTTGCGGCTATAGGa’,\n ‘CTTGCGGCTATa’,\n ‘CTTGCGGCTATAGGa’,\n ‘CTTGCGGCTa’]
\nAnd to repeat this for the remaining three nucleotides we will do the following crazy thing:
\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-21", "source": [ "new_strands = []\n", "for _ in range(10):\n", " synthesized_strand = ''\n", " for nucleotide in seq:\n", " if nucleotide == 'A' and random.random() > 0.5:\n", " synthesized_strand += nucleotide.lower()\n", " break\n", " synthesized_strand += nucleotide\n", " new_strands.append(synthesized_strand)\n", "\n", "for _ in range(10):\n", " synthesized_strand = ''\n", " for nucleotide in seq:\n", " if nucleotide == 'C' and random.random() > 0.5:\n", " synthesized_strand += nucleotide.lower()\n", " break\n", " synthesized_strand += nucleotide\n", " new_strands.append(synthesized_strand)\n", "\n", "for _ in range(10):\n", " synthesized_strand = ''\n", " for nucleotide in seq:\n", " if nucleotide == 'G' and random.random() > 0.5:\n", " synthesized_strand += nucleotide.lower()\n", " break\n", " synthesized_strand += nucleotide\n", " new_strands.append(synthesized_strand)\n", "\n", "for _ in range(10):\n", " synthesized_strand = ''\n", " for nucleotide in seq:\n", " if nucleotide == 'T' and random.random() > 0.5:\n", " synthesized_strand += nucleotide.lower()\n", " break\n", " synthesized_strand += nucleotide\n", " new_strands.append(synthesized_strand)" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-22", "source": "\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-23", "source": [ "len(new_strands)" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-24", "source": "40
\nRepeating the same code four times is just plain stupid so instead we will write a function called polymerase. Here we need to worry about the scope of variables. The scope of a variable refers to the regions\nof the code where the variable can be accessed or modified. Variables that are defined within a certain block of code (such as a function or a loop) are said to have a local scope, meaning that they can\nonly be accessed within that block of code. Variables that are defined outside of any block of code are said to have a global scope, meaning that they can be accessed from anywhere in the code.
In most programming languages, a variable defined within a function has a local scope, and it can only be accessed within that function. If a variable with the same name is defined outside the function,\nit will have a global scope and can be accessed from anywhere in the code. However, if a variable with the same name is defined within the function, it will take precedence over the global variable and will be used within the function.
\nThere are also some languages that have block scope, where a variable defined within a block (such as an if statement or a for loop) can only be accessed within that block and not outside of it.
\nIn Python, variables defined in the main module have a global scope and can be accessed from any function or module. Variables defined within a function have local scope, and they can only be\naccessed within that function. Variables defined within a loop or a block can be accessed only within the scope of the loop or block.
\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-25", "source": [ "def ddN(number_of_iterations, template, base, ddN_ratio):\n", " new_strands = []\n", " for _ in range(number_of_iterations):\n", " synthesized_strand = ''\n", " for nucleotide in template:\n", " if nucleotide == base and random.random() > ddN_ratio:\n", " synthesized_strand += nucleotide.lower()\n", " break\n", " synthesized_strand += nucleotide\n", " new_strands.append(synthesized_strand)\n", " return(new_strands)" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-26", "source": "\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-27", "source": [ "ddN(10,seq,'A',0.5)" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-28", "source": "[‘CTTGCGGCTATAGGa’,\n ‘CTTGCGGCTATa’,\n ‘CTTGCGGCTATa’,\n ‘CTTGCGGCTATAGGAATa’,\n ‘CTTGCGGCTa’,\n ‘CTTGCGGCTa’,\n ‘CTTGCGGCTa’,\n ‘CTTGCGGCTATa’,\n ‘CTTGCGGCTATAGGa’,\n ‘CTTGCGGCTATa’]
\nTo execute this function on all four types of ddNTPs with need to wrap it in a for loop iterating over the four possibilities:
To store the sequences being generated in the previous loop we will create and initialize a list called seq_run:
you will see that the seq run is a two-dimensional list:
\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-33", "source": [ "seq_run" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-34", "source": "[[‘CTTGCGGCTa’,\n ‘CTTGCGGCTATAGGa’,\n ‘CTTGCGGCTa’,\n ‘CTTGCGGCTATa’,\n ‘CTTGCGGCTATAGGa’,\n ‘CTTGCGGCTATa’,\n ‘CTTGCGGCTa’,\n ‘CTTGCGGCTATAGGAATAa’,\n ‘CTTGCGGCTa’,\n ‘CTTGCGGCTATa’],\n [‘CTTGCGGCt’, ‘Ct’, ‘Ct’, ‘Ct’, ‘Ct’, ‘CTt’, ‘Ct’, ‘Ct’, ‘CTTGCGGCt’, ‘Ct’],\n [‘CTTGc’,\n ‘c’,\n ‘c’,\n ‘c’,\n ‘c’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGc’,\n ‘c’,\n ‘CTTGCGGCTATAGGAATAAAAGGc’,\n ‘c’,\n ‘CTTGc’],\n [‘CTTg’,\n ‘CTTGCg’,\n ‘CTTg’,\n ‘CTTg’,\n ‘CTTGCg’,\n ‘CTTGCGGCTATAGGAATAAAAg’,\n ‘CTTGCGGCTATAGGAATAAAAg’,\n ‘CTTg’,\n ‘CTTg’,\n ‘CTTg’]]
\nas you will read in your next home assignment list elements can be addressed by “index”. The first element has the number 0:
\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-35", "source": [ "seq_run[0]" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-36", "source": "[‘CTTGCGGCTa’,\n ‘CTTGCGGCTATAGGa’,\n ‘CTTGCGGCTa’,\n ‘CTTGCGGCTATa’,\n ‘CTTGCGGCTATAGGa’,\n ‘CTTGCGGCTATa’,\n ‘CTTGCGGCTa’,\n ‘CTTGCGGCTATAGGAATAa’,\n ‘CTTGCGGCTa’,\n ‘CTTGCGGCTATa’]
\nAnother way to store these data is in a dictionary, which is a collection of key:value pairs where a key and value can be anything:
{‘A’: [‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTCGTCGCGCCCTCACa’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTa’,\n ‘CTTGCGGCTATAGGAATAa’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTa’,\n ‘CTTGCGGCTa’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTCGTCGCGCCCTCACAACTGCa’,\n ‘CTTGCGGCTATAGGa’,\n ‘CTTGCGGCTATAGGAATAAAa’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTCGTCGCGCCCTCACAACTGCAa’,\n ‘CTTGCGGCTATAGGAa’],\n ‘T’: [‘CTTGCGGCTAt’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTt’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTCGTCGCGCCCTCACAACt’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTCGTCGCGCCCTCACAACTGCAAGGTCGTGGCACC’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATt’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTAt’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCAt’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTCGTCGCGCCCTCACAACTGCAAGGTCGTGGCACC’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCAt’,\n ‘CTt’],\n ‘C’: [‘CTTGCGGCTATAGGAATAAAAGGc’,\n ‘CTTGCGGCTATAGGAATAAAAGGc’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTc’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGAc’,\n ‘c’,\n ‘CTTGCGGCTATAGGAATAAAAGGc’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTc’,\n ‘CTTGCGGc’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGAc’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGc’],\n ‘G’: [‘CTTGCGGCTATAGg’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTg’,\n ‘CTTGCg’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTCGTCGCGCCCTCACAACTGCAAGGTCGTGGCACC’,\n ‘CTTg’,\n ‘CTTGCGGCTATAGGAATAAAAGg’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTg’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCg’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGg’,\n ‘CTTGCGGCTATAGGAATAAAAg’]}
\ndictionary elements can be retrieved using a key:
\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-41", "source": [ "seq_run['A']" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-42", "source": "[‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTCGTCGCGCCCTCACa’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTa’,\n ‘CTTGCGGCTATAGGAATAa’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTa’,\n ‘CTTGCGGCTa’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTCGTCGCGCCCTCACAACTGCa’,\n ‘CTTGCGGCTATAGGa’,\n ‘CTTGCGGCTATAGGAATAAAa’,\n ‘CTTGCGGCTATAGGAATAAAAGGCTTTGCGGGTAGTGACCGCGCCGCGTATGTAATTCATGGGTGTCGTCGCGCCCTCACAACTGCAa’,\n ‘CTTGCGGCTATAGGAa’]
\nNow that we can simulate and store newly synthesized sequencing strands terminated with ddNTPs let us try to draw a realistic representation of the sequencing gel. For this, we will use\nseveral components that will be discussed in much greater detail in the upcoming lectures. These components are:
\npandas - a dataframe manipulation toolaltair - a statistical visualization library built on top of vega-light visualization grammarThese two libraries will be used in almost all lectures concerning Python in this class.
\nGel electophoresis separates molecules based on mass, shape, or charge. In the case of DNA all molecules are universally negatively charged and thus will always migrate\nto (+) electrode. All our molecules are linear single-stranded pieces (our gel is denaturing) and so the only physical/chemical characteristic that distinguishes them is length.\nTherefore the first thing we will do is to convert our sequences into their lengths. For this, we will initialize a new dictionary called seq_lengths:
{‘base’: [‘A’,\n ‘A’,\n ‘A’,\n ‘A’,\n ‘A’,\n ‘A’,\n ‘A’,\n ‘A’,\n ‘A’,\n ‘A’,\n ‘T’,\n ‘T’,\n ‘T’,\n ‘T’,\n ‘T’,\n ‘T’,\n ‘T’,\n ‘T’,\n ‘T’,\n ‘T’,\n ‘C’,\n ‘C’,\n ‘C’,\n ‘C’,\n ‘C’,\n ‘C’,\n ‘C’,\n ‘C’,\n ‘C’,\n ‘C’,\n ‘G’,\n ‘G’,\n ‘G’,\n ‘G’,\n ‘G’,\n ‘G’,\n ‘G’,\n ‘G’,\n ‘G’,\n ‘G’],\n ‘length’: [81,\n 54,\n 19,\n 34,\n 10,\n 87,\n 15,\n 21,\n 88,\n 16,\n 11,\n 27,\n 84,\n 100,\n 57,\n 51,\n 60,\n 100,\n 60,\n 3,\n 24,\n 24,\n 67,\n 39,\n 1,\n 24,\n 67,\n 8,\n 39,\n 47,\n 14,\n 65,\n 6,\n 100,\n 4,\n 23,\n 28,\n 43,\n 31,\n 22]}
\nnow let’s import pandas:
and inject seq_lengths into a pandas dataframe:
it looks pretty:
\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-51", "source": [ "sequences" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-52", "source": "| \n | base | \nlength | \n
|---|---|---|
| 0 | \nA | \n81 | \n
| 1 | \nA | \n54 | \n
| 2 | \nA | \n19 | \n
| 3 | \nA | \n34 | \n
| 4 | \nA | \n10 | \n
| 5 | \nA | \n87 | \n
| 6 | \nA | \n15 | \n
| 7 | \nA | \n21 | \n
| 8 | \nA | \n88 | \n
| 9 | \nA | \n16 | \n
| 10 | \nT | \n11 | \n
| 11 | \nT | \n27 | \n
| 12 | \nT | \n84 | \n
| 13 | \nT | \n100 | \n
| 14 | \nT | \n57 | \n
| 15 | \nT | \n51 | \n
| 16 | \nT | \n60 | \n
| 17 | \nT | \n100 | \n
| 18 | \nT | \n60 | \n
| 19 | \nT | \n3 | \n
| 20 | \nC | \n24 | \n
| 21 | \nC | \n24 | \n
| 22 | \nC | \n67 | \n
| 23 | \nC | \n39 | \n
| 24 | \nC | \n1 | \n
| 25 | \nC | \n24 | \n
| 26 | \nC | \n67 | \n
| 27 | \nC | \n8 | \n
| 28 | \nC | \n39 | \n
| 29 | \nC | \n47 | \n
| 30 | \nG | \n14 | \n
| 31 | \nG | \n65 | \n
| 32 | \nG | \n6 | \n
| 33 | \nG | \n100 | \n
| 34 | \nG | \n4 | \n
| 35 | \nG | \n23 | \n
| 36 | \nG | \n28 | \n
| 37 | \nG | \n43 | \n
| 38 | \nG | \n31 | \n
| 39 | \nG | \n22 | \n
In our data, there is a number of DNA fragments that have identical length (just look at the dataframe above). We can condense these by grouping dataframe entries first by nucleotide (['base']) and then by length (['length']). For each group we will then compute count and put it into a new column named, …, count:
| \n | base | \nlength | \ncount | \n
|---|---|---|---|
| 0 | \nA | \n10 | \n1 | \n
| 1 | \nA | \n15 | \n1 | \n
| 2 | \nA | \n16 | \n1 | \n
| 3 | \nA | \n19 | \n1 | \n
| 4 | \nA | \n21 | \n1 | \n
| 5 | \nA | \n34 | \n1 | \n
| 6 | \nA | \n54 | \n1 | \n
| 7 | \nA | \n81 | \n1 | \n
| 8 | \nA | \n87 | \n1 | \n
| 9 | \nA | \n88 | \n1 | \n
| 10 | \nC | \n1 | \n1 | \n
| 11 | \nC | \n8 | \n1 | \n
| 12 | \nC | \n24 | \n3 | \n
| 13 | \nC | \n39 | \n2 | \n
| 14 | \nC | \n47 | \n1 | \n
| 15 | \nC | \n67 | \n2 | \n
| 16 | \nG | \n4 | \n1 | \n
| 17 | \nG | \n6 | \n1 | \n
| 18 | \nG | \n14 | \n1 | \n
| 19 | \nG | \n22 | \n1 | \n
| 20 | \nG | \n23 | \n1 | \n
| 21 | \nG | \n28 | \n1 | \n
| 22 | \nG | \n31 | \n1 | \n
| 23 | \nG | \n43 | \n1 | \n
| 24 | \nG | \n65 | \n1 | \n
| 25 | \nG | \n100 | \n1 | \n
| 26 | \nT | \n3 | \n1 | \n
| 27 | \nT | \n11 | \n1 | \n
| 28 | \nT | \n27 | \n1 | \n
| 29 | \nT | \n51 | \n1 | \n
| 30 | \nT | \n57 | \n1 | \n
| 31 | \nT | \n60 | \n2 | \n
| 32 | \nT | \n84 | \n1 | \n
| 33 | \nT | \n100 | \n2 | \n
The following chart is created using the alt.Chart() function and passing the data as an argument. The mark_tick() function is used to create a tick chart with a thickness of 4 pixels.
The chart is encoded with two main axis:
\n'length' field of the data.'base' field of the data.\nThe chart also encodes a color, it encodes the 'count' field of the data, sets the legend to None, and uses the 'greys' scale from the Altair library.Finally, the chart properties are set to a width of 100 pixels and a height of 800 pixels.
\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-57", "source": [ "import altair as alt\n", "alt.Chart(sequences_grouped_by_length).mark_tick(thickness=4).encode(\n", " y = alt.Y('length:Q'),\n", " x = alt.X('base'),\n", " color=alt.Color('count:Q',legend=None,\n", " scale=alt.Scale(scheme=\"greys\"))\n", ").properties(\n", " width=100,\n", " height=800)" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-58", "source": "And here is a color version of the same graph using just one line of the gel:
\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-59", "source": [ "import altair as alt\n", "alt.Chart(sequences_grouped_by_length).mark_tick(thickness=4).encode(\n", " y = alt.Y('length:Q'),\n", " color=alt.Color('base:N',#legend=None,\n", " scale=alt.Scale(scheme=\"set1\"))\n", ").properties(\n", " width=20,\n", " height=800)" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-60", "source": "‘GTCGATGCCTGTTTGACCTAACTGGCGTGAAGGCTATATCAGTTATCCCAAGCGTAGGCTTTCAATTCGCCCGGTTGCGTCGCCCGATTATCAATCGCGGAAGGTGGGTGCGATTGGAAGTCCAAAACCTTTATCCTGACACACTTTCTGACTCGGCTTGGCAATGGGAAGTGTAGAACGTAGCGGGGACCTACATCATATCGTACATAACTGAGACGTGCTCACCCGCAGAGATAAGAACTGCAATACCCGGGTGAATACTTGGGGAGTCTCACCCAGATGGTTGGCCTGATCCTCCCC’
\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-65", "source": [ "# Function simulating a single run of a single polymerase molecule\n", "\n", "def ddN(number_of_iterations, template, base, ddN_ratio):\n", " new_strands = []\n", " for _ in range(number_of_iterations):\n", " synthesized_strand = ''\n", " for nucleotide in template:\n", " if nucleotide == base and random.random() > ddN_ratio:\n", " synthesized_strand += nucleotide.lower()\n", " break\n", " synthesized_strand += nucleotide\n", " new_strands.append(synthesized_strand)\n", " return(new_strands)" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-66", "source": "\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-67", "source": [ "# Generating simulated sequencing run\n", "\n", "seq_run = {}\n", "for nt in 'ATCG':\n", " seq_run[nt] = ddN(100000,seq,nt,0.95)" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-68", "source": "\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-69", "source": [ "# Computing lengths\n", "\n", "seq_lengths = {'base':[],'length':[]}\n", "for key in seq_run.keys():\n", " for sequence in seq_run[key]:\n", " seq_lengths['base'].append(key)\n", " seq_lengths['length'].append(len(sequence))" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-70", "source": "\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-71", "source": [ "# Converting dictionary into Pandas dataframe\n", "\n", "sequences = pd.DataFrame(seq_lengths)" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-72", "source": "\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-73", "source": [ "# Grouping by nucleotide and length\n", "\n", "sequences_grouped_by_length = sequences.groupby(\n", " ['base','length']\n", ").agg(\n", " count=pd.NamedAgg(\n", " column='length',\n", " aggfunc='count'\n", " )\n", ").reset_index()" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-74", "source": "\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-75", "source": [ "# Plotting (note the quadratic scale for realism)\n", "\n", "import altair as alt\n", "alt.Chart(sequences_grouped_by_length).mark_tick(thickness=4).encode(\n", " y = alt.Y('length:Q',scale=alt.Scale(type='sqrt')),\n", " x = alt.X('base'),\n", " color=alt.Color('count:Q',legend=None,\n", " scale=alt.Scale(type='log',scheme=\"greys\")),\n", " tooltip='count:Q'\n", ").properties(\n", " width=100,\n", " height=800)" ], "cell_type": "code", "execution_count": null, "outputs": [], "metadata": { "attributes": { "classes": [ "# The problem" ], "id": "" } } }, { "id": "cell-76", "source": "