0:02
Now that we've implements native exact matching and use it to match artificial
reads against the genome, let's do the same thing with real sequencing reads.
So we've already downloaded the 5X genome.
The other link we've given you is to a set of actual
sequencing reads from that same genome.
So you can download that, and I'm going to do that with Wget.
1:14
As this has downloaded the sequencing reads from that file and
that's a set of a 1000 reads.
So, now let's do the same thing we did above with artificial reads
in practical six and actually match these, align these against the genome and
see how many of them match.
1:53
And n is just a count of the total number of reads that we've processed.
So now I can say, if the length of the set of matches returned is greater than zero,
that means that they read aligned in at least one place, so we found a match.
So I can increment numMatched.
2:34
And only seven out of 1,000 reads match the genome.
>> Okay, so that might be surprising.
These are sequencing reads actually from this genome, and yet
only 7 out of 1,000 of them actually match the genome.
So why would that be?
Well, we actually learned a couple of reasons why.
A sequence obtained with a DNA sequencer might not match
back to the genome exactly.
One reason is because of sequencing errors, so
it could be that some of the bases were read incorrectly by the sequencer.
And then another reason is because there might be differences between
the individual that was sequenced,
in this case it's a virus, and the reference genome sequence.
So in this case it's probably mostly the sequencing errors that are the reason why
many of these reads are not aligning.
So one of the things that we're going to try next is to just take
part of the read, right?
It could be that the whole read, all 100 bases, don't necessarily match.
But if we just take a piece, like maybe the first 30 bases or so.
That has a greater chance of matching exactly.
All right, so we'll try that next.
>> Okay. So, I just can make a little change to get
the first 30 bases.
We'll set r equal to, just the 30 base prefix of r.
I'm going to run it again.
3:54
And this time, we're 459 out of 1000 reads match.
So, a little under one half.
>> So this still isn't very close to a thousand.
This is less than half of the reads aligning to our reference genome.
So you might be able to think of another reason why maybe not so
many of the reads are aligning.
In this case its because the genome is double stranded and
the reads can come from one strand or they can come from the other.
But the exact matching problem we've set up is only looking on
one of those two strands.
It's not looking on both.
So that could explain in large part why we're seeing
less than half of the reads actually matching back to the genome.
So the next thing to try would be not just to try to match the read to the genome,
but also to try to match the reverse complement of the read to the genome.
That will cover the case where the read comes from the other strand of the genome.
5:07
But I'm going to do two steps here.
So I'm setting matches to be the naive exact match between r and the genome.
But I'm going to do something similar, I'm going to say matches.extend
with naive of the reverse compliment
of r in the genome.
So what dot extend means is it just takes that list and
adds on to it the results of another list.
So the first line here, matches equals naive of r in the genome
is going to test whether that read matches in the forward direction on the genome.
In the next line, we'll add onto it the results
of any matches of that read to the reverse complement of the genome.
6:00
>> Oh [LAUGH].
So the problem here is that our.
Our dictionary in the reverse complement function doesn't allow for
there to potentially be an n in a sequencing read.
>> Right. >> And yet, like we said before,
sometimes a sequencing read has an n meaning no competence.
So we just have to add one more key value pair to our dictionary there.
>> So I'll just add the base n and this will just match to an n.
Since in both cases we don't know what it is.
>> Yeah. >> So I mean update that and run it again.
6:40
>> It's taking a little longer because it has to do more matching now.
>> There we go 932 out 1,000 reads have matched.
>> Yeah. So still not 100% But
this is to be expected, and because of sequencing errors and things like that.
But once we figured out that we had missed one whole strand of the genome, and
once we took a prefix of the read instead of the entire read,
we were able to match nearly all the reads in our data set.
But the fact that we had to take the first 30 base pairs or so
is a pretty good indication that exact matching
is not necessarily going to be what we want to do, right?
We want to allow for sequencing errors,
we want to allow for differences between individuals.
So that's why later in the course we'll
examine some approximate matching algorithms.
>> Yeah.
Ben Langmead, PhDAssistant Professor
Jacob Pritt
