View source for Canavan Task 2 - Sequence alignments

==Protocol==

Further information can be found in the [[CD_task2_protocol|protocol]].


== <div id="go_terms">GO Term Enrichment</div> ==

In the following we are performing different sequence searches with the protein Aspartoacylase (UniProt ID: P45381). In order to validate the found hits, we are looking for common GO classifications of the hits with the query sequence.

For our protein Aspartoacylase there are 17 annotated GO terms (using [http://www.ebi.ac.uk/QuickGO/ EMBLs QuickGO]):

<figtable id="aspa_go_terms">
<table  cellspacing=0 align="center" cellpadding=5>
<caption align="center"><b><xr nolink id="aspa_go_terms"/></b> <br> The 17 annotated GO terms for Aspartoacylase(P45381) </caption>

<tr><td style="border-right:solid;border-bottom:solid" align="left"></td><td style="border-bottom:solid;" align="left"><b>GO ID</b></td><td style="border-bottom:solid;" align="left"><b>GO Name</b></td></tr>

<tr><td style="border-right:solid;" align="left"><b>Cellular Compartment</b>
</td><td> GO:0005634 </td>
<td>nucleus (3X)</td>
<tr><td style="border-bottom:solid;border-right:solid;" align="left"></td>
<td style="border-bottom:solid;" align="left">GO:0005634 </td>
<td style="border-bottom:solid;" align="left">cytoplasm (4X)</td></tr>

<tr><td style="border-right:solid;" align="left"><b>Biological Process</b></td><td>GO:0006533</td><td>	aspartate catabolic process</td></tr>
<tr><td style="border-right:solid;" align="left"></td><td>GO:0008152</td><td>	metabolic process</td></tr>
<tr><td style="border-right:solid;" align="left"></td><td>GO:0022010</td><td>	central nervous system myelination</td></tr>
<tr><td style="border-bottom:solid;border-right:solid;" align="left"></td><td style="border-bottom:solid;" align="left">GO:0048714</td><td style="border-bottom:solid;" align="left">	positive regulation of oligodendrocyte differentiation</td></tr>

<tr><td style="border-right:solid;" align="left"><b>Molecular Function</b></td><td>GO:0046872	</td><td>metal ion binding</td></tr>
<tr><td style="border-right:solid;" align="left"></td><td>GO:0004046</td><td>	aminoacylase activity</td></tr>
<tr><td style="border-right:solid;" align="left"></td><td>GO:0016787</td><td>	hydrolase activity</td></tr>
<tr><td style="border-right:solid;" align="left"></td><td>GO:0016788</td><td>	hydrolase activity, acting on ester bonds</td></tr>
<tr><td style="border-right:solid;" align="left"></td><td>GO:0016811</td><td>	hydrolase activity, acting on carbon-nitrogen (but not peptide) bonds, in linear amides</td></tr>
<tr><td style="border-right:solid;" align="left"></td><td>GO:0019807</td><td>	aspartoacylase activity</td></tr>
</table>
</figtable>

== Pairwise Sequence Search==

=== BLASTP ===

<table>
<tr>
<td width="70%">
Even when using the non-restricitve default EValue cutoff, we only find 196 hits. Using a stricter cutoff, we find 94 significant hits. Most of the resulting proteins are Aspartoacylases of other species. Most of the results with EValue > e-15 are Succinylglutamate Desuccinylases, which are in the same protein  family ([http://pfam.sanger.ac.uk/family/PF04952|Succinylglutamate Desuccinylase / Aspartoacylase family]) and catalyze a reaction similar to Aspartoacylase.
</td>

<td>
<figtable id="aspa_blastp">
<table  cellspacing=0 align="center" cellpadding=5>
<caption align="center"><b><xr nolink id="aspa_blastp"/></b> <br> Results from a BlastP search with Aspartoacylse using different EValues.</caption>
<tr>
<td style="border-bottom:solid;border-right:solid;" align="left" width="160px"> <b>E-Value</b></td>
<td style="border-bottom:solid;border-right:solid;" align="center"> default E-Value (10)</td>

</tr>
<tr>
<td style="border-right:solid;" align="left"><b>results (EVal: 10)</b></td><td style="border-right:solid;" align="center">196</td>
</tr>
<tr>
<td style="border-right:solid;" align="left"><b>results (EVal: 2e-3)</b></td><td style="border-right:solid;" align="center">104</td>
</tr>
<tr>
<td style="border-right:solid;" align="left"><b>best E-Value</b></td><td style="border-right:solid;" align="center">1e-155</td>
</tr>

</table>
</figtable>

</td></tr></table>

===PSIBLAST===

We run PSIBlast with four different parameter combinations: 
* 2 iterations (j=2), default E-Value cutoff for inclusion of .. (h=0.002) 
* 2 iterations (j=2), strict E-Value cutoff for inclusion of .. (h=10e-10) 
* 10 iterations (j=10), default E-Value cutoff for inclusion of .. (h=0.002) 
* 10 iterations (j=10), strict E-Value cutoff for inclusion of .. (h=10e-10) 

As a cutoff value for significance of hits we chose 2e-10.

The Psiblast run with default parameters (2 iterations, EValue 0.002) results in an amount of 597 hits. The most significant hit has an E-Value of e-142. 
When restricting the search to a lower E-Value of 10e-10, we still get 502 hits. 

only get 93 hits. This is similar to the result from the BlastP search, only that the worst hit still has a significant EValue of 2e-29. Hits with the good EValues are mostly Aspartoacylases, while hits with worse EValues contain more and more Succinylglutamate Desuccinylases.

Increasing the number of iterations obviously results in many more hits. Using different EValue cutoffs has little influence on the found results though. After 8 rounds, the search converged using the default EValue, while it run all 10 rounds using the more restrictive EValue cutoff. Interestingly, the best hits are less significant (~e-70) than for the psiblast search with only 2 iterations (~e-140). Also, we feel that the most significant results include more Succinylglutamate Desuccinylases than Aspartoacylases.


<figtable id="aspa_psiblast">
<table  cellspacing=0 align="center" cellpadding=5>
<caption align="center"><b><xr nolink id="aspa_psiblast"/></b> <br> Results from a psiblast search with Aspartoacylse using different EValues and iterations.</caption>

<tr>
<td style="border-bottom:solid;border-right:solid;" align="left"> <b>Parameters</b></td> 
<td style="border-bottom:solid;border-right:solid;" align="center"><b> it2, def E-Value (h=2e-3)</b></td> 
<td style="border-bottom:solid;border-right:solid;" align="center"><b> it2, E-Value h=10e-10</b></td> 
<td style="border-bottom:solid;border-right:solid;" align="center"><b> it10 def E-Value (h=2e-3)</b></td>
<td style="border-bottom:solid;border-right:solid;" align="center"><b> it10 E-Value h=10e-10</b></td>
</tr>


<tr><td style="border-right:solid;" align="left"><b>time</b></td>       <td style="border-right:solid;" align="center">~2m30</td>                    <td style="border-right:solid;" align="center">~2m30</td>              <td style="border-right:solid;" align="center">~10m</td>                    <td style="border-right:solid;" align="center">~10m</td></tr>

<tr><td style="border-right:solid;" align="left"><b>best E-Value</b></td><td style="border-right:solid;" align="center">1e-142</td>                  <td style="border-right:solid;" align="center">1e-145</td>            <td style="border-right:solid;" align="center">5e-70</td>               <td style="border-right:solid;" align="center">9e-64</td></tr>

<tr><td style="border-right:solid;" align="left"><b>results (EVal: 10)</b></td>    <td style="border-right:solid;" align="center">915</td>                      <td style="border-right:solid;" align="center">835</td>                <td style="border-right:solid;" align="center">500</td>                     <td style="border-right:solid;" align="center">500</td></tr>

<tr><td style="border-right:solid;" align="left"><b>results (EVal: 2e-3)</b></td>    <td style="border-right:solid;" align="center">597</td>                      <td style="border-right:solid;" align="center">502</td>                <td style="border-right:solid;" align="center">500</td>                     <td style="border-right:solid;" align="center">500</td></tr>


</table>
</figtable>

===HHBLITS===


<tr><td><b>comment</b></td><td width=300>mixed results with Aspartoacylases and Succi</td><td>very varying results: Aspartoacylasen, Succinylasen, Zinc Proteins</td></tr>

<figtable id="aspa_hhblits">
<table  cellspacing=0 align="center" cellpadding=5>
<caption align="center"><b><xr nolink id="aspa_psiblast"/></b> <br> Results from a HHblits search with Aspartoacylse using different EValues and iterations.</caption>

<tr><td style="border-bottom:solid;border-right:solid;" align="left"><b>Parameters</b></td>
<td style="border-bottom:solid;border-right:solid;" align="center">it 2</td>
<td style="border-bottom:solid;border-right:solid;" align="center"> it 8</td></tr>

<tr><td style="border-right:solid;" align="left"><b>time</b></td>
<td  style="border-right:solid;" align="center">2m50</td>
<td  style="border-right:solid;" align="center">~6m</td></tr>

<tr><td style="border-right:solid;" align="left"><b>results</b></td>
<td  style="border-right:solid;" align="center">274</td>
<td  style="border-right:solid;" align="center">500</td></tr>

<tr><td style="border-right:solid;" align="left"><b>best E-Value</b></td>
<td  style="border-right:solid;" align="center">2e-110</td>
<td  style="border-right:solid;" align="center">2.9e-68</td></tr>

<tr><td style="border-right:solid;" align="left"><b>worst E-Value</b></td>
<td  style="border-right:solid;" align="center">0.0011</td>
<td  style="border-right:solid;" align="center">9.5e-09</td></tr>

</table>
</figtable>

==Validation and Comparison==

Along with the expactations one can find more hits with Psi-Blast than with a simple Blast search. 

In general, one can distinguish between two kinds of proteins, that frequently are identified by the sequence searches:
* Aspartoacylases 
* Succinylglutamate Desuccinylases 



===BlastP===

<figtable id="blastp_comp">
<table>
<tr><td width="60%">
A simple blast search yields only about 90 significant hits if one considers a threshold of 10e-10 as a significance cutoff.
As one can see in Figure <xr id="blastp_comp"/>, the restriction of the E-Value results in less hits with a low sequence similarity.


192 out of the 196 hits have annotated GO terms. As one can see in <xr id="go_blastp_def" /> all hits share many GO terms with Aspartacylase. The most common GO terms shared with Aspartoacylase are "metabolic process" (185x) and "hydrolase activity, acting on ester bonds" (184x). The GO term "Zinc binding" is not officially associated with Aspartoacylase, yet we know that is does bind zinc. As already mentioned, some of the found proteins have a succinylglutamate desuccinylase activity. They belong to the same family as Aspartoacylase and their occurence is not surpising.

For all 94 hits found with an E-Value cutoff of 10e-10, there are annotated GO terms. Furthermore all founds hits share the GO term "hydrolase activity, acting on ester bonds" and "metabolic process". Also, as one can see in <xr id="go_blastp_10e10" /> all hits share the most GO terms with Aspartacylase. Again, "Zinc binding" could also be associated with Aspartoacylase.  Therefore, all GO terms that are found more than 5 times, are also associated with Aspartoacylase. The results are more accurate concerning shared GO terms with Aspartoacylase. This is what one would expect when restricting the EValue for finding the closer related proteins.

</td>
<td>
<figure id="blastp_comp">
[[File:CD_blast_seqid_distri.png|thumb|350px|center|<b><xr nolink id="blastp_comp"/></b><br>Comparison of distribution of Sequence Identity between the two BlastP runs ]]</figure>
</td></tr></table>
</figtable>



<figtable id="go_blastp">
<table cellspacing=0 align="center" cellpadding=5>
<caption align="center"><b><xr nolink id="aspa_blastp"/></b> <br> Results from a BlastP search with Aspartoacylse using different EValues.</caption>
<tr>
<td style="border-bottom:solid;border-right:solid;" align="left" width="100px"> <b>E-Value</b></td>
<td style="border-bottom:solid;border-right:solid;" align="center"> default (10)</td>
<td style="border-bottom:solid;border-right:solid;" align="center"> 10e-10</td>
</tr>

<tr>
<td style="border-right:solid;" align="left" valign="top"><b>GO terms</b></td>
<td style="border-right:solid;" align="left" valign="top">
<figure id="go_blastp_def">
[[File:cd_go_blastp_def.png|thumb|center|450px|<b><xr nolink id="go_blastp_def"/></b><br>Go Term enrichment for the hits found with BlastP using the default E-Value cutoff of 10. All Go terms with occurence more than once are shown. Go Terms that are identical with the Go annotation for Aspartoacylase (P45381) are colored red.]]</figure>
</td>
<td valign="top" align="left">
<figure id="go_blastp_10e10">
[[File:cd_go_blastp_10e10.png|thumb|center|350px|<b><xr nolink id="go_blastp_10e10"/></b><br>Go Term enrichment for the hits found with BlastP using an EValue cutoff of 10e-10. All Go terms with occurence more than once are shown. Go Terms that are identical with the Go annotation for Aspartoacylase (P45381) are colored red.]]</figure>
</td>
</tr>
</table>
</figtable>

===Psi Blast===

Increasing the amount of iterations performed in a PSI-Blast search, obviously increases the running time. One can see, that the best ranked hits of the runs with 10 iterations have lower E-Values than the best hits of the runs with less iterations. Yet, the result includes a larger amount of significant hits with higher E-Values. This means, increasing the iterations finds further distantly related sequences, which is the expected outcome. 
This outcome is also represented in the distribution of sequence identities. As one can see in figure ??, running PSI-Blast with 10 iterations results in hits with a lower sequence identity to our query sequence than the hits from the run with 2 iterations.



When restricting the E-Value Cutoff for the profile built-up, we found that more hits are classified as Aspartoacylases than as Succinylglutamate Desuccinylases. The running time, as well as the E-Values of the resulting hits did not change significantly. The majority of the results from the runs with only two iterations, has moderate sequence identities with a broad distribution between 10% and 50%. In contrast, the results from the run with 10 iterations split up into two groups of hits which form cluster at about 15% and 35% sequence identity. This difference is also represented in the E_Value distribution. The runs with 10 iterations result in Hits with moderate E_Values between -200 and -40 log(E_Values). The runs with 10 iterations in contrast result in many low significant hits (log(E_Value > -20)) and a variety of high significant hits. 

<center>
<table>
<tr>
<td>
<figure id="PSI_10e10_seqd">
[[File:CD_PSI_10e10_seqd_distri.png|thumb|250px|center|<b><xr nolink id="PSI_10e10_seqd"/></b><br> Distribution of Sequence Id between Psi-Blast runs with 2 iterations vs 10 iterations (using E-Value 10e-10)]]</figure>
</td>
<td>
<figure id="PSI_2it_seqid">
[[File:CD_PSI_2it_seqid_distri.png |thumb|260px|center|<b><xr nolink id="PSI_2it_seqid"/></b><br> Distribution of Sequence Id for Psi-Blast runs with 2 iterations with different E-Values (def E-Value vs E-Value of10e-10)]]
</figure>
</td>
<td>
<figure id="PSI_10it_seqid">
[[File:CD_PSI_10it_seqid_distri.png |thumb|250px|center|<b><xr nolink id="PSI_10it_seqid"/></b><br> Distribution of Sequence Id for Psi-Blast runs with 10 iterations with different E-Values (def E-Value vs E-Value of10e-10)]]
</figure>
</td>
<td>
<figure id="eval_distri_psiblast">
[[File:cd_eval_distri_psiblast.png|thumb|310px|center|<b><xr nolink id="eval_distri_psiblast"/></b><br> Distribution of logarithmic E_Values for the four different PSIBlast runs]]
</figure>
</td>
</tr>
</table>
</center>


For the run with 2 iterations and the default cutoff value of 0.002, we received 915 hits. We considered 597 hits as significant (E-Value cutoff 2e-3). 586 of these significant proteins have GO terms annotated.
For the run with 2 iterations and a EValue cutoff of 10e-10, we received 835 hits out of which we considered 502 proteins as significant. 496 proteins have GO terms annotated.

===HHBlits===

Running HHBlits with 2 iterations yields a small amount of hits (270) with very low (2e-110) and very high (0.0011) E-Values. To increase the amount of hits, we repeated the HHBlits search with the maximum amount of 8 iterations which resulted in a broader output with more Hits with lower averaged E-Values (compare figure ??).  
Regarding the Sequence Identity distribution, running HHBlits with 8 iterations results in more distant related Hits (see Figure ??).

<table>
<tr>
<td>
[[File:CD_HHBlits_seqid_distri.png|thumb|300px|<b>Figure ??</b><br> Sequence identity distributions of HHBlits run with 2 and with 8 iterations.]]
</td>
<td>[[File:cd_eval_distri_HHblits.png|thumb|300px|<b>Figure ??</b><br> logarithmic E_Value distributions of HHBlits run with 2 and with 8 iterations.]]</td>
</tr>
</table>

===Overlap===
As one can see in Figure ??, roughly 40 percent of the resulting hits are unique to each method. From our considerations, about 25 percent of the hits are significant hits, that could be further investigated (overlap of 50 percent). 
[[File:seq_overlap.png|thumb|300px|center|<b>Figure ??</b> <br> Distribution of overlapping Hits for the eight different used Sequence Searches.]]


<table>
<tr>
<td>[[File:defaults.png|200px|thumb|Default E-Values: as could be expected, the normal BLAST search is mostly contained in the PsiBLAST search with two iterations. HHBlits found a large number of different hits, with only 48 out of 274 common hits in common with the BLAST searches.]]</td>
<td>[[File:Psi2_psi10_blast_defs.png|200px|thumb|Taking PsiBLAST with 10 iterations into account brings in a large number of common sequences among the three searches (110), which could be interesting since there seems to be high conversation among them.]]</td>
<td>[[File:blast_psi2_e10_hhblits.png|200px|thumb|Strict E-Values for PsiBLAST and default E-Value for HHBlits with 2 iterations: The number of common hits among all three is now substantially lower, while PsiBLAST with two and ten iterations share a great number of their hits. ]]</td>
<td>[[File:Psi2_psi10_e10_hhblits_n10.png|200px|thumb|Increasing the number of HHBlits-iteration yields more hits for HHBlits, but does not increase the number of common hits with PSI-Blast in 2 or 10 iterations. However, 10 sequences are common and could be interesting for further investigation.]]</td>
</tr></table>




===Further Evaluation===
We tried to further validate the sequence search hits via structural similarity. Unfortunately none of the resulting Hits was a PDB Hit. Furthermore we tried to map the sequence identifiers against the UniProtKB/Swiss-Prot PDB cross-references (http://www.uniprot.org/docs/pdbtosp.txt). Again, this mapping yielded no results, which is why we cannot include any structural information for our ongoing research. When inspecting the annotation for the sequence hits, we already found, that the majority of the hits codes for Aspartoacylases or respectively the highly related protein Succinylglutamate Desuccinylases. Since there already exists a crystal structure of the human Aspartoacylase, it is only reasonable that one will not find other structures for this class of proteins. Additionally, a huge amount of hits codes for not yet characterized proteins, which also will hardly be an interesting target for crystallization.

== Multiple Sequence Alignments ==

For generating our dataset for the MSA we clustered all Hits into Sequence Identity groups:
* >90%: 1 
*60-89%: 59
*40-59%: 197
*20-39%: 1141

Since we only got one hit with an sequence Identity >90% we decided to group out hits as follows:
three groups of sequences with eight members each:
*60-99%
*40-59%
*20-39%

We chose those hits from the respective groups, that have been found by at least 4 methods (overlap of 50%).


 id       eVal  identity
 
 # 60-99% sequence identity
 Q8BZC2	1.7e-25	90
 E1BVP5	e-140	72
 H2RVG4	e-141	63
 G3VM93	e-105	72
 F6ZFQ0	e-139	78
 F8WFU8	e-145	86
 Q28C61	e-132	68
 H2M5L4	e-133	64
 
 # 40-59% sequence identity
 G5BTW1	e-133	43
 G6FRX8	e-103	39
 F7NV91	e-112	39
 G1Q6P7	e-120	42
 H0WH68	e-135	44
 F2PFG6	e-119	40
 H2MX25	5e-81	40
 Q1Z2X2	e-115	38
 
 # 20-39% sequence identity
 Q2F9Q7	e-109	31
 Q8YQC1	e-117	41
 E1SMZ8	e-108	39
 D7E1T3	e-110	36
 A5GQV1	7e-92	33
 E8LP14	e-107	31
 F9TUZ3	e-106	30
 A6VUE4	e-101	35


===General Results===

All in all the three Alignment methods yield comparable results. One can identify several conserved regions. Especially the two groups with sequence identities <60% show very similar MSAs.

There are three strongly conserved motivs located in the first half of the sequences:
*GGTHGNE
*DLNR
*DLHNT

For the second half of the sequence alignments there is no clear concensus about reserved motifs, but several residues are strongly conserved and may be of  functional or structural importance.

In the alignment of the >60% group the first two motifs are not colored in the alignment. This is due to two very short sequences which produce gaps in the alignment and thus lower the consensus.

===ClustalW===
command: <code> clustalw -align -infile=./db_over60.fa -outfile=./clustalw_msa_60.aln </code>



<table>
<tr>
<td>[[File:clustal_20.png|thumb|300px|Jalview Representation of the ClustalW Alignment with the dataset with 20-39% sequence identity. Colored are conserved residues (>65%).]]</td>

<td>[[File:clustal_40.png|thumb|300px|Jalview Representation of the ClustalW Alignment with the dataset with 40-59% sequence identity. Colored are conserved residues (>65%).]]</td>

<td>[[File:clustal_60.png|thumb|300px|Jalview Representation of the ClustalW Alignment with the dataset with 60-100% sequence identity.Colored are conserved residues (>65%).]]</td>
</tr>
</table>

===TCoffee===
<table>
<tr>
<td>[[File:tcoffee_20.png|thumb|300px|Jalview Representation of the T-Coffee Alignment with the dataset with 20-39% sequence identity. Colored are conserved residues (>65%).]]</td>

<td>[[File:tcoffee_40.png|thumb|300px|Jalview Representation of the T-Coffee Alignment with the dataset with 40-59% sequence identity. Colored are conserved residues (>65%).]]</td>

<td>[[File:tcoffee_60.png|thumb|300px|Jalview Representation of the T-Coffee Alignment with the dataset with 60-100% sequence identity.Colored are conserved residues (>65%).]]</td>
</tr>
</table>

===Muscle===

<table>
<tr>
<td>[[File:muscle_20.png|thumb|300px|Jalview Representation of the Muscle Alignment with the dataset with 20-39% sequence identity. Colored are conserved residues (>65%).]]</td>

<td>[[File:muscle_40.png|thumb|300px|Jalview Representation of the Muscle Alignment with the dataset with 40-59% sequence identity. Colored are conserved residues (>65%).]]</td>

<td>[[File:muscle_60.png|thumb|300px|Jalview Representation of the Muscle Alignment with the dataset with 60-100% sequence identity.Colored are conserved residues (>65%).]]</td>
</tr>
</table>

===Concerning the wildtype human Aspartoacylase===

The three identified motifs can also be foung in the wildtype protein (compare coloring of sequence on top of page). We colored the respective residues in the structure. They all are positioned in the same region of the protein and thus might implicate an important functional region.

<table>
<tr>
<td>[[File:DLHNT.png|thumb|300px|Human Aspartoacylse (PDB 2O53) with the motif residues DLHNT]]</td>

<td>[[File:DLNR.png|thumb|300px|Human Aspartoacylse (PDB 2O53) with the motif residues DLNR]]</td>

<td>[[File:GGLNR.png|thumb|300px|Human Aspartoacylse (PDB 2O53) with the motif residues GGLNR.]]</td>
</tr>
</table>