Practical: Working with sequencing reads

Loading...

Do curso por Johns Hopkins University

Algoritmos para Sequenciamento de DNA

331 classificações

Johns Hopkins University

Algoritmos para Sequenciamento de DNA

331 classificações

Curso 4 de 8 em Specialization Genomic Data Science

We will learn computational methods -- algorithms and data structures -- for analyzing DNA sequencing data. We will learn a little about DNA, genomics, and how DNA sequencing is used. We will use Python to implement key algorithms and data structures and to analyze real genomes and DNA sequencing datasets.

Na lição

DNA sequencing, strings and matching

This module we begin our exploration of algorithms for analyzing DNA sequencing data. We'll discuss DNA sequencing technology, its past and present, and how it works.

Module 1 Introduction1:32

Lecture: Why study this?4:29

Lecture: DNA sequencing past and present3:08

Lecture: Genomes as strings, reads as substrings5:02

Lecture: String definitions and Python examples3:33

Practical: String basics 7:14

Practical: Manipulating DNA strings 7:00

Practical: Downloading and parsing a genome 6:29

Lecture: How DNA gets copied3:37

Optional lecture: How second-generation sequencers work 7:10

Optional lecture: Sequencing errors and base qualities 6:37

Lecture: Sequencing reads in FASTQ format4:39

Practical: Working with sequencing reads 11:11

Practical: Analyzing reads by position 6:58

Lecture: Sequencers give pieces to genomic puzzles5:41

Lecture: Read alignment and why it's hard3:09

Lecture: Naive exact matching10:18

Practical: Matching artificial reads 6:53

Practical: Matching real reads 7:24

Conheça os instrutores

Ben Langmead, PhD
Assistant Professor
Computer Science
Jacob Pritt

Department of Computer Science

In this practical, we'll be describing how to download and

work with a set of sequencing reads.

So, we've given you a URL to a data set.

If you download that file and open it, you should see something like this.

This is a fast q file, which is a common format used for

storing reads from a sequencing experiment.

Now if you'll look closely at these lines, you'll see that they come in sets of four.

The first line in each set of four starts with an @ symbol, and

then just has a name.

This is just a tag to indicate which read this is.

The following line is a string of DNA bases.

The next line is just a plus sign,

and then the third line is just a bunch of symbols.

This is the quality sequence from that read.

So if you remember from lecture, each quality score is Phred 33 encoded.

This converts into a ASCII character.

And that's what's printed in the fourth line here.

>> So this file is sequencing reads from a real human being.

So these are actual human sequencing reads we're looking at.

>> Okay, now let's download this dataset in Python.

So, once again I'm going to use wget.

I'm going to paste in the URL, now our data set has been downloaded.

Let's write a function to read this *.fastq file and parse it.

So I'm going to just copy over a set of lines here just to refer to

when writing this function.

So I'm going to create a new, new function, reed fast Q,

pass in a file name, and we want to store a list of sequences from this file,

and we want to store a list of quality string associated with each sequence.

So I'm going to say with, I'm going to open the file.

I give it a file handle.

And then, I'm just going to put a while loop, I'm going to say while True, so

we'll loop forever, but we'll break when we reach the end of the file.

So the first thing I need to do is, so each time I go through this loop,

I'm going to need four lines.

So the first line, if you look back at this block here.

the first line is just the tag for the read.

We don't care what's in this line, so we're just going to read it, but

we're not going to store that to anything.

The second line is the sequence, so

I'm going to say, Save this to the variable seek.

The third line is just a plus sign.

We don't care about that.

So once again, we're not going to store that.

And then, the last line is the string of quality scores, so

we're going to save that to qua.

And you want to strip white space again here, right?

So the seek and qual lines are going to have new lines at the end of them, so

we want to get rid of those.

>> Right, yeah.

So I need to add an rstrip to the end of each line.

And now that I've read these lines, I have to check if I reached the end of the file.

So if I reach the end of the file at a previous state,

all of these strings will just be initialized to an empty string.

So I want to check.

I'm going to make sure that the length of the sequence is greater than zero.

So if that sequence has a length zero, then we've reached the end of the file and

there's nothing left.

So in this case, I'm going to break out of my while loop.

If that's not the case, if we're still in the file,

still reading lines, then I want to append these values to my list.

So I'm going to do sequences.append(seq) the sequence that I read, and

I'm going to do a similar thing for the string of qualities course that I read.

And then when I'm done, I'm just going to return the sequences and quality strings.

And if you're new to Python, you can return multiple values from one function

just like this, by separating them with commas.

Okay, now we've written this function,

let's run it on the data set we downloaded.

Worded.

So I'm going to save these to variables for sequences and qualities.

going to run readFastq.

And I need to paste in there, the name of the file that we downloaded.

So now you have two parallel lists, one that has all the sequences in it, and

one that has all the quality strings in it.

>> Right. So let's make sure that these are filled

out correctly.

Let me print the first five sequences.

There you go.

And then we'll do something similar.

Let's just print out the first five quality strings.

And there are the strings of quality scores associated with each read.

Now that we've downloaded and parsed this data set, let's do some analysis of it.

Let's create a histogram of the different quality scores,

to see what quality scores are most common and which are least common.

So before I do this, I want to write a quick helper function

to convert the ASCIIi symbols to the quality score.

So this will convert phred 33 encoded value to just a quality score.

And we'll pass in a quality value.

The function word just takes a character and

converts it to it's asking value in integer number.

So I'm going to convert that ASCII value to a character,

or to a number, and then subtract 33, and then return that.

So I see a bunch of hash symbols in those quality strings.

What happens when you run that function on a hash symbol?

>> Let's try that.

So I'm going to do this, put in a hash symbol.

And that's a quality score of two.

>> That's low.

>> Right. A low quality score,

remember, means that we have a very low confidence in our value.

It's not that likely to be correct.

In this case, I think a quality score of two, corresponds to about a 30,

maybe a 32% chance that we're incorrect.

So, you notice, most of these strings of hashtags occur near the end of the reads.

And as was discussed in lecture, as you step down the reading and

get towards the end, your confidence gets lower.

>> So try now, I see also some capital Js in there.

Just try a capital J.

>> Okay. Oh.

>> Oopsy.

>> So J is a score of 41, so that's very high.

That's a very high confidence in our sent value.

>> Yeah, that's, less than one is.

In what, 10,000 chance that the base call was incorrect.

>> Yeah.

Okay.

Now let's write our function to create a histogram of the quality scores.

So I'm going to write a function, createHist.

We'll pass in a list of quality strings.

I'm going to create a histogram.

This will just keep track of the frequency of each quality score that we've seen.

I think the highest quality score you can see is actually 41, which we saw.

Just to be safe,

I'm going to go up to 50, so that should give us some extra padding on the end.

And I'm going to say for each string of quality values,

qual, in our list qualities.

And then, in each string I want to step through and look at each quality score.

So for each phred quality score in that string.

So, firstly, I need to convert that to a number.

So I'll say q equals.

And I'll call the function I wrote above, phred33ToQ.

So now I have the number associated with, or the actuality quality score for

that position, so I'm going to increment that value in the histogram.

And then when I'm done, I'll just return the histogram.

So let's run this on our quality scores that we downloaded.

And print out the histogram.

So you can see the lowest quality score we have is a value of two.

That's actually the lowest that can be returned.

And the highest is, this should correspond to 41, I think for this data set.

But it's really hard to visualize this, when it's just printed out in text format.

So, let's graph it.

We're going to do this using Nat Plot Lib.

>> So that's just a special command that you can use in I Python,

that tells Nat Plot Lib to actually put the plots right there in the same

document with the text that you're typing.

>> Right ,and this line we'll just import >> The module that we need.

We'll call it matpotlib.pyplot.

And I'm going to plot this as a bar chart.

And I need to pass in the X values and the Y values.

So the Y values are going to be the values of our histogram.

The X values are just going to be the numbers from zero

up to the length of the histogram.

>> So I'm going to use range up to the length of h,

and then the y values will just be h.

And then to actually display it, I need to call plt.show.

And here's the result.

So you can see we have kind of a big spike of low quality value,

then it drops down really low.

And then it increases gradually and we have a lot of high quality values.

So Ben, why is there such a big spike at two,

but doesn't seem like there's any quality scores at three?

>> Yeah, well I'm guessing that what's happening at a quality value of two,

is that these are clusters where the base caller was.

>> Not confident at all at what the base ought to be there.

For example, if it was an equal mix of colors coming from that cluster,

it will have very low confidence.

And so that probability that the call is incorrect, it's going to be high.

It's going to be like a half ,or something like that, and

that corresponds to a very low Q value.

So that corresponds to maybe two or so.

Because we see that fairly often toward the ends of the reads, for

reasons that we discussed in the previous optional lecture,

we have a lot of occurrences with a quality value of two.

>> Cool.

Explore nosso catálogo

Registre-se gratuitamente e obtenha recomendações, atualizações e ofertas personalizadas.

O Coursera proporciona acesso universal à melhor educação do mundo fazendo parcerias com as melhores universidades e organizações para oferecer cursos on-line.

© 2018 Coursera Inc. Todos os direitos reservados.

Baixar na App Store

Baixar no Google Play