• Computational Thinking
  • Decomposition
  • Pattern Recognition
  • Abstraction
  • Algorithm Design
• Not about memory
• About Problem solving

FASTA files show only one strand

We call it the forward strand

The other strand is not shown, but can be calculated

We call it the reverse strand

According to Watson & Crick, the reverse strand is complementary

• Change a by t and vice-versa
• Change c by g and vice-versa

Moreover, the reverse strand is reversed

• 5’ becomes 3’ and vice-versa
• i.e. we have to read it backwards

opposite_strand(dna)

Error in opposite_strand(dna): object 'dna' not found

We can use the debugger to see what is happening

dna <-c('t','a','a','c','g','t')
dna

[1] "t" "a" "a" "c" "g" "t"

opposite_strand(dna)

[1] "a" "c" "g" "t" "t" "a"

The relation between %G and %C is not constant over all the genome

In some parts of the genome there are more G than C

We can evaluate this using the GC skew \[\frac{G-C}{G+C}\]

When we examine a small part of the genome, it can happen that there are no G or C. Then, the previous program fails

This version will work better

nG <- sum(dna=="g")
nC <- sum(dna=="c")
if( (nG+nC) == 0) {
  ans <- 0
} else {
  ans <- (nG-nC)/(nG+nC)
}

Local statistics look only through a small window

A Window is a part of the genome, from an initial position, and with a fixed size

For example, we can analyze a region of 500 bp starting at position 3000 of the genome

position 3000, size 500

window <- dna[3000:3499]

In the general case, we can write

window <- dna[start:(start+size-1)]

Remember the function seq() for numeric sequences

window <- dna[seq(from=3000, length.out=500)]

In the general case, we can write

window <- dna[seq(from=start, length.out=size)]

seq(from, to, by, length.out)

All inputs are optional

Choose carefully which ones to use

seq(from=5, to=14)

 [1]  5  6  7  8  9 10 11 12 13 14

seq(from=5, length.out=10)

 [1]  5  6  7  8  9 10 11 12 13 14

seq(from=5, to=30, by=3)

[1]  5  8 11 14 17 20 23 26 29

We traverse the genome with several windows of fixed size

Each window starts after the end of the previous one

In other words, the windows do not overlap

We start at the first DNA position, until the last

Be careful to not step out of DNA

win_pos <- seq(from=1, to=length(dna)-size, by=size)

Now we apply gc_skew() to each element of win_pos

For example, for windows of size 100000

size <- 1E5
win_pos <- seq(from=1, to=length(dna)-size, by=size)
skew <- rep(NA, length(win_pos))
for(i  in 1:length(win_pos)) {
    skew[i] <- gc_skew(dna, win_pos[i], size)
}

The result depends on window’s size

• the GC skew is positive in one part of chromosome
• and negative on the other, on average
• it changes sign in two places

The reason for this asymmetry is the DNA replication process

The origin of replication is in one of the sign changing places

We want to write a program so we can find Ori on many different genomes

We can answer questions like:

• Are both replichores of the same size?
• If they are different, what is the biological reason?
• Is there some insertion?
• What are the extra genes? What do they do?

The values can change a lot between a window and the next one

skew data is noisy. A lot of small random changes

Since they have different values every day, we represent the income on day i by income[i]

We also represent the balance on day i by balance[i]

The vector income has the values for all days

The element income[i] has the value for a single day

We will see several “days”, represented by i

It is useful to think that

• i represents “today”
• i-1 represents “yesterday”

The balance on first day is your initial money

balance[1] <- income[1]

The balance of today is the balance of yesterday plus the income of today

for(i in 2:length(income)) {
    balance[i] <- balance[i-1] + income[i]
}

The GC skew in each window is like the income in each day

And we want to see the “balance”: cumulative skew

acc_skew[1] <- skew[1]
for(i in 2:length(skew)) {
    acc_skew[i] <- acc_skew[i-1] + skew[i]
}

Same pattern, same solution

How can we calculate acc_skew?

It is easy to see that

acc_skew[i]=acc_skew[i-1] + skew[i]

So we can write something like this:

acc_skew[1] <- skew[1]
for(i in 2:length(skew)) {
    acc_skew[i] <- acc_skew[i-1] + skew[i]
}

Can you write this as a function?

The maximum value of acc_skew is

max(acc_skew)

[1] 4.9877

It is not really important

We are looking for position, not value of maximum

This is the index of the maximum in the acc_skew vector

which.max(acc_skew)

[1] 155

This is always an integer number, between 1 and length(acc_skew)

It is not the position in the genome

Since acc_skew vector has the same size as win_pos, we can find the Ori position with this

win_pos[which.max(acc_skew)]

[1] 1540001

win_pos contains the genome positions of each window

This code gives us the position of the window with maximum GC-skew

First plot is GC-skew versus Position

Second plot is AT-skew versus Position

What is the relationship between GC-skew and AT-skew?

How can you show it?

• If we know the income, we can get balance
• If we know the skew, we can get acc_skew
• If we know the change, we can know the state

How many cells we have after 100 hours?

ncells[1] <- 1
growth[1] <- 0
for(i in 2:100) {
    growth[i] <- ncells[i-1]*0.1
    ncells[i] <- ncells[i-1] + growth[i]
}
ncells[100]

[1] 12527.83

In practice cells never grow that much

There is a limited amount of food

So we need to include food in the system

We represent it with the vector food

Let’s say that initially food[1] <- 20 [ng/uL]

More cells eat more food, limited by the eat_rate:
food_change[i] <- ncells[i-1]*food[i-1]*eat_rate

We assume that when a cell gets food, it divides:
growth[i] <- ncells[i-1]*food[i-1]*eat_rate

What happens now?

Make a program to simulate this system

What happens when food finish?

Is that realistic?