HMMs what they are good at....

As biologists most of us feel little bit awkward with mathematics and statistics, and more so with the application of mathematics and statistics in interpreting biological data. Recently we see a meteoric rise in availability of biological sequences. However, interpreting them properly still is a bottleneck.
HMMs are there for a long time now, very valuable in speech recognition, weather prediction etc. But recently, we find a surge in application of HMMs in interpreting biological sequences such as gene prediction, profile searches, domain prediction etc. A beginner however gets awed by the complexity and finds it difficult to understand completely how this amazing method/tool can help interpret sequences. Here is my attempt to simplify things a bit.



What is HMM?

Profile HMMs are based on position specific scoring matrices. Lets take the example of the following diagram where we have multiple sequences aligned.
The following are the characteristics:
1. Columns in a multiple sequence alignment is called as a state.
2. Number of alphabets taken together is called as order of the markov chain e.g; zero order, first order and so on.
3. PSSM(Position specific scoring matrix) is calculated for calculating the match states probabilities.

In a HMM, there is a transition probability and an emission probability and there is a state. Lets look at the simplest HMM heuristic:

Suppose say we have a multiple sequence alignment, as represented in Figure - 1.

Figure 1

[ state -> is the column in a multiple sequence alignment.
Transition probabilities are determined using deletions and insertions.  Emission probabilities at each state are determined by counting the occurrence of each character in each column. ]
 

 So, lets try to generate the Emission and Transition probabilities from this alignment from the figure.

Emission Probability:
                  A                 T                   G                  C
Pos1          1.0               0                    0                   0

Pos2          0                  0                    0.2                0.8
Pos3          0.2               0                    0.4                0
Pos4          0.6               0                    0.2                0

Transition Probability:( There are 3 transition states, match, deletion, insertion. So, there can be 9 + 2 transition states:
                     m->m     m->i     m->d     i->m    i->i    i->d    d->m      d->i    d->d      s->m       m->e
Trans 1         0             0           0            0         0        0         0             0         0           1             0
Trans 2         1             0           0            0         0        0         0             0         0           0             0
Trans 3         0.6          0           0.4         0         0        0         0             0         0           0             0
Trans 4         0.8          0           0.2         0         0        0         0             0         0           0             0
Trans 5         0             0           0            0         0        0         0             0         0           0             1
[In this table m stands for match, i stands for insertion and d for deletion, s for start and e for end]

Constructing the HMM:

There are multiple algorithms that can be used to build and train HMMs. 

i) EM, or Baum-Welsch, algorithm.  

This algorithm uses a random MSA, such as the one constructed by ClustalW.  The number of match states is defined as the number of conserved columns.  This algorithm uses a two stage, iterative process; first, it uses the aligned residues to build a model, and then it realigns the sequences to the model.  These two steps are repeated until the model no longer changes.  

ii)Viterbi algorithm:

Uses dynamic programming to calculate the best estimate of the emission probabilities. 

iii)“Surgery” algorithm:

Uses dynamic adjustments of the HMM length to calculate the best models.

HMMER software comes with a nice user guide. I recommend you read this guide from cover to cover and if you still did not get it, then continue reading my post.

Use of HMM in aligning against another sequence or profile:  
After we create the HMM, we can use it to search for sequences that are similar to the ones used to create the HMM.  So, we use the model to estimate the probability that a sequence belongs to the model.  This probability is often expressed as negative log likelihood score and is dependent on the length of the sequence and the length of the model.  Once we have a score, we need a way to assess its significance, just as we did for pairwise sequence alignments.

HMMs have some distinct advantages.  They are built on a formal probabilistic basis.  They also allow for more sensitive searching than pairwise alignments do.  The probabilistic theory can be used to guide the scoring parameters; the theory also allows us to train an HMM using unaligned sequences if we don’t know or trust the alignment.  There is a consistent theory behind gap and insertion penalties.  Finally, less skill and human intervention is needed to train a good HMM than to train a hand-constructed profile.  Because of this, we can create libraries of hundreds of profile HMMs and apply them on a large scale. 

Its one major drawback is that its probabilistic basis does not capture higher-order correlations between sequence positions.  It assumes that the identity of a particular position is independent of the identity of all of the other positions in the sequence, and we know that not to be true. 

Publicly available program HMMER3 just does this.
 
Components of Hmmer3 Program:
phmmer: Makes a blastp like search.
jackhmmer: Iteratively searches a sequence against a sequence database (something similar to Psi-blast).
hmmbuild: Builds a profile HMM from multiple sequence alignment.
hmmsearch: Searches a profile hmm against a sequence database(Query: HMM, target: AA).
hmmscan: Searches a sequence against a profile hmm database(Query: AA; target: HMM e.g; pfam).
hmmalign: Make a multiple sequence alignment of many sequences to a existing profile.

hmmconvert: Converts profiles into hmm format and vice versa.
hmmemit: Generate sample sequences from profiles.
hmmfetch: Get a profile HMM by name or accession number from the HMM database
hmmpress: Format an hmm database into a binary format to be used by hmmscan.
hmmstat: Show summary statistics for each profile in an hmm database.

In the next post, we will discuss about Applications of Profile search and command line options of Hmmer3 with examples.

Job opportunity as an environment intern!!!

I got this email to my official email account. Thought will share it here..

www.TerraScapes.org

Environmental Organizer Intern

Job Description:

Do you want to incorporate the environment into your everyday life through environmental activism, advocacy, social justice and event planning? Interning for Terra-Scapes Environmental Consulting is a unique opportunity for a student with a passion for the environment.  Engaging the community through educational initiatives about current environmental concerns will be your goal.  The Environmental Organizer will work on several pertinent environmental campaigns ranging from protecting our oceans and fisheries to the prohibition of additional chicken farms on the Eastern Shore.

The Environmental Organizer will have the opportunity to gain experience and forge connections in the environmental and green business communities. The position entails working with the following prominent environmental organizations:  Pew Environment Group, Sierra Club, Chesapeake Bay Foundation, Lynnhaven River Now, Nature Conservancy, and more.

Responsibilities:

· Research, organize, and execute logistics for environmental workshops and events.

· Recruit and communicate with new and existing volunteers.

· Strengthen and build coalitions with new and existing groups, scientists, students, elected officials, and businesses.

· Work with the media on press releases, advisories, letters to the editor, and opinion editorials.

· Schedule and participate in meetings with local, state, and federal legislators regarding environmental policy initiatives.

· Database organization and entry using Google docs and Excel spreadsheets.

· Research and preparation of reports, documents, and memorandums using Microsoft Office.

· Assist in creation of materials for Terra-Scapes and its coalition partners.

· Develop plan to reach goals with each individual campaign.

Qualifications:

· Grasp office-related computer (Mac and PC) and internet applications.

· Have a passion for the environment and the growth of the green economy.

· Ability to prioritize, work independently, and follow through on concurrent projects and goals.

· Superior organizational skills and ability to meet deadlines required.

· Ability to communicate easily with a diverse range of people.

· Demonstrated experience working independently.

· Owns a vehicle and willing to travel throughout VA for specific workshops or conferences.

· Willingness to meet with Congressional members in Washington, D.C.

· Ability to work from your home/dorm.

Location: Primarily in Hampton Roads, but requires some travel throughout Virginia, North Carolina, Maryland, and Washington, D.C.

To apply, please email resume and cover letter to terra@terrascapes.org

Thanks again!

Terra Pascarosa Duff

TerraScapes Environmental Director

VA Pew Environment Group Representative

Chesapeake Bay Group Sierra Club Chair

www.TerraScapes.org

http://www.pewenvironment.org/

p: 757-291-4736 | f: 866-543-6696 | e: terra@terrascapes.org

"You must be the change you wish to see in the world."
— Mahatma Gandhi

Extracting Fastq files from bam/sam files

Earlier I have blogged about converting colorspace files into fastq format, which is pretty easy [here]. However, this method does not go without the risk of getting completely wrong reads if there is an error in the read number in the middle.

There is another approach to it which is reliable - that is to extract the fastq files from the sam/bam alignment files. There are quite a few different ways you can achieve this. One of the easiest is to use samtools to extract matching and unmatching reads separately and merging them together. Something like this:

# Extracting mapping and non-mapping reads from bam files
samtools view -f4 input.bam > /data/bowtieOutput/unmapped.sam
# This is for extracting mapping reads from bam files
samtools view -F4 input.bam > /data/bowtieOutput/mapped.sam

Remember, this method extracts files in sam format, so you may have to get into another step to extract fastq files from sam files

But the only problem with this approach is, you end up getting reads that are mapped to the reference sequence multiple times in multiple numbers. For unmapped reads this problem does not exist.

There is another easier approach is also to use bamtofastq . This is a cool software available here.

Simply use command:

# Extracts aligned reads
./bam2fastq --no-unaligned -o out.fastq input.bam
#Extracts un-aligned reads
./bam2fastq --no-aligned -o out.fastq1 input.bam

This is one of the fastest and easiest methods to extract mapping and unmapping reads from bam files into fastq format. However, this also has similar problem like  the sam extractor - that is, it extracts reads mapped to multiple locations multiple times. In order to get around this the fastest way would be to make a uniq set from the aligned.fastq file. Here are the steps that needs to be followed to get it done faster.

1. Sort your bam file before running bam2fastq on them. This way the output fastq file will have all the redundant reads in a sorted order and hence easier to remove redundant.
2. Run bam2fastq on them as described above.
3. Write a small perl script to fish out the uniq reads from the fastq output file as below:

open FH, $ARGV[0] or die "Can't open file for reading $!\n";

my $prev="";

while(<FH>){


        if(/\@(\d+_\d+_\d+_)/){

                if($prev eq $1){
                        # Skips reading 3 lines
                        for(0..2){
                         <FH>
                        }
                next;
                }
        $prev = $1;
        }

        print $_;
}

My Results:

                         Reads recovered from bam2fastq               uniq reads
                         ----------------------------------------               -------------
expt 1                  23921306                                                    8909517
expt 2                  22685196                                                    8924548
expt 3                  40484226                                                    15248005
expt 4                  13659968                                                    5141253
expt 5                  17185935                                                     6862160
expt 6                  19855362                                                     8978295      
expt 7                  13189432                                                     7067910
expt 8                  14423653                                                     6713690


This is quite consistent with what I wanted to achieve!!                                                      

How to code position fixed content on a browser.

While representing sequence text alignment on a browser, it is very useful to have a static genome reference on the top of the page, that descends as you scroll down. So, you are actually not loosing track of the reference but at the same time you see the read alignments and their position with respect to the genome. The genome references can be long and the browser width is not enough to contain it in a single view. So, there will be a lot of room for horizontal scrolling. In this case you want the reference on the top behave normally as with any horizontal scroll, but vertically it should be position fixed. Here is an example of what I am talking about:


 I have created this application after realizing its importance in dealing with nextgen sequence alignment viewer. Here are the few fixed libraries, you may find incredibly useful, so better install/download these in your local server.

A large part of the application have been downloaded from mootools site. The following are the files I have downloaded from mootools site:

• mootools-yui-compressed.js 

From http://code.google.com/p/joomla-scripts-optimized/downloads/detail?name=mootools-yui-compressed.js&can=2&q=

• DemosMain.php from mootools web site:  (If you don't find it in the web site, the source is copy pasted below)

• Style.php from mootools site (Source code is copy pasted below)

Then bundle up these 3 scripts into a directory say  textStyle.

Now create a Php file and write functions that will print your object of interest Some like this:

 open FH, $fileName or die "Can't open file for reading $! \n";  
 my $count=0;  
 my $start;  
 &printHTMLcode;  
 print "<body>";  
 print "<pre>";  
 while(<FH>){  
     if(/(\d+)\t(\d+)\t(\S+)/){  
         if($count == 0){  
             if($genomeStart >= $2){  
             print "<div id=\"DemoItem2\">";  
                 for(my $i=0; $i<=$genomeStart - $2; $i++){  
                     print "&nbsp;";  
                 }  
             &colorCode($genomeSeq);  
             print "</div><br><br>";  
             $start = $2;  
             }  
             else{  
             print "<div id=\"DemoItem2\">";  
             &colorCode($genomeSeq);  
             $start = $genomeStart+1;  
             print "</div><br><br>";  
             }  
         }  
         for(my $i=0; $i<=($2-$start);$i++){  
             print "&nbsp;";  
         }  
     &colorCode($3);  
     }  
 print "\n";  
 $count++;  
 }  
 print "<script src=\"http://www.eumicrobedb.org/develop/textStyle/mootools-yui-c  
 ompressed.js\" type=\"text/javascript\">   \  
     </script>                  \  
     <script src=\"http://www.eumicrobedb.org/develop/textStyle/demosMain.php  
 \" type=\"text/javascript\">                     \  
     </script>                  \  
 ";  
 print "</pre>";  
 print "</body>";  
 close(FH);  
 }  
 sub printHTMLcode  
 {  
 print "<html>";  
 print "<head>";  
 print "<link rel=\"stylesheet\" href=\"http://www.eumicrobedb.org/develop/textSt  
 yle/style.php\" type=\"text/css\" media=\"screen\">";  
 print "<style type=\"text/css\">  
     #DemoItem2{  
         position: fixed;  
         top   : 10px;  
         left  : 6px;  
         width  : 200px;  
         height : 50px;  
     }  
     #DemoItem2{  
         top: 10px;  
     }  
     </style>";  
 print " <script text=\"text/javascript\">  
         function loadDemo(){  
          var demoTwo = new ScrollSpy({  
                  min: 0,  
                  mode: 'horizontal',  
                  onEnter: function(position,enters) {  
                  },  
                  onLeave: function(position,leaves) {  
                  },  
                  onTick: function(position,state,enters,leaves)  
 {  
                 \$(\"DemoItem2\").style.left = -position.x+ 6+\"  
 px\";  
                  },  
                  container: window  
                 });  
         }  
     </script>  
     </head>  
     <body>  
 \n";  
 }  

I will explain now what all these means.

The file it opens in the while loop is a data file that has sequence and their starting point info. So ignore that part.

There is a function call to printHTML code (marked in red). This function is a style sheet text wrapped. You don't have to wrap it in a function, you could just write it on the top if you like. But it looks messy. Please note that inside printHTML code, I have defined demoItem2 div class that actually does the text/content to scroll down. The main text that scrolls down is within &colorcode function, that actually colors and prints the genome reference that is suppose to be static on the screen.

In the end I am calling the moo tools js files(marked in yellow) that actually lets the content scroll horizontally.

here are the source code for style.php and demomain.php, if you did not find it from the mootools:

Source for style.php

 <?php  
 print "<br>  
 *{<br>  
     margin: 0px;<br>  
     padding: 0px;<br>  
 }<br>  
 <br>  
     body, html{<br>  
         width: 100%;<br>  
         height: 100%;<br>  
 <br>  
     }<br>  
 <br>  
 body{<br>  
     font-family:Gill Sans,sans-serif;<br>  
     background: #2a2a2a;  <br>  
 }<br>  
 <br>  
 #wrap{<br>  
     width: 100%;<br>  
     height: 100%;<br>  
     min-width: 1000px;<br>  
 }<br>  
 <br>  
 #head{<br>  
 height:156px;<br>  
 position:fixed;<br>  
 width:100%;<br>  
 top: 0px;<br>  
 }<br>  
 <br>  
 #main{<br>  
     padding-top:148px;<br>  
 }<br>  
 <br>  
 #logo{<br>  
     background-image:url(images/logo.png);<br>  
     background-position:0 0;<br>  
     height:108px;<br>  
     left:50px;<br>  
     position:absolute;<br>  
     top:26px;<br>  
     width:90px;<br>  
 }<br>  
 <br>  
 #return{<br>  
     color:#FAC500;<br>  
     font-size:13px;<br>  
     position:absolute;<br>  
     right:50px;<br>  
     text-decoration:none;<br>  
     top:122px;<br>  
 } <br>  
 #return:hover{<br>  
     color: #fff;<br>  
 }<br>  
 <br>  
 img, a{<br>  
     outline: none;<br>  
 }<br>  
 <br>  
 h1{<br>  
     color:#FFFFFF;<br>  
     font-family:Helvetica,\"Helvetica Neue\",Arial,sans-serif;<br>  
     font-size:18px;<br>  
     font-weight:normal;<br>  
     letter-spacing:-1px;<br>  
     position:absolute;<br>  
     right:50px;<br>  
     text-shadow:2px 2px 0 #111111;<br>  
     top:95px;<br>  
 }<br>  
 <br>  
 .hr{<br>  
     border-top:1px dashed #555555;<br>  
     height:1px;<br>  
     position:absolute;<br>  
     top:153px;<br>  
     width:100%;<br>  
 }<br>  
 <br>  
 <br>  
 .TabHolder{<br>  
     height: 26px;<br>  
     width: auto;<br>  
     position: fixed;<br>  
     z-index: 100000;<br>  
 }<br>  
 <br>  
 .leftTab, .rightTab{<br>  
     display: inline-block;<br>  
     width: 17px;<br>  
     background-image: url(\"images/leftTabDarker.png\");<br>  
     height: 26px;<br>  
     margin: 0px;<br>  
     float: left;<br>  
 }<br>  
 <br>  
 .rightTab{<br>  
         background-image: url(\"images/rightTabDarker.png\");<br>  
 }<br>  
 <br>  
 .midTab{<br>  
     display: inline-block;<br>  
     width: 17px;<br>  
     background-image: url(\"images/midTabDarker.png\");<br>  
     height: 26px;<br>  
     margin: 0px;<br>  
     min-width: 200px;<br>  
     width: auto;<br>  
     float: left;<br>  
     font-size: 11px;<br>  
     font-family: 'Lucida Sans Unicode', 'Lucida Grande', sans-serif;<br>  
     line-height: 23px;<br>  
     text-shadow: 1px 1px 0px #fff;<br>  
     padding: 0px 2px;<br>  
 }<br>  
 <br>  
 .midTab .midTabcol{<br>  
     text-decoration: none;<br>  
     color: #444;<br>  
     border-right: 1px solid #999;<br>  
     border-left: 1px solid #fff;<br>  
     padding: 0px 10px;<br>  
     display: inline-block;<br>  
     line-height: 12px;<br>  
 }<br>  
 <br>  
 .midTab .midTabcol:first-child{<br>  
     border-left: 0px;<br>  
 }<br>  
 <br>  
 .midTab .midTabcol:last-child, .midTab .lastmidTabcol{<br>  
     border-right: 0px;<br>  
 }<br>  
 <br>  
 .midTab a:hover{<br>  
      color: #111;<br>  
 }<br>  
 <br>  
 .midTab a:focus{<br>  
      position: relative;<br>  
      top: 1px;<br>  
 }<br>  
 <br>  
 .goodBrowsers{<br>  
         display:inline-block;<br>  
         position: relative;<br>  
         top: 3px;<br>  
 }<br>  
 <br>  
 #footerTab{<br>  
     bottom: 0px;<br>  
     right: 100px;<br>  
     bottom: -27px;<br>  
 }<br>  
 <br>  
 <br>  
 ";  
 ?> 

Source code for DemosMain.php

 var Page = {  
     initialize: function() {  
          // nice one David Walsh  
          var ss = new ScrollSpy({  
                  min: 0,  
                  mode: 'vertical',  
                  onEnter: function(position,enters) {  
                  },  
                  onLeave: function(position,leaves) {  
                  },  
                  onTick: function(position,state,enters,leaves)  
 {  
                     $("head").style.top = -position.y+"px";  
                  },  
                  container: window  
                 });  
                 if($("footerTab")){  
                     $("footerTab").set('morph', {  
                     duration: 500,  
                     transition: 'Sine:out'  
                 });  
                     setTimeout(function(){  
                             $("footerTab").morph({  
                     'bottom': ['-27px','0px']  
                     });  
                     },750);  
                 }  
                 loadDemo();  
     }  
  }  
  window.addEvent("domready", Page.initialize);  
 var ScrollSpy = new Class({  
     /* implements */  
     Implements: [Options,Events],  
     /* options */  
     options: {  
         min: 0,  
         mode: 'vertical',  
         max: 0,  
         container: window,  
         onEnter: $empty,  
         onLeave: $empty,  
         onTick: $empty  
     },  
     /* initialization */  
     initialize: function(options) {  
         /* set options */  
         this.setOptions(options);  
         this.container = $(this.options.container);  
         this.enters = this.leaves = 0;  
         this.max = this.options.max;  
         /* fix max */  
         if(this.max == 0)  
         {  
             var ss = this.container.getScrollSize();  
             this.max = this.options.mode == 'vertical' ? ss.y : ss.x  
 ;  
         }  
         /* make it happen */  
         this.addListener();  
     },  
     /* a method that does whatever you want */  
     addListener: function() {  
         /* state trackers */  
         this.inside = false;  
         this.container.addEvent('scroll',function() {  
             /* if it has reached the level */  
             var position = this.container.getScroll();  
             var xy = this.options.mode == 'vertical' ? position.y :  
 position.x;  
             /* if we reach the minimum and are still below the max..  
 . */  
             if(xy >= this.options.min && xy <= this.max) {  
                     /* trigger Enter event if necessary */  
                     if(!this.inside) {  
                         /* record as inside */  
                         this.inside = true;  
                         this.enters++;  
                         /* fire enter event */  
                         this.fireEvent('enter',[position  
 ,this.enters]);  
                     }  
                     /* trigger the "tick", always */  
                     this.fireEvent('tick',[position,this.ins  
 ide,this.enters,this.leaves]);  
             }  
             else {  
                 /* trigger leave */  
                 if(this.inside)  
                 {  
                     this.inside = false;  
                     this.leaves++;  
                     this.fireEvent('leave',[position,this.le  
 aves]);  
                 }  
             }  
         }.bind(this));  
     }  
 });