Computational Biology Blog in fasta format

If the typing monkeys have met Mr Markov: probabilities of spelling "omglolbbq" after the digitial monkeys have read Dracula

2014-09-15T13:54:00.001-05:00

On the weekend, randomly after watching Catching Fire, I remember the problem of the typing monkeys (Infinite monkey theorem) in which basically could be defined as (Thanks to Wiki):

# *******************
# INTRODUCTION
# *******************

The infinite monkey theorem states that a monkey hitting keys at random on a typewriter keyboard for an infinite amount of time will almost surely type a given text, such as the complete works of William Shakespeare.

There is a straightforward proof of this theorem. As an introduction, recall that if two events are statistically independent, then the probability of both happening equals the product of the probabilities of each one happening independently. For example, if the chance of rain in Moscow on a particular day in the future is 0.4 and the chance of an earthquake in San Francisco on that same day is 0.00003, then the chance of both happening on that day is 0.4 * 0.00003 = 0.000012, assuming that they are indeed independent.

Suppose the typewriter has 50 keys, and the word to be typed is banana. If the keys are pressed randomly and independently, it means that each key has an equal chance of being pressed. Then, the chance that the first letter typed is 'b' is 1/50, and the chance that the second letter typed is a is also 1/50, and so on. Therefore, the chance of the first six letters spelling banana is

less than one in 15 billion, but not zero, hence a possible outcome.

# *******************
# METHODS
# *******************

In my implementation, I will only consider 26 characters of the alphabet (from a to z, excluding the whitespace). The real question I would like to ask is the following:

Given a target word, say "banana", how many monkeys would be needed to have at least one successful event (a monkey typed the target) after the monkeys have typed 6 characters.

To solve this, first calculate the probability of typing the word banana:

Now, just compute the number of monkeys that might be needed:

The model that assigns the same probability for each character is labeled as "uniform model" in my simulation.

My goal is to optimize n (minimize the number of monkeys needed because I am on a tight budget). So I decided to use a Markov Chain model of order 1 to do so. If you are unfamiliar with Markov Chains here is a very nice explanation of the models here.

The training set of the emission probability matrix, consist on a parsed version of Dracula (chapters 1 to 3, no punctuation signs, lowercase characters only)

The emission probability matrix of the Markov Chain ensures that the transition from one character to another character is constrained by previous character and this relation is weighted based on the frequencies obtained in the training text.

It is like having a keyboard with lights for each key, after "a" is pressed, the light intensity of each key would be proportional of what characters are more likely to appear after an "a". For example "b" would have more light than "a", because it is more common to find words having *a-b* than *a-a*.

# *******************
# RESULTS
# *******************

1) Plot the distribution of characters in the uniform model

Distribution of characters after 10,000 iterations using the uniform model

2) Plot the emission matrices

A) As expected, the transition from one character to another character is constrained by previous character and this relation is weighted based on the frequencies obtained in the training text. B) in the uniform model each character has the same probability to be typed and does not depend on the previous character.

3) Compare the performance of the two models

In this plot I am comparing the number of monkeys (log10(x)) required to type the target words (indicated in red text) using the Markov Chain model and the uniform model. In general the Markov Chain model requires less monkeys in words that are likely to appear in the training set, like "by", "the", "what" , "where" and "Dracula". On the other hand, words that only have one character like "a", given that there's no prior information the models perform equally. Now another interesting example is the word "linux", in which is not very likely to appear in the training set and therefore the models perform likely equally. The extreme case example is the word "omglolbbq", in which the Markov Chain model performs worse than the uniform model due of the very low probability of this word to happen, so it is penalized and I will need more monkeys to get this target word

# *******************
# SOURCE AND FILES
# *******************

Source and files

Benjamin

To cooperate of defect (besides of coding): Prisoners dilemma, a game theory example in R

2014-08-10T13:15:00.001-05:00

Hello Computer Science and/or R enthusiasts.

This week I had the opportunity to try something that was in my To-Do list a while ago. The idea came almost instantly after reading Dr. Richard Dawkins book, The Selfish Gene (which was a BD gift, thanks Andy).

I feel the obligated necessity to program my own implementation of the prisoners dilemma and make my own version of the contest. To try different parameters of greediness and vindictiveness, just coding in the name of science.

You can download the code here

*****************

INTRODUCTION

*****************

Here's a basic definition of prisoners dilemma that can be found in this Wikipedia article:

The structure of the traditional Prisoners’ Dilemma can be generalized from its original prisoner setting. Suppose that the two players are represented by the colors, red and blue, and that each player chooses to either "Cooperate" or "Defect".

If both players cooperate, they both receive the reward, R, for cooperating. If Blue defects while Red cooperates, then Blue receives the temptation, T payoff while Red receives the "sucker's", S, payoff.

Similarly, if Blue cooperates while Red defects, then Blue receives the sucker's payoff S while Red receives the temptation payoff T. If both players defect, they both receive the punishment payoff P.

And to be a prisoner's dilemma game in the strong sense, the following condition must hold for the payoffs:

T > R > P > S

*****************

MATERIALS

*****************

I decided to program three basic strategy functions:

tit.for.tat.bot: this simple strategy repeats opponent's last choice

This bot/strategy is parameter-free

greedy.bot: this strategy is affected by the parameter of greedy.level, which increases the probability of the bot to "defect" when this parameter is close to 1.0. If this parameter is set to 0.5, then the bot will choose randomly.

This bot does not care about the previous action of the other player, it only decides to "defect" or "cooperate" based on its own greediness.

vindictive.bot: this strategy is affected by the parameter of vindictive.level, which is only used when the previous action of the other player is "defect" in the same sense as the greedy.level affects the greedy.bot, otherwise it will "cooperate".

This means that if the other player plays "cooperate", this bot will "cooperate" too, but at any moment when the other player plays "defect", this bot will trigger its own vindictive.level to decide if it will play "defect" or to "cooperate".

All bots have three parameters as arguments:

action: last action of the other player
greedy.level*
vindictive.level

*The greedy.level parameter in the tit.for.tat.bot and the vindictive.bot is only used when this bots are playing as player.1 and have to decide the first action/move against player.2. After that event, this parameter is no longer used for these bots.
Each point in the plots means 1000 (iter=1000) plays of player.1 vs player.2
Parameters for the payoff matrix: T=20, R=5, P=1 and S=0

*****************

METHODS

*****************

Try different pairwise comparisons of strategies varying the greediness and the vindictiveness of the bots
Make our clandestine fight club with the bots (1st RULE: You do not code about FIGHT CLUB)

*****************

RESULTS

*****************

1.1) GREEDY BOT VS TIT FOR TAT BOT

Using different thresholds of greediness in a greedy.bot vs tit.for.tat.bot. we can see that the tit.for.tat.bot handles very well when this bot increases its greediness. this bot cooperates when the adversary cooperates, but reacts when the other player decides to "defects"

1.2) VINDICTIVE BOT VS TIT FOR TAT BOT

Using different thresholds of vindictive.level in a vindictive.bot vs tit.for.tat.bot. A) We can see that the vindictive.bot does not "defect" unless provoked. Therefore as these bots "cooperate" along all the simulations, they get the score for 5000 which is the maximum score for mutual cooperation along 1000 iterations. B) For this experiment, we are increasing the greedy.level and the vindictive.level by steps of 0.05. The greedy.level in the vindictive.bot only determines the probability to "defect" in the first move only and only if the vindictive.bot is playing as player.1. We can see that the tit.for.tat.bot handles very well when this bot increases its greediness and the vindictiveness. This bot cooperates when the adversary cooperates, but reacts when the other player decides to "defects"

1.3) GREEDY BOT VS VINDICTIVE BOT

Using different thresholds of greedy.level in greedy.bot vs the effect of vindictive.level in a vindictive.bot. A) greedy.level and vindictive.level increased proportionally equal among the two bots. we can see that the greedy.bot when becomes more greedy it outperform the vindictive.bot because the vindictive.bot only triggers its self-defense strategy when provoked, and the probability to take revenge is based on the vindictive.level. B) greedy.level and vindictive.level increased proportionally inverse among the two bots. This means that the most greedy bot will face with the most forgiving bot.

2.1) FIGHT CLUB: ROUND 1 (0.7 to 1.0 in greedy and vindictive.level)

Tournament of all vs all using a range of parameters from 0.7 to 1.0 in greedy and vindictive.level A) Scatter-plot of each strategy. in x-axis is the mean score of the strategy versus all other strategies when it is playing as player.1 and y.axis when is playing as player.2. We can see that the greedy.bot.1.0 (greedy.level=1.0) has the lowest score followed by the vindictive.bots and the tit.for.tat.bot. The bots that gets higher scores are the ones that have a greedy personality, but not "so-greedy" after all, for example the most successful bot only has 0.7 of probability of "defect". This experiment is in someway biased because the greedy.bots are taking advantage by exploiting the "nice bots" and getting higher scores for that reason. B) This barplot shows how is the effect of the greediness in the greedy.bot, as it becomes more greedy, the scores tend to be lower. The opposite behavior is shown for the vindictive.bot

2.2) FIGHT CLUB: ROUND 2 (0.8 to 1.0 in greedy and vindictive.level)

To avoid that bias, now we show the experiment with more fierce adversaries

Tournament of all vs all using a range of parameters from 0.8 to 1.0 in greedy and vindictive.level A) Scatter-plot of each strategy. in x-axis is the mean score of the strategy versus all other strategies when it is playing as player.1 and y.axis when is playing as player.2. We can see that the greedy.bot.1.0 (greedy.level=1.0) has the lowest score followed by the vindictive.bots and the not-so-greedy bots. Now the bot that wins the contest is the simplest strategy, tit.for.that. B) This barplot shows how is the effect of the greediness in the greedy.bot, as it becomes more greedy, the scores tend to be lower. The opposite behavior is shown for the vindictive.bot. In general we can conclude that the the vindictive.bots performs better than the greedy.bots because they tend to cooperate unless provoked, so, they are nicer bots.

*****************

DISCUSSION AND CONCLUSIONS

*****************

NICE "BOTS" FINISH FIRST, as Dr Dawkins concluded in his book.

Dr. Dawkins illustrates the four conditions of tit for tat.

Unless provoked, the agent will always cooperate.
If provoked, the agent will retaliate.
The agent is quick to forgive.
The agent must have a good chance of competing against the opponent more than once.

In conclusion, it is better to be a nice human than a greedy one.

Benjamin

UPDATE: 100 Prisoners, 100 lines of code, a game theory example in R

2014-08-01T18:53:00.000-05:00

I was looking for some Game-Theory post in R-bloggers and found this very cool post, obviously I decided to read the post, to look at the code and implement a version of it.

The idea of the post was very cool for me that I thought it would be a nice idea to get more "juice" from it.

All the credit for this post, is for its original author.

The code is available here.

As Matt Asher describes in his original post:

# ***************
# INTRODUCTION
# ***************

In this game, the warden places 100 numbers in 100 boxes, at random with equal probability that any number will be in any box. Each convict is assigned a number. One by one they enter the room with the boxes, and try to find their corresponding number. They can open up to 50 different boxes. Once they either find their number or fail, they move on to a different room and all of the boxes are returned to exactly how they were before the prisoner entered the room.

The prisoners can communicate with each other before the game begins, but as soon as it starts they have no way to signal to each other. The warden is requiring that all 100 prisoners find their numbers, otherwise she will force them to listen to hundreds of hours of non-stop, loud rock musician interviews. Can they avoid this fate?

The first thing you might notice is that if every prisoner opens 50 boxes at random, they will have a 0.5 probability of finding their number.

There's actually a strategy that can guarantee all the prisoners to get their numbers.

For this post, we will call that strategy, model, and the strategy with no strategy at all, we will call it random.model.

The model strategy consist in that at the beginning, each prisoner will open the box based on their number, for example, the prisoner 42 will open the box at the position 42, after open it, say it contains the number 12, and given that 42 is not equal to 12, he will have to open the box at position 12. This will happen until he had opened 50 boxes or found his number.

The random.model strategy consist in no strategy at all, only opening boxes at random until he had opened 50 boxes or found his number.

Every time a prisoner found his number we will call this event a success.

Our interest is to find the strategy that maximizes the success rate. We will run each simulation 1000 times.

# ***************
# RESULTS
# ***************

Let's take a look at the performance of each model individually with this histogram:

We can see that the random.model behaves as we were expecting assuming that by chance we have a 0.5 chance of getting our number by selecting 50 random boxes, On the other hand, the other strategy has a very different distribution, there's a gap between 51 and 99 with 0 of success rate!, any suggestion to explain how is this possible?

Another way to plot the distribution of the success rates obtained from each model is by looking at a density plot:

This histogram is telling the same story than the previous plots, but in the same plot:

And as usual, I always left the best plot for the grand finale.

What we can say about this plot?, before going there, I will continue using Matt's explanation:

The theoretical success rate of the model strategy is about 31%

One way to look at the distribution of numbers in boxes is to see it as a permutation of the numbers from 1 to 100. Each permutation can be partitioned into what are called cycles. A cycle works like this: pick any number in your permutation. Let's say it's 23. Then you look at the number the 23rd place (ie the number in the 23rd box, counting from the left). If that number is 16, you look at the number in the 16th place. If that number is 87, go open box number 87 and follow that number. Eventually, the box you open up will have the number that brings you back to where you started, completing the cycle. Different permutations have different cycles.

The key for the prisoner is that by starting with the box that is the same place from the left as his number, and by following the numbers in the boxes, the prisoner guarantees that if he is in a cycle of length less than 50, he will eventually open the box with his number in it, which would complete the cycle he began. One way to envision cycles of different lengths is to think about the extreme cases. If a particular permutation shifted every single number over one to the left (and wrapped number 1 onto the end), you would have a single cycle of length 100. Box 1 would contain number 2, box 2 number 3 and so on. On the other hand, if a permutation flipped every pair of consecutive numbers, you would have 50 cycles, each of length 2: box 1 would have number 2, box 2 would have number 1. Of course if your permutation doesn't change anything you have 100 cycles of length 1.

As you can see from the histogram of the model strategy, you can never have between 50 and 100 winning prisoners. Anytime you have a single cycle of length greater than 50, for example 55, then all 55 prisoners who start on that cycle will fail to find their number. If no cycles are longer than 50, everyone wins.

Now, looking at the boxplot, OK we would say, using the model strategy we have like a 30% percent of chance that every one of the prisoners will get their number, but on the other hand, the random strategy guarantees that at least every prisoner will have 50% chance of getting his number, and therefore win.

So, the model strategy is like an all of nothing strategy, you can win the lottery or die trying. This is observable looking at the large standard deviation in the boxplot, on the other hand, the random.model strategy is very consistent with low dispersion values. This means that if the prisoners decides to use the model strategy, they will have to face tree possible outcomes, or everybody wins, or nobody does, or only a few does. On the other hand, with the random.model strategy, they will expect that at least 50 prisoners will win, thinking as a prisoner as a logical individual, I think this is the best strategy he can use.

Write your opinions.

Thanks to Matt Asher for the idea

Benjamin

What are the names of the school principals in Mexico?, If your name is Maria, probably this post will interest you. Trends and cool plots from the national education census of Mexico in 2013

2014-07-07T12:28:00.000-05:00

I will start this post with a disclaimer:

The main intention of the post is to show how is the distribution of the school principal names in Mexico, for example, to show basic trends regarding about what is the most common nation-wide first name and so on, also to show trends delimited by state and regions.

These trends in data would answer questions such:

1. Are the most common first names distributed equally among the states?
2. Does the states sharing the same region also share the same "naming" behavior?

Additionally, this post includes cool wordclouds.

Finally and the last part of my disclaimer is that, I am really concerned about the privacy of the persons involved. I am not by any sense promoting the exploitation of this personal data, if you decide to download the dataset, I would really ask you to study it and to generate information that is beneficial, do not join the Dark side.

Benjamin

##################
# GETTING THE DATASET AND CODE
##################

The database is located here
The R code can be downloaded here
Additional data can be downloaded here

All the results were computed exploring 202,118 schools across the 32 states of Mexico from the 2013 census

##################
# EXPLORING THE DATA
# WITH WORDCLOUDS
##################

Here is the wordcloud of names (by name, I am referring to first name only), it can be concluded that MARIA is by far the most common first name of a school principal in all Mexican schools, followed by JOSE and then by JUAN

The following wordcloud includes every word in the responsible_name column (this includes, first name, last names). Now the plot shows that besides the common first name of MARIA, also the last names of HERNANDEZ, MARTINEZ and GARCIA are very common.

##################
# EXPLORING THE FREQUENCY
# OF FIRST NAMES (TOP 30 | NATION-WIDE)
##################

Looking at this barplot, the name MARIA is by far the most common name of the Mexican school's principals, with a frequency ~25,000. The next most popular name is JOSE with a frequency of ~7,500

Looking at the same data, just adjusted to represent the % of each name inside the pool of first names we have that MARIA occupy ~11% of the names pool.

##################
# HEATMAPS OF THE DATA
##################

With this heatmap, my intention is to show the distribution of the top 20 most common first names across all the Mexican states

It can be concluded that there is a small cluster of states which keep the most number of principals named MARIA(but no so fast!, some states, for example Mexico and Distrito Federal are very populated, so I will reduce this effect in the following plot). In summary the message of this plot is the distribution of frequency of the top 20 most frequent first-names across the country.

##################
# CLUSTERS OF THE DATA
##################

For me, a young data-science-padawan, this is my favorite analysis: "hunting down the trends".

The setup of the experiment is very simple, map the top 1,000 most frequent nation-wide names across each state to create a 32 x 1000 matrix (32 states and 1,000 most nation-wide frequent names).

With this matrix, normalize the values by diving each row by the sum of it (this will minimize the effect of the populated states vs the non populated while maintaining the proportion of the name frequencies per state). Then I just computed a distance matrix and plotted it as a heatmap.

What I can conclude with this plot is that, there are clusters of states that seems to maintain a geographical preference to be clustered within its region, this would be concluded that it is likely that states sharing the same regions would be more likely to share the "naming" trends due to some cultural factors (like the cluster that includes Chihuahua, Sonora and Sinaloa). But this effect is not present in all the clusters.

All images can be downloaded in PDF format here, just don't do evil with them!

Plot 1 here
Plot 2 here
Plot 3 here
Plot 4 here
Plot 5 here
Plot 6 here

Benjamin

[SOLVED] Problems compiling Rcpp dependent R packages in Crunchbang Linux 11 (Debian 7.5 (wheezy) 64-bit)

2014-06-18T12:26:00.000-05:00

I had a very weird issue when I tried to compile (install) certain R packages like "wordcloud", "RSNNS" or "GOSemSim" (a Bioconductor package). The installation always ended with a compiling error at the end.

At first, I tried to solve the problem by finding anyone that had faced the same issue, so I Duck-Duck-go-it, Google-it and at least for me, I did not find anyone.

At the end, and which is the core of my post, I solved the problem by removing or commenting all my customizations in the /usr/lib/R/etc/Rprofile.site file (for example, loading libraries automatically).

In summary, using the default Rprofile.site worked for me.

I wrote this post to make a record of the issue, so maybe it would help someone to solve the same in the future.

Benjamin

Convert Ensembl, Unigene, Uniprot and RefSeq IDs to Symbol IDs in R using Bioconductor

2014-02-25T10:35:00.000-06:00

Hello, I have programmed a function that converts different sources of IDs to Symbol IDs.

The input ID types allowed are (at the moment): Ensembl, Unigene, Uniprot and RefSeq.

The code is available clicking here

NOTE: The function depends on the Bioconductor package "org.Hs.eg.db" available here

For example, lets show 10 Ensembl IDs:

> id[1:10]
[1] "ENSG00000121410" "ENSG00000175899" "ENSG00000256069" "ENSG00000171428"
[5] "ENSG00000156006" "ENSG00000196136" "ENSG00000114771" "ENSG00000127837"
[9] "ENSG00000129673" "ENSG00000090861"

And their Symbol IDs:

> res[1:10]

[1] "A1BG" "A2M" "A2MP1" "NAT1" "NAT2" "SERPINA3"

[7] "AADAC" "AAMP" "AANAT" "AARS"

This is a running example of the function to convert Unigene IDs to Symbol IDs (For all the other IDs types, just replace "unigene" to "ensembl" or "refseq" or "uniprot"):

# USAGE EXAMPlE: UNIGENE

require(org.Hs.eg.db)

unigene <- toTable(org.Hs.egUNIGENE)

# extract 100 random unigene entries

id <- unigene[sample(1:length(unigene[,2]),100),2]

id.type <- "unigene"

res <- get.symbolIDs(id,id.type)

Benjamin

Upgrade and update R 2.15 to R 3.0 in Debian Wheezy

2014-02-10T16:51:00.000-06:00

Following the instructions from CRAN, you need to add the R backports in your source list.

FIRST PART: ADD R BACKPORTS:

First, open a Terminal and open the sources.list file:


$ gksudo gedit /etc/apt/sources.list

Then, add these lines at the bottom of the file (Note, I use the Revolution Analytics Dallas, TX server, but this can be easily changed taking a look here for the mirrors):

## R BACKPORTS FOR WHEEZY
deb http://cran.revolutionanalytics.com/bin/linux/debian wheezy-cran3/
#deb-src http://cran.revolutionanalytics.com/bin/linux/debian wheezy-cran3/

SECOND PART: RENAME THE R PACKAGES FOLDER:

There's a folder where R uses to store the packages we download, just rename it to the current version of R. For example, mine was "2.15" and then I just renamed it to "3.0" and was inside this path:

Before:

/home/benjamin/R/x86_64-pc-linux-gnu-library/2.15

After:

/home/benjamin/R/x86_64-pc-linux-gnu-library/3.0

Remember that some packages also needs to install some files in folders that belongs to the root, so, I would recommend to open R in sudo mode (only if you're sure about what you're doing :P) just by executing R this way: "sudo R" and then, in the R console type :

update.packages(checkBuilt=TRUE, ask=FALSE)

THIRD PART: SECURE APT:

The Debian backports archives on CRAN are signed with the key ID 381BA480, to add them, in a Terminal prompt type:

gpg --keyserver pgp.mit.edu --recv-key 381BA480
gpg -a --export 381BA480 > jranke_cran.asc
sudo  apt-key add jranke_cran.asc

FOURTH PART: UPDATE AND UPGRADE R:

Save the file and you can either enter to Synaptic, update the packages list and then just upgrade the packages or in a terminal type:


sudo apt-get update

sudo apt-get upgrade

And that's all.

Benjamin

Optimizing a multivariable function parameters using a random method, genetic algorithm and simulated annealing in R

2013-10-19T11:17:00.000-05:00

Say that you are implementing a non-linear regression analysis, which is shortly described by wikipedia as:

"In statistics, nonlinear regression is a form of regression analysis in which observational data are modeled by a function which is a nonlinear combination of the model parameters and depends on one or more independent variables."

For the training set, we have the following:

And the function to optimize the parameters is:

Which leads us to the following equality:

In other words, we want to optimize the value of theta in order to minimize the sum of the error among y and predicted.y:

Given theta (each parameter a0,..a3 has a range from 0 to 15):

And the error function:

finally, the goal function:

In other words, the goal function searches for the value of theta that minimizes the error.

COMPUTATIONS BEGIN

This is the scatter plot of the training set:

Here is the implementation in R, you can download the file clicking here

Here is a result plot using the genetic algorithm:

Benjamin

Lieutenant Dan You Got a New Interface..

2013-08-17T01:38:00.004-05:00

After some days of thinking, I realize that this blog deserved a little bit more attention, so I decided to change the interface and I'am happy about how it looks now.

Hope to see you again in my personal blog.

keep on programming!

Benjamin

Accuracy versus F score: Machine Learning for the RNA Polymerases

2013-08-16T22:19:00.002-05:00

Hello, today I'm going to show you the difference of using two different common performance measures (useful not only for Machine Learning purposes, is useful in every scientific field). Until now, I have found more the accuracy values than F scores in the performance measuring of some methods which ranges from metaheuristics (Genetic Algorithms fitness functions) to promoter recognition programs, diagnose methods and so on.

But, I would really recommend to avoid using the accuracy measure. The reason is shown below with a nice example in R programming language (all the functions used in the simulation are included, you can download them clicking here).

Case study 1:

Imagine that you are in a Computer Vision project and your task is to "teach" a program to recognize among electric guitars and acoustic guitars showing the program pictures of different guitars.

Suppose that you've already developed that program and now you want to measure the performance of this Boolean classifier (this is for example, you show the program a picture of a an electric guitar, and the program has to decide whether it will recognize and "classify" it as an electric or as an acoustic guitar).

For the function of this post, lets write down some useful concepts

Consider the following:

TP: a true positive is when the program classifies an electric guitar as an electric guitar, we will use the letter "E" to denote the electric guitar "class"

FP: a false positive is when the program classifies an acoustic guitar as an electric guitar, we will use the letter "A" to denote the acoustic guitar "class"

FN: a false negative is when the program classifies an electric guitar as an acoustic guitar

TN: a true negative is when the program classifies an acoustic guitar as an acoustic guitar

Now that we are ready, we shall begin with the calculations

In R, I have simulated the results of the program. Say, for 1,000 electric guitar pictures and 1,000 acoustic guitar pictures

The program prompt the following results:

PREDICTED.E PREDICTED.A
TRUE.E 485 515
TRUE.A 9 991

If you notice, from the 1000 electric guitar pictures, only 485 were labeled as electric (TP=485), the rest were labeled as acoustic (FN=515). I feel bad for the hypothetical programmer of this hypothetical example.

On the other hand, from the 1000 acoustic guitars, 991 were labeled as acoustic (TN=991) and only 9 of them were labeled as electric (FP=9). Well not bad!..... or it is?

The accuracy value of this program is = 0.738

And, for computing the F score is necessary to compute the precision and the recall first, where:

precision = 0.9817814 and recall = 0.485

Then, the F score is equal to 0.6492637

Well, the F scores seems to be more "strict", and in fact it is in comparison of the accuracy performance measure. But this example is not very "cool". Lets pass to the case study 2

Case study 2:

Now we have 1,000 electric guitar pictures and 100,000 acoustic guitar pictures, the confusion matrix of the results are:

PREDICTED.E PREDICTED.A
TRUE.E 493 507
TRUE.A 1017 98983

If you notice, from the 1,000 electric guitar pictures, only 493 were labeled as electric (TP=493), the rest were labeled as acoustic (FN=507)

On the other hand, from the 100,000 acoustic guitars, 98983 were labeled as acoustic (TN=98983) and only 1017 of them were labeled as electric (FP=1017)

Now (cha cha chan!), the performance values are:

Accuracy: 0.9849109
Precision: 0.3264901
Recall: 0.493
F score: 0.3928287

Now you see it?, how come or how is possible that missing almost the 50% of the labels of the electric guitars, the performance of the program in accuracy is almost 0.99?, despite of having a precision and recall not greater than 0.50. Then we have a winner and is the F score measure.

for references visit the following pages:

http://en.wikipedia.org/wiki/Accuracy
http://en.wikipedia.org/wiki/F1_score
http://en.wikipedia.org/wiki/Precision_and_recall

Make your LaTex presentations using the Beamer class

2013-05-25T20:14:00.000-05:00

Hello, I have been working with this tutorial. All you need to do is download the source files and take a look at the source code and the PDF file.

Download link here

Hope you find LaTex as useful as I do

Benjamin

Póster presentado en la XIV Escuela de Otoño en Biología Matemática, México

2012-12-02T12:33:00.001-06:00

Póster presentado en la XIV Escuela de Otoño en Biología Matemática - 8vo Encuentro de Biología Matemática celebrado en San Luis Potosí, S.L.P, México:

"Predicción de promotores RNA POL-II en Drosophila melanogaster utilizando propiedades de señal, contexto y estructura a partir de secuencias nucleotídicas"

Póster disponible aquí

Brownian motion simulation in R

2012-12-02T12:21:00.001-06:00

I have based this post on a very useful piece of code which basically is the core of my own implementation of a Brownian Motion simulation in R. The original reference code <- http://landshape.org/enm/r-code-for-brownian-motion/

According to Wikipedia the mathematical model for Brownian motion (also known as random walks) can also be used to describe many phenomena as well as the random movements of minute particles,
such as stock market fluctuations and the evolution of physical characteristics in the fossil record. The simple form of the mathematical model for Brownian motion has the form:

S_t = eS_t-1

where e is drawn from a probability distribution.

The source code is here

After loading the source code, there are two functions:

The first one, brownian will plot in an R graphics window the resulting simulation in an animated way.

The second function, export.brownian will export each step of the simulation in independent PNG files.

Example of running:

> source("brownian.motion.R")
> brownian(500)

The second function will produce this output
> export.brownian(500)
CODE:

# *******************************
# BROWNIAN MOTION SIMULATION
# December 2012 | Benjamin Tovar
# *******************************
#
#   REFERENCES
#   http://landshape.org/enm/r-code-for-brownian-motion/
#
#   According to Wikipedia the mathematical model for Brownian motion 
#   (also known as random walks) can also be used to describe many 
#   phenomena as well as the random movements of minute particles, 
#   such as stock market fluctuations and the evolution of physical 
#   characteristics in the fossil record. The simple form of the 
#   mathematical model for Brownian motion has the form:
#
#    S_t = eS_t-1
#
#    where e is drawn from a probability distribution.
#
#######################################################################

brownian <- function(n.times){
    x <- y <- x.new <- y.new <- x.new.p <- y.new.p <- vector()
    for(i in 1:n.times){
        # Initialize variables
        x <- rnorm(1)
        y <- rnorm(1)
        # concatenate variables 
        # to increase the vector size
        x.new <- c(x.new,x)
        y.new <- c(y.new,y)
        # sum the vector numbers
        x.new.p <- cumsum(x.new)
        y.new.p <- cumsum(y.new)  
        # plot the model
        plot(x.new.p,y.new.p,type="b",
             main=paste("Brownian motion simulation in R\nTime =",i,sep=" "),
             xlab="x coordinates",ylab="y coordinates",
             col=c(rep("gray",i-1),"red"),
             pch=c(rep(20,i-1),1))    
    }
}

# Test the function
# brownian(500)

# ****************************************
# EXPORT BROWNIAN MOTION SIMULATION IMAGES
# ****************************************

export.brownian <- function(n.times){
    x <- y <- x.new <- y.new <- x.new.p <- y.new.p <- vector()
    for(i in 1:n.times){
        # Initialize variables
        x <- rnorm(1)
        y <- rnorm(1)
        # concatenate variables to increase the
        # vector size
        x.new <- c(x.new,x)
        y.new <- c(y.new,y)
        # sum the vector numbers
        x.new.p <- cumsum(x.new)
        y.new.p <- cumsum(y.new)  
        # plot the model
        png(paste("image",i,"png",sep="."),width=600,height=600)
            plot(x.new.p,y.new.p,type="b",
                 main=paste("Brownian motion simulation in R\nTime =",
                 i,sep=" "),
                 xlab="x coordinates",ylab="y coordinates",
                 col=c(rep("gray",i-1),"red"),
                 pch=c(rep(20,i-1),1))
                 cat("image",i,"DONE",date(),"\n")
        dev.off()
    }
}

# Test the function
# export.brownian(500)

Extracting upstream regions of a RefSeq human gene list in R using Bioconductor

2012-07-29T13:22:00.002-05:00

Suppose that you want to do local mapping of upstream regions of a given RefSeq IDs in a particular genome in R using Bioconductor. Download the script here.

In this case, you may take a look at the Bioconductor AnnotationData Packages here: http://www.bioconductor.org/packages/release/data/annotation/

The goal of this post is that for example I have the following RefSeq IDs and want to extract 250 bases upstream of each gene in a single list with another useful information such the entrez.id, the symbol and the gene description.


# RefSeqs IDs:

gene.list.refseq <- c("NM_003588","NM_001145436", "NM_001135188","NM_020760","NM_173362", "NM_198393","NM_022736","NM_025074","NM_033449","NM_015726", "NM_022110","NM_016478","NM_020634","NM_002291","NM_000418", "NM_001862","NM_017752","NM_006591","NM_000124","NM_144610")



# How many bases upstream of each gene:    

bases.upstream <- 250

Before starting, please download the following packages in R:


       source("http://bioconductor.org/biocLite.R")

       biocLite("BSgenome.Hsapiens.UCSC.hg19")                               

       biocLite("Biostrings")

       biocLite("org.Hs.eg.db")

Load the function using:


source("extract.five.utr.sequence.R")

And finally make the computations:


output.sequences <- extract.five.utr.sequence(gene.list.refseq,bases.upstream)

CODE:

# ******************************************************************************
#       FUNCTION: extract.five.utr.sequence | 29/JULY/2012 | BENJAMIN TOVAR
# ******************************************************************************
extract.five.utr.sequence <- function(gene.list.refseq,
                                      number.bases.upstream=1000){
    
    
    # *****************************************************************
    # RUNNING NOTES: Please download this packages from Bioconductor
    # http://www.bioconductor.org/packages/release/data/annotation/
    # ***************************************************************** 
    #   source("http://bioconductor.org/biocLite.R")
    #   biocLite("BSgenome.Hsapiens.UCSC.hg19")                                
    #   biocLite("Biostrings")
    #   biocLite("org.Hs.eg.db")
                                      
    # *****************************************************************
    # RUNNING EXAMPLE
    # *****************************************************************                              
    #
    ## Extract 250 bases upstream of each gene in gene.list.refseq 
    #
    # gene.list.refseq <- c("NM_003588","NM_001145436", "NM_001135188","NM_020760","NM_173362",   
    #                    "NM_198393","NM_022736","NM_025074","NM_033449","NM_015726",   
    #                    "NM_022110","NM_016478","NM_020634","NM_002291","NM_000418",   
    #                    "NM_001862","NM_017752","NM_006591","NM_000124","NM_144610") 
    #
    # bases.upstream <- 250
    #  
    # result <- extract.five.utr.sequence(gene.list.refseq,bases.upstream)                             
                             
    # *****************************************************************
    # LOAD THE LIBRARIES
    # ***************************************************************** 
    cat("Loading libraries",date(),"\n")   
    # human genome DNA sequences
        library(BSgenome.Hsapiens.UCSC.hg19) 
    # human genome wide annotations 
        library(org.Hs.eg.db) 
    # load IDs, symbol and gene descriptions
        refseq.id <- toTable(org.Hs.egREFSEQ)
        symbol <- toTable(org.Hs.egSYMBOL)
        gene.description <- toTable(org.Hs.egGENENAME)
    
    # *****************************************************************
    # LOAD THE UPSTREAM SEQUENCES
    # *****************************************************************

    if(number.bases.upstream <= 1000){
        # load the 1000 bases upstream of all genes in the human genome 
        # in the BSgenome.Hsapiens.UCSC.hg19 package
        upstream <- Hsapiens$upstream1000
    }
    
    if(number.bases.upstream > 1000 && number.bases.upstream <= 2000){
        # load the 1000 bases upstream of all genes in the human genome 
        # in the BSgenome.Hsapiens.UCSC.hg19 package
        upstream <- Hsapiens$upstream2000
    }
    
    if(number.bases.upstream > 2000 && number.bases.upstream <= 5000){
        # load the 1000 bases upstream of all genes in the human genome 
        # in the BSgenome.Hsapiens.UCSC.hg19 package
        upstream <- Hsapiens$upstream5000
    }        
    
    if(number.bases.upstream > 5000){
        cat("ERROR: The number of bases upstream of 5'UTR is too large (max value = 5000)\n")
        return(0);
    } 
    
    if(number.bases.upstream < 1){
        cat("ERROR: The number of bases upstream of 5'UTR cannot be less than 1 nucleotide (min value = 1)\n")
        return(0);
    }                   
        
    # extract the names of each gene
    names.genes.upstream <- names(upstream)

    # extract the refseq of each gene  
    refseq.upstream.id <- vector("character",length(names.genes.upstream))
    for(i in 1:length(names.genes.upstream)){
        refseq.upstream.id[i] <- gsub(", ","_",toString(unlist(strsplit(names.genes.upstream[i],
                                      "\\_"))[c(1,2)],sep=""))
    }    
    names(refseq.upstream.id) <- 1:length(refseq.upstream.id)

    # *****************************************************************
    # CREATE OUTPUT OBJECT
    # *****************************************************************    
    output.params <- c("entrez.id","symbol","gene.description","five.UTR.sequence")
                        
    output.list <- list()
    for(i in 1:length(gene.list.refseq)){
        output.list[[i]] <- list()
        for(j in 1:length(output.params)){
            output.list[[i]][[j]] <- list()
        }
        names(output.list[[i]]) <- output.params
    }
    names(output.list) <- gene.list.refseq

    # *****************************************************************
    # LET THE COMPUTER DECIDE
    # *****************************************************************     
    cat("Extracting data",date(),"\n")
    
    for(i in 1:length(gene.list.refseq)){
        cat(rep(".",1))

        # extract the entrez.id of each gene in geneList 
        output.list[[i]]$entrez.id <- refseq.id[which(refseq.id$accession==gene.list.refseq[i]),]$gene_id

        # extract the symbol of each gene in geneList 
        output.list[[i]]$symbol <- symbol[which(symbol==output.list[[i]]$entrez.id),]$symbol
        
        # extract the gene description of each gene in geneList 
        output.list[[i]]$gene.description <- gene.description[which(gene.description==output.list[[i]]$entrez.id),]$gene_name  
        
        # Return the index of the target genes
        index.target.gene <- sequence <- NA;
        index.target.gene <- as.numeric(names(refseq.upstream.id[which(refseq.upstream.id==gene.list.refseq[i])]))
        
        # NOTE: some genes have repeated RefSeq.id which corresponds to the same gene in other Chromosomes.
        # by using index.target.gene[1] we are selecting only the first RefSeq entry.
        sequence <- upstream[index.target.gene[1]]
        sequence.length <- nchar(sequence)
        start.position <- ((sequence.length-number.bases.upstream)+1)
        sequence <- substring(toString(sequence),1:sequence.length,1:sequence.length)
        sequence <- sequence[start.position:sequence.length]
        output.list[[i]]$five.UTR.sequence <- sequence
    }
    cat("\n")
    cat("Computations DONE",date(),"\n")
    return(output.list)
}

Benjamin.

Handling large FASTA sequence datasets in R: Shuffle and retrieve "n" number of sequences of fixed length from the whole FASTA file and export them in a new FASTA file.

2012-06-23T11:17:00.001-05:00

When you are working with large FASTA datasets is probable to find out that the sequences are in sort of a mixed quality (obviously, depending on your scientific question),

I mean for example, imagine that you retrieve the whole collection of exons of a given organism and suppose that the FASTA file is 50mb and there are included ~50,000 DNA sequences but, if you look at them, you may find that there are sequences much larger than others and others will be probably 50 bases long.

Well that's what I mean with mixed quality and therefore, a sort of filter might be very helpful because you may find more informative to put a threshold and say "hey, from the pool of ~50K sequences I only want 5K sequences randomly chosen (from the entire FASTA dataset) and every sequence must have a fixed length, say 1000 bases long".

That's why I wrote some code in R language, it's a function titled "shuffleAndExtract" and you can download the example set and the source here.

Here is the description of the function:

handleFastaDatasets.R

shuffleAndExtract: This function in R is designed to open a fasta file dataset, shuffle the sequences and extract the desired sequences wanted by the user to generate a new dataset of fixed size (number of required sequences) and with the same length for each sequence.

And, after you download the example set and the R source file, you can try to run a very simple example:

NOTE: my implementation depends on the function "seqinr", if is not installed, you may do this before all the magic begin with a simple install.packages("seqinr").

# run example:
source("handleFastaDatasets.R")
shuffleAndExtract("example.fasta",1000,200)

The arguments of the function are:

inputFastaFile: name of the input fasta file
numberOfoutputSeqs: number of desired sequences in the output file
lengthOutputSeqs: fixed length of every sequence in the output file
initialPos: Position where the new window sizing will begin, default = 1
outputFileName: name of the output file, by default will be (e.g): "inputFastaFile.fasta.output.fasta"

CODE:

################################################################################

shuffleAndExtract <- function(inputFastaFile,
                              numberOfoutputSeqs,
                              lengthOutputSeqs,
                              initialPos=1,
                              outputFileName = paste(inputFastaFile,"output",
                                                     "fasta",sep=".")
                              ){

    #
    #    handleFastaDatasets.R   
    #
    #    This function in R is designed to open a fasta file dataset, shuffle the 
    #    sequences and extract the desired sequences wanted by the user to generate 
    #    a new dataset of fixed size (number of required sequences) and with the 
    #    same length for each sequence. 
    #
    #    Author: Benjamin Tovar
    #    Date: 22/JUNE/2012
    #
    ################################################################################

    #    run example:
    #    source("handleFastaDatasets.R")
    #    shuffleAndExtract("example.fasta",1000,200)

    # inputFastaFile: name of the input fasta file
    # numberOfoutputSeqs: number of desired sequences in the output file 
    # lengthOutputSeqs: fixed length of every sequence in the output file
    # initialPos: Position where the new window sizing will begin, default = 1
    # outputFileName: name of the output file, by default will be (e.g):
    #                 "inputFastaFile.fasta.output.fasta"

    cat("*** Starting computations |",date(),"****\n")    
        
    # Load the seqinr package
    require(seqinr)

    # Load the large seq and not shuffled dataset
    inputSeqs <- read.fasta(inputFastaFile)

    cat("\tProcessing for",length(inputSeqs),"sequences |",date(),"\n")    

    # Extract the length of every sequence and store the results in a vector.
    inputSeqsSize <- rep(NA,length(inputSeqs))
    for(i in 1:length(inputSeqs)){ 
        inputSeqsSize[i] <- summary(inputSeqs[[i]])$length
    }

    # Extract the index of the sequences which are longer than threshold
    inputSeqsLongSeqIndex <- which(inputSeqsSize > (lengthOutputSeqs+initialPos))

    # randomly pick numberOfoutputSeqs indexes that will be used to create 
    # the output dataset 
    inputSeqsIndex <- sample(inputSeqsLongSeqIndex,numberOfoutputSeqs,rep=F)
    
    # Store the Fasta header of each selected sequence in a vector
    inputSeqsIndexNames <- rep(NA,numberOfoutputSeqs)
    
    # create output object 
    outputSeqs <- list()
    for(i in 1:numberOfoutputSeqs){
        # Extract the fasta headers 
        inputSeqsIndexNames[i] <- attr(inputSeqs[[inputSeqsIndex[i]]],"name")
        # Extract the sequence
        outputSeqs[[i]] <- inputSeqs[[inputSeqsIndex[i]]][initialPos:((initialPos+lengthOutputSeqs)-1)]
    }
 
    # Export the sequences in a new fasta file 
    write.fasta(outputSeqs,inputSeqsIndexNames,outputFileName)
    
    cat("\n*** DONE |",date(),"****\n")
}


# 22 Jun 2012 | Benjamin Tovar

Benjamin

First anniversary of my blog

2012-05-06T00:04:00.000-05:00

Exactly today (5th/May), a year ago I decided to start with a project, mostly motivated in a way that I feel that probably some of my post will help someone in the future (bots included).

My motivation comes from, because my knowledge about Biology is ~80% and the other ~20% is about Computer Science, I thought that it will be a very cool idea to share some code, ideas, tutorials and even random messages for people that have not been surrounded with bunches of code in front of a UNIX Terminal.

I promise to still publishing posts, and in the same way I still continue my own training, the awesomeness :P of the posts will increase too.

Wishing you the best for every person which has visited my blog. Thank you so much and keep coding!

Benjamin.

My poster "Mining Biosequences" (Spanish) at the Scientific Poster Contest due to the BSc in Genomic Biotechnology anniversary at Faculty of Biological Sciences in UANL

2012-05-02T20:05:00.000-05:00

My poster is an introduction to HMM using some awesome material from Eddy S (among other awesome authors)., hope you like it.

Click the image to download the poster.

Benjamin

Calculate the average distance between a given DNA motif within DNA sequences in R

2012-04-21T15:40:00.000-05:00

Suppose that we want to calculate the expected distance of a DNA motif within a DNA target sequence, if we know the composition bias or the probability distribution (multinomial model) we can compute it just fine.

Download the R code <- here

FIRST PART

For example, supose that we want to compute the expected distance of the motif "GATC" in a sequence composed of 10,000 bases given that the whole sequence follows a probability distribution of p(A) = 0.3, p(C) = 0.2, p(G) = 0.2, p(T) = 0.3.

So open an R prompt and load the functions:


source("motifOccurrence.R")

Then, enter the initial values:


lengthSeq <- 10000

motif <- c("G","A","T","C")

probDistr <- c(0.3,0.2,0.2,0.3)

And finally, compute the average expected distance of every occurrence of the motif inside the target sequence. Using the following formula (this computation corresponds to the "computeExpectedDistance" function of the R script "motifOccurrence.R"):

expectedDistance = lengthDNAseq/(lengthDNAseq*p))

where "p" stands for the joint probability of the motif, in other words: p = p(GATC) = p(G)*p(A)*p(T)*p(C) = 0.2*0.3*0.3*0.2 = 0.0036

To easily compute the expected distance in R, type:


expectedDist <- computeExpectedDistance(probDistr,lengthSeq,motif)

As you see, it returns the value of "277", this number means that, for the "GATC" motif inside a sequence of 10,000 bases with a composition bias of p(A) = 0.3, p(C) = 0.2, p(G) = 0.2, p(T) = 0.3 we may expect a distance of 277 bases between each "GATC".

Or, a little more graphic:

277 bases | GATC | 277 bases | GATC | .....

SECOND PART

Another way to compute this (even though it involves more computations) we can simulate X number of DNA sequences with a fixed length with an equal probability distribution per sequence and extract the coordinates of the motif within each sequence to finally compute the average distance of the motif.

There is a function titled "iterateComputeDistance" to do the calculations for you. Add the next parameter to the R environment:


iterations <- 100

Compute the average distance of the "GATC" motif within 100 DNA sequences (the other parameters remain equal)


expectedDistWithSimSeqs <- iterateComputeDistance(probDistr,lengthSeq,motif,iterations)

As we expected, the results of the two approaches are highly similar (ouuu yeah!)

THIRD PART

But, what happens when we already have a sequence and want to know the expected distance of that motif inside of it?.

Just like "Hey dude, I have an E.coli plasmid DNA sequence and want to know the average distance of the GATC motif".

Lets test using the sequence "Escherichia coli 2078 plasmid pQNR2078 complete sequence" <- http://www.ncbi.nlm.nih.gov/nuccore/HE613857.1

Ok, use the "ape" library to import the sequence to the R environment (if this library is not installed, type: install.packages("ape"))

NOTE: the GenBank sequences are in lowecase, so it will be needed to use a motif in lowercase to do the right computations.

Import the library:


require("ape")

Import the sequence:


plasmid <- read.GenBank("HE613857.1")

plasmidDNA <- as.character.DNAbin(plasmid)

plasmidDNA <- plasmidDNA[[1]]

motifEcoli <- c("g","a","t","c")

Get the coordinates:
plasmidDNAcoord <- coordMotif(plasmidDNA,motifEcoli)
Get the average distance between the motif occurrences.


plasmidDNAmotifDistance <- computeDistance(plasmidDNAcoord)

> plasmidDNAmotifDistance
[1] 270

Compute the number of occurrences of the motif among the plasmid (the result is 151 occurrences):

The number of occurrences of the motif in R is:


(length(plasmidDNAcoord)-1)

So, the main distance between the motif "gatc" finally is 270 bases and we are done :D

CODE:

#
#    Script: motifOccurrence.R
#    Author: Benjamin Tovar
#    Date: 21/April/2012
#
################################################################################

#                       ############
#                        FUNCTIONS:
#                       ############

##############################################################################
iterateComputeDistance <- function(multinomialDNAmodel,
                                   lengthDNAseq,
                                   motif,
                                   numberOfIterations){

    # This function returns the mean distance 
    # of a given motif given X number of DNA sequences given a multinomial model
    # (probability distribution of each base).
    
    # So, it will generate X number of DNA sequences using a given 
    # probability distribution and then it will compute the distance among 
    # that mofit within the total set of sequences to finally returns 
    # the average distance of the motif.
    
    result <- rep(NA,numberOfIterations)
    for(i in 1:numberOfIterations){
        currentGenome <- NA
        currentCoordinatesOfMotif <- NA        
        currentSequence <- sample(c("A","C","G","T"),
                                    lengthDNAseq,rep=T,
                                    prob=multinomialDNAmodel)
        currentCoordinatesOfMotif <- coordMotif(currentSequence,motif)
        result[i] <- computeDistance(currentCoordinatesOfMotif)
        cat(" *** Iteration number: ",i," completed *** | average distance = "
            ,result[i],"\n")  
    }
    result <- trunc(mean(result))
    cat(" \n*** Computation status: DONE ***\n\n")
    return(result)  
}
##############################################################################
coordMotif <- function(targetSequence,motif){

    # This function returns the coordinates of the motif of study in a target 
    # DNA sequence. In other words, if I found the motif, tell me exactly in
    # which position of the DNA sequence is.
    
    lengthMotif <- length(motif)
    lengthTargetSeq <- (length(targetSequence)-lengthMotif)
    motif <- toString(motif)
    motif <- gsub(", ","",motif)
    res <- 1    
    for(i in 1:lengthTargetSeq){
        currentTargetSeq <- targetSequence[i:(i+(lengthMotif)-1)]
        currentTargetSeq <- toString(currentTargetSeq)
        currentTargetSeq <- gsub(", ","",currentTargetSeq)
        if(currentTargetSeq == motif){
            res[(length(res)+1)] <- i
        }
    }
    return(res)
}
##############################################################################
computeDistance <- function(coordinatesOfMotif){
    
    # This function returns the mean distance 
    # of a given motif given its coordinates within a target DNA sequence.
    # In other words, If I already got a list with the coordinates where the 
    # motif is inside a DNA sequence, tell me the average distance between
    # this coordinates to get the expected distance of that motif.

    currentDistance <- rep(NA,(length(coordinatesOfMotif)-1))
    lengthCoord <- length(currentDistance)
    for(i in 1:lengthCoord){
        currentDistance[i] <- coordinatesOfMotif[i+1]-coordinatesOfMotif[i]
    }
    res <- trunc(mean(currentDistance))
    return(res)   
}
##############################################################################
computeExpectedDistance <- function(multinomialModel,
                                    lengthDNAseq,
                                    motif){
                                    
    # This function computes the expected distance of a given motif in a DNA
    # sequence given its multinomial model (probability distribution of 
    # each base)
    
    # Convert the motif into an index                                       
    motifIndex <- gsub("A",1,motif); motifIndex <- gsub("C",2,motifIndex)
    motifIndex <- gsub("G",3,motifIndex); motifIndex <- gsub("T",4,motifIndex)
    motifIndex <- as.numeric(motifIndex)
    # Compute p value of the motif given the multinomial model
    p <- rep(NA,length(motif))
    for(i in 1:length(motifIndex)){
        p[i] <- multinomialModel[motifIndex[i]]
    }    
    p <- prod(p)
    result <- trunc(lengthDNAseq/(lengthDNAseq*p))
    return(result)
}                                   
##############################################################################

# Benjamin

Benjamin

Introduction to Markov Chains and modeling DNA sequences in R

2012-04-13T22:57:00.001-05:00

Markov chains are probabilistic models which can be used for the modeling of sequences given a probability distribution and then, they are also very useful for the characterization of certain parts of a DNA or protein string given for example, a bias towards the AT or GC content.

For now, I am going to introduce how to build our own Markov Chain of zero order and first order in R programming language. The definition of a zero order Markov Chain relies in that, the current state (or nucleotide) does not depends on the previous state, so there's no "memory" and every state is untied.

For the first order Markov Chain the case is different because the current state actually depends only on the previous state. Given that points clear, a second order Markov Model will be a model that reflects that the current state only depends on the previous two states before it (This model will be useful for the study of codons, given that they are substrings of 3 nucleotides long).

Download the scripts (R script and graphviz scripts) here

Example of a Markov Chain of zero order (the current nucleotide is totally independent of the previous nucleotide).

The multinomial model is:
p(A)+p(C)+p(G)+p(T) = 1.0
0.4 +0.1 +0.1 +0.4 = 1.0

Example of the structure of the zero order model:

Example of a Markov Chain of first order (the current nucleotide only depends on the previous nucleotide).

The multinomial model per base is:
p(A|A)+p(C|A)+p(G|A)+p(T|A) = 1.0
p(A|C)+p(C|C)+p(G|C)+p(T|C) = 1.0
p(A|G)+p(C|G)+p(G|G)+p(T|G) = 1.0
p(A|T)+p(C|T)+p(G|T)+p(T|T) = 1.0

So:
0.6 + 0.1 + 0.1 + 0.2 = 1.0
0.1 + 0.5 + 0.3 + 0.1 = 1.0
0.05+ 0.2 + 0.7 + 0.05 = 1.0
0.4 + 0.05+0.05 + 0.5 = 1.0

Example of the structure of the first order model:

Code:

#    Author: Benjamin Tovar
#    Date: 13/April/2012
#
#    Example of a Markov Chain of zero order (the current nucleotide is
#    totally independent of the previous nucleotide).

#    The multinomial model is:
#    p(A)+p(C)+p(G)+p(T) = 1.0
#    0.4 +0.1 +0.1 +0.4  = 1.0

# Define the DNA alphabet that will be used to put names to the objects
alp <- c("A","C","G","T")
# Create the vector that represents the probability distribution of the model
zeroOrder <- c(0.4,0.1,0.1,0.4)
# Put the name of reference of each base 
names(zeroOrder) <- alp 
# Create a sequence of 1000 bases using this model.
zeroOrderSeq <- sample(alp,1000,rep=T,prob=zeroOrder)

# ***** Study the composition bias of the sequence *****
# We wil use the "seqinr" package.
# For the installation of the package, type:
# install.packages("seqinr")
# Then, load the package:
require("seqinr")

# Count the frequency of each base 
# in the sequence using the "count" function
zeroOrderFreq <- count(zeroOrderSeq,1,alphabet=alp,freq=TRUE)

# Count the frequency of dinucleotides 
# in the sequence using the "count" function
zeroOrderFreqDin <- count(zeroOrderSeq,2,alphabet=alp,freq=TRUE) 
 
# Now, plot the results in the same plot:
layout(1:2)
barplot(zeroOrderFreq,col=1:4,main="Compositional bias of each nucleotide",
        xlab="Base",
        ylab="Base proportion")
barplot(zeroOrderFreqDin,col=rainbow(16),
        main="Compositional bias of each dinucleotide",
        xlab="Base",
        ylab="Base proportion")

# &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

#    Example of a Markov Chain of first order (the current nucleotide only
#    depends on the previous nucleotide).

#    The multinomial model per base is:
#    p(A|A)+p(C|A)+p(G|A)+p(T|A) = 1.0
#    p(A|C)+p(C|C)+p(G|C)+p(T|C) = 1.0
#    p(A|G)+p(C|G)+p(G|G)+p(T|G) = 1.0
#    p(A|T)+p(C|T)+p(G|T)+p(T|T) = 1.0        

#    So:
#    0.6 + 0.1 + 0.1 + 0.2  = 1.0
#    0.1 + 0.5 + 0.3 + 0.1  = 1.0
#    0.05+ 0.2 + 0.7 + 0.05 = 1.0
#    0.4 + 0.05+0.05 + 0.5  = 1.0
    
    
# Create the matrix that will store the probability distribution given 
# a certain nucleotide:
firstOrderMat <- matrix(NA,nr=4,nc=4)
# Put names to the 2 dimensions of the matrix 
colnames(firstOrderMat) <- rownames(firstOrderMat) <- alp 
# Add the probability distribution per base:
firstOrderMat[1,] <- c(0.6,0.1,0.1,0.2)  # Given an A in the 1st position
firstOrderMat[2,] <- c(0.1,0.5,0.3,0.1)  # Given a  C in the 1st position
firstOrderMat[3,] <- c(0.05,0.2,0.7,0.05)# Given a  G in the 1st position
firstOrderMat[4,] <- c(0.4,0.05,0.05,0.5)# Given a  T in the 1st position

# Now we got a matrix
#    > firstOrderMat
#         A    C    G    T
#    A 0.60 0.10 0.10 0.20
#    C 0.10 0.50 0.30 0.10
#    G 0.05 0.20 0.70 0.05
#    T 0.40 0.05 0.05 0.50

# In order to continue, we need an initial probability distribution to know
# which base is the most probable to start up the sequence.
inProb <- c(0.4,0.1,0.1,0.4); names(inProb) <- alp
# So, the sequence will have a 40% to start with an A or a T and 10% with C or G

# Create a function to generate the sequence.
# NOTE: To load the function to the current environment, just copy
# the entire function and paste it inside the R prompt.

generateFirstOrderSeq <- function(lengthSeq,
                                  alphabet,  
                                  initialProb,
                                  firstOrderMatrix){
#    lengthSeq = length of the sequence
#    alphabet = alphabet that compounds the sequence 
#    initialProb   = initial probability distribution
#    firstOrderMatrix = matrix that stores the probability distribution of a 
#                       first order Markov Chain

    # Construct the object that stores the sequence
    outputSeq <- rep(NA,lengthSeq)
    # Which is the first base:
    outputSeq[1]  <- sample(alphabet,1,prob=initialProb) 
    # Let the computer decide:
    for(i in 2:length(outputSeq)){
        prevNuc <- outputSeq[i-1]    
        currentProb <- firstOrderMatrix[prevNuc,]
        outputSeq[i] <- sample(alp,1,prob=currentProb)
    } 
    cat("** DONE: Sequence computation is complete **\n")
    return(outputSeq)
}

# Use the generateFirstOrderSeq function to generate a sequence of 1000 bases 
# long
firstOrderSeq <- generateFirstOrderSeq(1000,alp,inProb,firstOrderMat)

# ***** Study the composition bias of the sequence *****
# We wil use the "seqinr" package.
# For the installation of the package, type:
# install.packages("seqinr")
# Then, load the package:
require("seqinr")

# Count the frequency of each base 
# in the sequence using the "count" function
firstOrderFreq <- count(firstOrderSeq,1,alphabet=alp,freq=TRUE)

# Count the frequency of dinucleotides 
# in the sequence using the "count" function
firstOrderFreqDin <- count(firstOrderSeq,2,alphabet=alp,freq=TRUE) 
 
# Now, plot the results in the same plot:
layout(1:2)
barplot(firstOrderFreq,col=1:4,main="Compositional bias of each nucleotide",
        xlab="Base",
        ylab="Base proportion")
barplot(firstOrderFreqDin,col=rainbow(16),
        main="Compositional bias of each dinucleotide",
        xlab="Base",
        ylab="Base proportion")
        
## Lets plot the 4 plots in one window
    layout(matrix(1:4,nr=2,nc=2))
    # Results from the zero order 
    barplot(zeroOrderFreq,col=1:4,
            main="Compositional bias of each nucleotide\nZero Order Markov Chain",
            xlab="Base",
            ylab="Base proportion")
    barplot(zeroOrderFreqDin,col=rainbow(16),
            main="Compositional bias of each dinucleotide\nZero Order Markov Chain",
            xlab="Base",
            ylab="Base proportion")
    # Results from the first order 
    barplot(firstOrderFreq,col=1:4,
            main="Compositional bias of each nucleotide\nFirst Order Markov Chain",
            xlab="Base",
            ylab="Base proportion")
    barplot(firstOrderFreqDin,col=rainbow(16),
            main="Compositional bias of each dinucleotide\nFirst Order Markov Chain",
            xlab="Base",
            ylab="Base proportion")
        
        
# end.

Plot of the zero order model:

Plot of the first order model

All plots

Benjamin

Generate artificial DNA or protein sequences in R in a single line of code.

2012-04-11T21:10:00.003-05:00

To generate an artificial DNA sequence of "n" bases long with a fixed composition bias in just one line of code, just open your R prompt and type:


seqX <- sample(c("A","C","G","T"),10000,rep=TRUE,prob=c(0.4,0.1,0.1,0.4))

As you see, the alphabet is the four letter alphabet of the DNA (so, if desired, you can replace that alphabet with any other alphabet that you might need, like the amino acids 20 letter code), the next parameter is the length of the desired output sequence. rep=TRUE means that each letter of the alphabet could be repeated and finally, I think that the most useful parameter of the function sample is the prob parameter because it allows you to select the multinomial model (the proportion or "bias" of each base). For example our simulated sequence is 80% AT rich and 20% GC rich given that for the "A" base we got a probability of 0.4, for the "C" base 0.2 and so on.

Now, to check it out that our artificial sequence follows that bias, a simple plot will tell us more than thousand words:

Lets extract the proportion of each base along the sequence using the table and the length function:


freqX <- table(seqX)/length(seqX)

Now, lets plot it doing:


barplot(freqX,col=1:4,main="Compositional bias of seqX",xlab="Base",ylab="Base proportion")

And finally we got this awesome plot.

So, the barplot shows that effectively each base is well represented by the multinomial model and our artificial sequence is loyal to it.

Benjamin

Epic R is Epic <- Beginners command reference card

2012-04-11T20:33:00.001-05:00

This reference card has been written by Tom Short.

Click the image to see the R magic:

Thank you for your support Tom.
Benjamin

Using the Blogger app for Android..

2012-04-06T23:54:00.001-05:00

Well, the application has a simple and elegant UI.

Let's check it out and see how it works :P

Install R 2.15 and further versions in Debian Squeeze

2012-04-06T16:11:00.000-05:00

The last Friday, March 30th, the last stable version of R, the version 2.15.0 was released.

So, to install it in Debian Squeeze, or in another Distro powered by Debian (I actually use CrunchBang Linux), just follow the same instructions described here for the installation of R 2.14.2 by clicking here.

Do not forget to upgrade the packages too.

Benjamin.

Install octave 3.6.1 in Debian Squeeze

2012-04-06T10:55:00.002-05:00

Following the nice instructions found it here <- http://verahill.blogspot.mx/2012/02/debian-testing-wheezy-64-compiling.html , I have successfully installed Octave 3.6.1 in Squeeze (I had to rename some packages names to fit with the Squeeze packages).

FIRST STEP: Install the required libraries and dependencies.



sudo apt-get install gfortran build-essential

sudo apt-get install libqhull-dev libpcre++-dev libblas-dev liblapack-dev libreadline-dev

sudo apt-get install libcurl4-openssl-dev

sudo apt-get install libfltk1.1-dev

sudo apt-get install libgraphicsmagick++-dev

sudo apt-get install libhdf5-serial-dev

sudo apt-get install libqrupdate-dev

sudo apt-get install libsuitesparse-dev

sudo apt-get install  glpk gperf flex bison libfontconfig1-dev

SECOND STEP: Download Octave and extract it.



wget ftp://ftp.gnu.org/gnu/octave/octave-3.6.1.tar.gz

tar -xvf octave-3.6.1.tar.gz

THIRD STEP: Compile it.

 cd octave-3.6.1/

./configure

make -j3

where 3 is the number of cores +1 (for me 2 cores)

FOURTH STEP: Validate the compiled version:

make check

FIFTH STEP: Install Octave



sudo make install

That's all. Thanks to Lindqvist - a blog about Linux and Science. Mostly.
Benjamin

UPDATE: Install R 2.14.2 or R 2.15 in Debian Squeeze

2012-03-29T23:37:00.000-06:00

Since not so long, I upgraded the version of R that is included in the Debian repositories (2.11.1 in Squeeze) and I was very happy with it since 2010. But there was some packages that I could not install and then, I decided to upgrade R despite of having to upgrade every package too (I have more than 1 GB in my home R packages folder).

There's nothing new or tricky about how to upgrade, the awesome guys of CRAN always brings an answer for our needs.

FIRST PART: ADD R BACKPORTS:

First, open a Terminal and open the sources.list file:



$ gksudo gedit /etc/apt/sources.list

Then, add these lines at the bottom of the file (Note, I use the UCLA server, but this can be easily changed taking a look here for the mirrors):

## R BACKPORTS

deb http://cran.stat.ucla.edu/bin/linux/debian squeeze-cran/

#deb-src http://cran.stat.ucla.edu/bin/linux/debian squeeze-cran/

The original set of instructions can be seen here.

SECOND PART: RENAME THE R PACKAGES FOLDER:

There's a folder where R uses to store the packages we download, just rename it to the current version of R.

For example, mine was "2.11" and then I just renamed it to "2.14" (or "2.15", without quotes) and was inside this path:

Before:

/home/benjamin/R/x86_64-pc-linux-gnu-library/2.11

After:

/home/benjamin/R/x86_64-pc-linux-gnu-library/2.15

update.packages(checkBuilt=TRUE, ask=FALSE)

THIRD PART: SECURE APT:

The Debian backports archives on CRAN are signed with the key ID 381BA480, to add them, in a Terminal prompt type:


gpg --keyserver pgp.mit.edu --recv-key 381BA480

gpg -a --export 381BA480 > jranke_cran.asc
sudo  apt-key add jranke_cran.asc

FOURTH PART: UPDATE AND UPGRADE R:

Save the file and you can either enter to Synaptic, update the packages list and then just upgrade the packages or in a terminal type:


sudo apt-get update

sudo apt-get upgrade

And that's all.

Benjamin