A way to get absolute frequency

If we want to know how may times we rolled a 7, we can do

sum(deney == 7)

[1] 6658

Remember that:

• Comparisons produce logic vectors
• sum() of a logic vector is the number of TRUE

Relative frequency

An event is a logic vector.

sum(event) is a number between 0 and length(event)

The relative frequency is

sum(event)/length(event)

or, better

mean(event)

GC content of DNA fragment

What is the GC content of a random DNA fragment?

You can see that GC content is a random variable

dna <- sample(c("A", "C", "T", "G"), size, replace=TRUE)
gc_dna <- sum(gene=="C" | gene=="G")/size

size is the DNA fragment size. Symbol | means or

Simulate GC content

gc_of_fragment <- function(size) {
    dna <- sample(c("A", "C", "G", "T"), size=size, replace=TRUE)
    return(sum(dna=="C" | dna=="G")/size)
}
gc_several <- replicate(10000, gc_of_fragment(900))
hist(gc_several)

Histogram of relative frequencies

hist(gc_several, probability=TRUE)

Better GC content

We can have a better simulation if we use the real proportions for ("A","C","T","G")

Let’s evaluate them for Carsonella rudii

library(seqinr)
genome <- toupper(read.fasta("AP009180.fna")[[1]])
prob <- c(mean(genome=="A"), mean(genome=="C"),
          mean(genome=="G"), mean(genome=="T"))
prob

[1] 0.41797046 0.08455988 0.08108379 0.41638587

Simulate again

gc_of_fragment <- function(size, prob) {
    dna <- sample(c("A", "C", "G", "T"), size=size, prob=prob, replace=TRUE)
    return(sum(dna=="C" | dna=="G")/size)
}
gc_several <- replicate(10000, gc_of_fragment(900, prob))
hist(gc_several, probability=TRUE)

Range of GC content depends on fragment size

small_fragments <- replicate(10000, gc_of_fragment(200, prob))
medium_fragments <- replicate(10000, gc_of_fragment(600, prob))
big_fragments <- replicate(10000, gc_of_fragment(1200, prob))

Another view of GC content v/s size

What other Random Variables can you imagine?

• Average
• Number of trials until success
• Number plants that grow in the lab
• Number of people that answer the homework
• Number of people that passes the course

What other random variable do you see?

How many times will you do this course?

Depending on how many questions you ask and how many homework you do, you will have different chances of passing this course

To be neutral, let’s try with different probabilities

How many times do you need to do this course until you pass?

Counting until success

Assuming we have a function called fail(), that returns TRUE when the experiment fails; then this code gives us number of trials until success

k <- 1
while(fail()) {
    k <- k+1
}

Success & Failure: two sides of a coin

A simple way to model fail() is to use a coin

fail <- function() {
    coin <- sample(c("H","T"), size=1)
    return(coin != "H")
}

or, better

fail <- function() {
    return(sample(c(FALSE, TRUE), size=1))
}

A biased coin

This function takes one input: the probability of success

fail <- function(p) {
    return(sample(c(FALSE, TRUE), size=1, prob=c(p, 1-p)))
}

Tossing a coin until Heads

roll_until_head <- function(p) {
    i<-1
    while(fail(p)) {
        i <- i+1
    }
    return(i)
}

Results for 50%

sim <- replicate(1000, roll_until_head(0.5))
barplot(table(sim)/1000)

Results for 90%

sim <- replicate(1000, roll_until_head(0.9))
barplot(table(sim)/1000)

Results for 10%

sim <- replicate(1000, roll_until_head(0.1))
barplot(table(sim)/1000)