• The goal of Science is to understand Nature
• We find the laws that explain Patterns in Nature
• The first step is to find Patterns
  • Day & Night, Winter & Summer are Patterns
  • Understanding seasons allows us to survive
  • Some food makes us healthy
  • Understanding which food makes us sick is essential to survive

When we see a pattern we can make

• loops, if we need to repeat the pattern in one place
  • Like drawing each angle in a star
• functions, if the pattern is used in several places
  • Like drawing several stars in different places.

There is another kind of pattern

One useful mathematical function is factorial

The factorial of a \(n\) is \(n!\), defined as \[n! = n\cdot(n-1)!\]

That is, the factorial of n is n times the factorial of n-1

How do we write it in R?

factorial <- function(n) {
    ans <- n * factorial(n-1)
    return(ans)
}

There is an error here. Can you find it?

Use the debugger to see what is happening

Each recursive function needs an exit condition

That is, some rule to decide when to finish

There is always an easy case that does not require recursion

For example, the factorial of 1 is \(1! = 1\).

factorial <- function(n) {
    if(n <= 1) {
        ans <- 1
    } else {
        ans <- n * factorial(n-1)
    }
    return(ans)
}

if(cond) {
    code if true
}

if(cond) {
    code if true
} else {
    code if false
}

library(TurtleGraphics)
turtle_init(400, 400, mode="clip")
turtle_setpos(100, 250)
turtle_setangle(60)
for(i in 1:40) {
  turtle_move(200)
  turtle_right(10)
  if(i>20) {
      turtle_col("blue")
  }
}

• use if() to make decisions
  • use if(cond){A} to do A only when cond is TRUE
  • use if(cond){A} else {B} to do either A or B but not both
• a recursive function is used inside itself
  • there must be an exit condition
• We have all the “lego” pieces we will need
  • loops
  • decisions
  • functions

Despite being large molecules, DNA and proteins are made with pieces of only a few types

DNA is made of four bases. We can represent it with four letters

AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC

Proteins have 20 aminoacids and three stop codons. We can represent them with 23 letters.

MKRISTTITTTITITTGNGAG

In molecular biology we often work with sequences

• DNA sequences use 4 letters to represent the nucleotides in one of the two strands
• Protein sequences use 20 letters to represent the amino-acids, from amino to carboxyl terminal
• Other sequences are sometime used:
  • RNA,
  • DNA with ambiguous nucleotides,
  • amino-acid sequences with stop codons

The main reason why computing is useful for molecular biology

• DNA is discrete data
• Either “A”, “C”, “G” or “T”
• No middle-values
• All other things we measure are in a range
  • Like temperature, concentration, expression
• We can forget all the chemistry and work with symbols

• For each nucleotide
  • if the nucleotide is A, turtle moves left
  • if the nucleotide is T, turtle moves right
  • if the nucleotide is G, turtle moves up
  • if the nucleotide is C, turtle moves down

There are several ways to store DNA or protein data

Most of the times they are stored in FASTA format

FASTA files are text files, with some rules

You should never use word to store sequences

NCBI stores all public biological sequences at
https://www.ncbi.nlm.nih.gov/nuccore

Anybody can upload sequences, and they may be wrong

NCBI has a curation process to validate the sequences

If a sequence is good enough to be a reference, then it is stored in the RefSeq collection https://www.ncbi.nlm.nih.gov/refseq

Accession numbers are the best way to identify a biological sequence

• Candidatus Carsonella ruddii PV DNA has accession AP009180
• Escherichia coli str. K-12 substr. MG1655 has accession NC_000913

Different sequences can have the same name, but never the same accession

NCBI search box looks for patterns everywhere, not only on the title

and many sequences can have the same name

To be sure of using the correct sequence, go to NCBI Taxonomy

https://www.ncbi.nlm.nih.gov/taxonomy

Try with “Escherichia coli”. There are many

To handle sequence data in R, we use the seqinr library

You have to install it once.

install.packages("seqinr")

Then you have to load it on every session

library(seqinr)

read.fasta(file, seqtype = "DNA", ...)

Inputs

file

The name of the file which the sequences in FASTA format are to be read from

seqtype

the nature of the sequence: DNA or AA, defaulting to DNA

Output

A list of vectors of chars. Each element is a sequence object.

• Some organisms have been sequenced
• We can read the sequences from FASTA files
• read.fasta returns a list of vectors of characters
• We can create functions to do things with the sequence
• Functions are black boxes with inputs and outputs
• For example we do statistics of the sequence

Statistics is a way to tell a story that makes sense of the data

In genomics, we look for biological sense

That story can be global: about the complete genome

Or can be local: about some region of the genome

We will start with global properties

Measuring the melting temperature of the DNA double helix using spectrophotometry

• The absorbance of DNA at a wavelength of 260nm increases when DNA separates into two single strands at melting temperature

If the DNA has been sequenced then the GC-content can be accurately calculated by simple arithmetic.

GC-content percentage is calculated as \[\frac{G+C}{A+T+G+C}\]

Determine the GC content of E.coli

Is this GC content uniform through all genome?

The result should depend on:

• the genomic sequence
• the position of the window
• the size of the window

gc_content <- function(sequence, position, size) {
    # count letters
    return(ans)
}

We want to evaluate gc_content on different positions of the genome

positions <- seq(from=1, to=length(dna), by= window_size)

Now we apply gc_content to each element of positions?

for(pos in seq(from=1, to=length(dna), by= window_size)) {
    skew[pos] <- gc_skew(pos, dna, window_size)
}

What is the difference \(G-C\) respect to the total \(G+C\)?

Calculate the ratio \[\frac{G-C}{G+C}\] using sliding windows and plot it

What do you see?

There is a richness of G over C and T over A in the leading strand

There is a richness of C over G and A over T in the lagging strand

GC skew changes sign at the boundaries of the two replichores

This corresponds to DNA replication origin or terminus

GC skew changes sign at the boundaries of the two replichores

This corresponds to DNA replication origin or terminus

The replication origin is usually called ori.

If you want to learn more about the science, read:

Roten C-AH, Gamba P, Barblan J-L, Karamata D. Comparative Genometrics (CG): a database dedicated to biometric comparisons of whole genomes. Nucleic Acids Research. 2002;30(1):142-144